types.hh

/*
 * types.hh
 *
 * License: Artistic License, see file LICENSE.TXT or 
 *          https://opensource.org/licenses/artistic-license-1.0
 */
 
#ifndef _TYPES_HH
#define _TYPES_HH
 
// project includes
#include "lldouble.hh"
 
 
#ifdef TESTING
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/impl/basic_text_oarchive.ipp>
#include <boost/archive/impl/text_oarchive_impl.ipp>
#include <boost/serialization/list.hpp>
#include <boost/serialization/vector.hpp>
 
/*
 * author Stefanie Koenig (included here and adapted by Giovanna Migliorelli on 14.05.2020)
 * append slash to the end of a directory name 
 */
void expandDir(string& filename);
#endif
 
// standard C/C++ includes
#include <string>
#include <exception>
#include <sstream>
#include <cstring>
#include <cstdlib>
#include <vector>
#include <limits>
#include <algorithm> // for std:sort
#include <map>
#ifdef COMPGENEPRED
#include <unordered_map>
#endif
 
using namespace std;
 
enum Strand {STRAND_UNKNOWN=-1, plusstrand, minusstrand, bothstrands};
char strandChar (Strand s);
ostream& operator<< (ostream& strm, const Strand s);
 
#define ASS_MIDDLE      2
#define DSS_MIDDLE      2
#define STARTCODON_LEN  3
#define STOPCODON_LEN   3
#define OCHRECODON      "taa"
#define OPALCODON       "tga"
#define AMBERCODON      "tag"
#define RCOCHRECODON    "tta"
#define RCOPALCODON     "tca"
#define RCAMBERCODON    "cta"
#define TRUNC_LEFT      1
#define TRUNC_RIGHT     2 
 
#define SPECIES_SUBDIR "species/"
#define MODEL_SUBDIR "model/"
#define EXTRINSIC_SUBDIR "extrinsic/"
 
#define VERSION "3.5.0"
 
#define PREAMBLE "# This output was generated with AUGUSTUS (version " VERSION ").\n\
# AUGUSTUS is a gene prediction tool written by M. Stanke (mario.stanke@uni-greifswald.de),\n\
# O. Keller, S. König, L. Gerischer, L. Romoth, Katharina Hoff, Henry Mehlan and Daniel Honsel.\n\
# Please cite: Mario Stanke, Mark Diekhans, Robert Baertsch, David Haussler (2008),\n\
# Using native and syntenically mapped cDNA alignments to improve de novo gene finding\n\
# Bioinformatics 24: 637-644, doi 10.1093/bioinformatics/btn013"
 
#define GREETING "AUGUSTUS (" VERSION ") is a gene prediction tool.\n\
Sources and documentation at https://github.com/Gaius-Augustus/Augustus"
 
#define SPECIES_LIST "usage:\n\
augustus [parameters] --species=SPECIES queryfilename\n\
\n\
where SPECIES is one of the following identifiers\n\
\n\
identifier                               | species\n\
-----------------------------------------|----------------------\n\
pea_aphid                                | Acyrthosiphon pisum\n\
aedes                                    | Aedes aegypti\n\
amphimedon                               | Amphimedon queenslandica\n\
ancylostoma_ceylanicum                   | Ancylostoma ceylanicum\n\
adorsata                                 | Apis dorsata\n\
honeybee1                                | Apis mellifera\n\
arabidopsis                              | Arabidopsis thaliana\n\
aspergillus_fumigatus                    | Aspergillus fumigatus\n\
aspergillus_nidulans                     | Aspergillus nidulans\n\
(anidulans)                              | Aspergillus nidulans\n\
aspergillus_oryzae                       | Aspergillus oryzae\n\
aspergillus_terreus                      | Aspergillus terreus\n\
bombus_impatiens1                        | Bombus impatiens\n\
bombus_terrestris2                       | Bombus terrestris\n\
botrytis_cinerea                         | Botrytis cinerea\n\
brugia                                   | Brugia malayi\n\
b_pseudomallei                           | Burkholderia pseudomallei\n\
caenorhabditis                           | Caenorhabditis elegans\n\
(c_elegans_trsk)                         | Caenorhabditis elegans\n\
(elegans)                                | Caenorhabditis elegans\n\
elephant_shark                           | Callorhinchus milii\n\
camponotus_floridanus                    | Camponotus floridanus\n\
candida_albicans                         | Candida albicans\n\
candida_guilliermondii                   | Candida guilliermondii\n\
candida_tropicalis                       | Candida tropicalis\n\
chaetomium_globosum                      | Chaetomium globosum\n\
chiloscyllium                            | Chiloscyllium punctatum\n\
chlamy2011                               | Chlamydomonas reinhardtii\n\
(chlamydomonas)                          | Chlamydomonas reinhardtii\n\
chlorella                                | Chlorella sp.\n\
ciona                                    | Ciona intestinalis\n\
coccidioides_immitis                     | Coccidioides immitis\n\
coccidioides_immitis                     | Coccidioides immitis\n\
Conidiobolus_coronatus                   | Conidiobolus coronatus\n\
coprinus_cinereus                        | Coprinus cinereus\n\
coprinus_cinereus                        | Coprinus cinereus\n\
(coprinus)                               | Coprinus cinereus\n\
coprinus                                 | Coprinus cinereus\n\
cryptococcus_neoformans_gattii           | Cryptococcus neoformans gattii\n\
cryptococcus_neoformans_neoformans_B     | Cryptococcus neoformans neoformans\n\
(cryptococcus)                           | Cryptococcus neoformans\n\
culex                                    | Culex pipiens\n\
zebrafish                                | Danio rerio\n\
debaryomyces_hansenii                    | Debaryomyces hansenii\n\
fly                                      | Drosophila melanogaster\n\
(fly_exp)                                | Drosophila melanogaster\n\
encephalitozoon_cuniculi_GB              | Encephalitozoon cuniculi\n\
eremothecium_gossypii                    | Eremothecium gossypii\n\
E_coli_K12                               | Escherichia coli K12\n\
fusarium_graminearum                     | Fusarium graminearum\n\
(fusarium)                               | Fusarium graminearum\n\
galdieria                                | Galdieria sulphuraria\n\
chicken                                  | Gallus gallus domesticus\n\
heliconius_melpomene1                    | Heliconius melpomene\n\
histoplasma_capsulatum                   | Histoplasma capsulatum\n\
(histoplasma)                            | Histoplasma capsulatum\n\
human                                    | Homo sapiens\n\
kluyveromyces_lactis                     | Kluyveromyces lactis\n\
laccaria_bicolor                         | Laccaria bicolor\n\
leishmania_tarentolae                    | Leishmania tarentolae\n\
japaneselamprey                          | Lethenteron camtschaticum\n\
lodderomyces_elongisporus                | Lodderomyces elongisporus\n\
magnaporthe_grisea                       | Magnaporthe grisea\n\
mnemiopsis_leidyi                        | Mnemiopsis leidyi\n\
nasonia                                  | Nasonia vitripennis\n\
nematostella_vectensis                   | Nematostella vectensis\n\
neurospora_crassa                        | Neurospora crassa\n\
(neurospora)                             | Neurospora crassa\n\
coyote_tobacco                           | Nicotiana attenuata\n\
rice                                     | Oryza sativa\n\
parasteatoda                             | Parasteatoda sp.\n\
sealamprey                               | Petromyzon marinus\n\
phanerochaete_chrysosporium              | Phanerochaete chrysosporium\n\
(pchrysosporium)                         | Phanerochaete chrysosporium\n\
pichia_stipitis                          | Pichia stipitis\n\
pisaster                                 | Pisaster ochraceus\n\
pfalciparum                              | Plasmodium falciparum\n\
pneumocystis                             | Pneumocystis jirovecii\n\
rhincodon                                | Rhincodon typus\n\
rhizopus_oryzae                          | Rhizopus oryzae\n\
rhodnius                                 | Rhodnius prolixus\n\
saccharomyces_cerevisiae_rm11-1a_1       | Saccharomyces cerevisiae\n\
saccharomyces_cerevisiae_S288C           | Saccharomyces cerevisiae\n\
(saccharomyces)                          | Saccharomyces cerevisiae\n\
schistosoma                              | Schistosoma mansoni\n\
(schistosoma2)                           | Schistosoma mansoni\n\
schizosaccharomyces_pombe                | Schizosaccharomyces pombe\n\
scyliorhinus                             | Scyliorhinus torazame\n\
tomato                                   | Solanum lycopersicum\n\
s_aureus                                 | Staphylococcus aureus\n\
s_pneumoniae                             | Streptococcus pneumoniae\n\
strongylocentrotus_purpuratus            | Strongylocentrotus purpuratus\n\
sulfolobus_solfataricus                  | Sulfolobus solfataricus\n\
tetrahymena                              | Tetrahymena thermophila\n\
cacao                                    | Theobroma cacao\n\
thermoanaerobacter_tengcongensis         | Thermoanaerobacter tengcongensis\n\
toxoplasma                               | Toxoplasma gondii\n\
tribolium2012                            | Tribolium castaneum\n\
trichinella                              | Trichinella sp.\n\
wheat                                    | Triticum sp.\n\
ustilago_maydis                          | Ustilago maydis\n\
(ustilago)                               | Ustilago maydis\n\
verticillium_albo_atrum1                 | Verticillium albo atrum\n\
verticillium_longisporum1                | Verticillium longisporum\n\
volvox                                   | Volvox sp.\n\
Xiphophorus_maculatus                    | Xiphophorus maculatus\n\
yarrowia_lipolytica                      | Yarrowia lipolytica\n\
maize                                    | Zea mays\n\
(maize5)                                 | Zea mays\n\
\n"
  
#define HELP_USAGE "usage:\n\
augustus [parameters] --species=SPECIES queryfilename\n\
\n\
'queryfilename' is the filename (including relative path) to the file containing the query sequence(s)\n\
in fasta format.\n\
\n\
SPECIES is an identifier for the species. Use --species=help to see a list.\n\
\n\
parameters:\n\
--strand=both, --strand=forward or --strand=backward\n\
--genemodel=partial, --genemodel=intronless, --genemodel=complete, --genemodel=atleastone or --genemodel=exactlyone\n\
  partial      : allow prediction of incomplete genes at the sequence boundaries (default)\n\
  intronless   : only predict single-exon genes like in prokaryotes and some eukaryotes\n\
  complete     : only predict complete genes\n\
  atleastone   : predict at least one complete gene\n\
  exactlyone   : predict exactly one complete gene\n\
--singlestrand=true\n\
  predict genes independently on each strand, allow overlapping genes on opposite strands\n\
  This option is turned off by default.\n\
--hintsfile=hintsfilename\n\
  When this option is used the prediction considering hints (extrinsic information) is turned on.\n\
  hintsfilename contains the hints in gff format.\n\
--AUGUSTUS_CONFIG_PATH=path\n\
  path to config directory (if not specified as environment variable)\n\
--alternatives-from-evidence=true/false\n\
  report alternative transcripts when they are suggested by hints\n\
--alternatives-from-sampling=true/false\n\
  report alternative transcripts generated through probabilistic sampling\n\
--sample=n\n\
--minexonintronprob=p\n\
--minmeanexonintronprob=p\n\
--maxtracks=n\n\
  For a description of these parameters see section 2 of RUNNING-AUGUSTUS.md.\n\
--proteinprofile=filename\n\
  When this option is used the prediction will consider the protein profile provided as parameter.\n\
  The protein profile extension is described in section 5 of RUNNING-AUGUSTUS.md.\n\
--progress=true\n\
  show a progressmeter\n\
--gff3=on/off\n\
  output in gff3 format\n\
--predictionStart=A, --predictionEnd=B\n\
  A and B define the range of the sequence for which predictions should be found.\n\
--UTR=on/off\n\
  predict the untranslated regions in addition to the coding sequence. This currently works only for a subset of species.\n\
--noInFrameStop=true/false\n\
  Do not report transcripts with in-frame stop codons. Otherwise, intron-spanning stop codons could occur. Default: false\n\
--noprediction=true/false\n\
  If true and input is in genbank format, no prediction is made. Useful for getting the annotated protein sequences.\n\
--uniqueGeneId=true/false\n\
  If true, output gene identifyers like this: seqname.gN\n\
\n\
For a complete list of parameters, type \"augustus --paramlist\". A description of the important ones can be found in the file RUNNING-AUGUSTUS.md.\n"
 
#define HELP_USAGE_ETRAINING "usage:\n\
etraining --species=SPECIES trainfilename\n\
\n\
'trainfilename' is the filename (including relative path) to the file in genbank format containing the training sequences. These can be multi-gene sequences.\n\
SPECIES is a name for your species.\n\
\n\
parameters:\n\
--/genbank/verbosity=n\n\
  Choose one of 0,1,2 or 3. The larger the verbosity, the more (error) messages you get.\n"
 
 
#define POWER4TOTHE(x) ((int) (pow((double) 4,(double) (x)) + 0.1))
 
/*
 * bonus for (almost sure) hint with negative score
 * the user sets an (almost sure) anchor
 */
 
#define BONUSSUREHINT 1e100
 
typedef long int    Integer;
typedef int         Short;
typedef bool        Boolean;
typedef LLDouble Double;
typedef double      Float;
 
class Constant {
public:
    static void init();
    static string fullSpeciesPath() {return configPath + speciesDir;}
    static string modelPath() {return configPath + MODEL_SUBDIR;}
    static string extrinsicPath() {return configPath + EXTRINSIC_SUBDIR;}
 
    // size of splicesite patterns
    static int ass_size() { return ass_start + ass_end; }           // excluding "AG"
    static int ass_whole_size() { return ass_size() + ASS_MIDDLE; } // including "AG"
    static int dss_size() { return dss_start + dss_end; }           // excluding "GT"
    static int dss_whole_size() { return dss_size() + DSS_MIDDLE; } // including "GT"
    
    static string configPath;
    static string speciesDir;
    static int decomp_num_at;
    static int decomp_num_gc;
    static int decomp_num_steps;
    static int trans_init_window;
    static int tis_maxbinsize;
    static int ass_upwindow_size;
    static int init_coding_len;
    static int et_coding_len;
    static int ass_start;    //  size of ASS pattern before "AG" 
    static int ass_end;      //                      after  "AG"
    static int dss_start;    //  size of DSS pattern before "GT" 
    static int dss_end;      //                      after  "GT"
    static int dss_maxbinsize;
    static int ass_maxbinsize;
    static int tss_upwindow_size;   // win from start of model up to tss
    static int tss_start;
    static int min_coding_len;
    static int max_exon_len;
    static Integer d_polyasig_cleavage;
    static bool keep_viterbi;
    static double gc_range_min;
    static double gc_range_max;
    static double probNinCoding;
    static double opalprob;
    static double amberprob;
    static double ochreprob;
    static bool utr_option_on;
    static bool nc_option_on;
    static Integer augustus_verbosity;
    static bool alternatives_from_evidence;
    static double subopt_transcript_threshold;
    static Integer almost_identical_maxdiff;  // maximum allowed difference for the ends of the transcript
    static bool uniqueGeneId;
    static double max_contra_supp_ratio;
    static bool reportUtrOnlyGenes;
    static bool useCRFtraining;
    static bool CRFtrainTIS;
    static bool CRFtrainSS;
    static bool CRFtrainIntron;
    static bool CRFtrainIgenic;
    static bool CRFtrainCDS;
    static bool CRFtrainUTR;
    static bool dss_gc_allowed;
    static Boolean tieIgenicIntron; // whether to tie igenic model parameters to intron model parameters, i.e. use just one content model, that of the intron
    static Boolean exoncands;
    static Boolean proteinOutput;
    static Boolean codSeqOutput;
    static Boolean contentmodels; // whether to use content models, default: true
    static Integer min_intron_len; // minimal intron length (for hints)
    static bool MultSpeciesMode; // whether we do comparative gene prediction in multiple species
    static string  treefile; // file name in which a tree is specified in Newick format
    static string speciesfilenames; // file name to file which contains the names of species and the corresponding file names
    static string dbaccess; // comma separated string with database access (hostname, database name, table name, user, passwd
    static string alnfile; // name of file that contains MSA of genomes
    static string codonalnfile; // name of file that contains MSA of codon sequences
    static bool overlapmode; // whether overlapping exons are allowed in Viterbi algorithm
    static Boolean printEvidence; // "evidence for and against"
    static Boolean printExonCandsMSA; // print the OE nucleotide alignments
    static Boolean printOEs; // output ortho exons to file
    static Boolean printHints;
    static Boolean printMEA; // output .mea files (base genes) during CGP
    static Boolean printSampled; // output .sampled_GFs during CGP
    static Boolean printGeneRangesBED;
    static Boolean printGeneRangesGFF;
    static Integer maxOvlp; // parameters for overlapping coding regions in bacteria
    static vector<Double> head2tail_ovlp;
    static vector<Double> head2head_ovlp;
    static vector<Double> tail2tail_ovlp;
    static unsigned temperature; // heating the distribution for sampling, 0=cold, 7=hottest, take probs to the power of (8-temperature)/8
    static bool softmasking; // if true, lower-case character regions
                             // give rise to nonexonpart hints
    static bool softmasking_explicitly_requested; // i.e. it was not
                                                 // left at default value
    static bool dbhints;
    // scores from logistic regression  
    static bool logreg;
    static string referenceFile;
    static string refSpecies;
    static double GD_stepsize;
    static bool rLogReg;
    static double label_flip_prob;
    static double lambda;
    #ifdef COMPGENEPRED
    static unordered_map<string, pair<int, vector<double> > > logReg_feature;
    #endif
    static vector<double> ex_sc;
    static vector<double> in_sc;
    static int oeExtensionWidth;
    static vector<double> lg_es;
    static bool computeNumSubs;
    static bool useAArates;
    static bool useNonCodingModel;
    static bool rescaleBoni;
};
 
 
extern bool inCRFTraining;
 
extern const int power2id[31];
 
#define A_SET_FLAG(x) power2id[x]
 
struct Bitmask {
    Bitmask (int n=0) : value(n) {}
    Bitmask (const Bitmask& other) : value(other.value) {}
    bool operator[] (int n) const {
    return value & A_SET_FLAG(n);
    }
    void set(int n) {
    value |= A_SET_FLAG(n);
    }
    void unset(int n) {
    value &= ~A_SET_FLAG(n);
    }
    static Bitmask any() {
    return Bitmask(-1);
    }
    int value;
};
    
 
class ProjectError : public exception {
public:
    ProjectError( string msg ) throw () : exception(), message( msg ) {
    }
    ProjectError( ) throw () : exception() { }
 
    ~ProjectError() throw() {}
 
    string getMessage( ) { return message; }
private:
    string message;
};
 
class InvalidNucleotideError : public ProjectError {
public:
    InvalidNucleotideError( char t ) : ProjectError( "Invalid nucleotide '" + string(1,t) + "' encountered." ) { }
};
 
 
#define NUM_TYPES 86
 
enum StateType{TYPE_UNKNOWN = -1, igenic, 
           // forward strand
           singleG, initial0, initial1, initial2, internal0, internal1, internal2, terminal,
           lessD0, longdss0, equalD0, geometric0, longass0,  // The five intron states for frame 0
           lessD1, longdss1, equalD1, geometric1, longass1,  // The five intron states for frame 1
           lessD2, longdss2, equalD2, geometric2, longass2,  // The five intron states for frame 2
           utr5single, utr5init, utr5intron, utr5intronvar, utr5internal, utr5term, // 5'UTR states
           utr3single, utr3init, utr3intron, utr3intronvar, utr3internal, utr3term, // 3'UTR states
           // reverse strand
           rsingleG, rinitial, rinternal0, rinternal1, rinternal2, rterminal0, rterminal1, rterminal2,
           rlessD0, rlongdss0, requalD0, rgeometric0, rlongass0,  // The five intron states for frame 0
           rlessD1, rlongdss1, requalD1, rgeometric1, rlongass1,  // The five intron states for frame 1
           rlessD2, rlongdss2, requalD2, rgeometric2, rlongass2,  // The five intron states for frame 2
           rutr5single, rutr5init, rutr5intron, rutr5intronvar, rutr5internal, rutr5term, // 5'UTR states
           rutr3single, rutr3init, rutr3intron, rutr3intronvar, rutr3internal, rutr3term, // 3'UTR states
           
           intron_type, rintron_type, exon_type,
           // non-protein-coding states
           ncsingle, ncinit, ncintron, ncintronvar, ncinternal, ncterm, // forward strand
           rncsingle, rncinit, rncintron, rncintronvar, rncinternal, rncterm // reverse strand
};
 
extern const char* stateTypeNames[NUM_TYPES];
extern const char* stateTypeIdentifiers[NUM_TYPES];
extern const int stateReadingFrames[NUM_TYPES];
 
 
 
/*
 * mod3, returns the number in {0,1,2} being congruent to the argument modulo 3
 */
inline int mod3 (int k) {
    return (k>=0)? k%3 : (k%3+3)%3;
}
 
inline int modm(int k, int m) {
    return (k>=0)? k%m : (k%m+m)%m;
}
 
 
/* 
 * state type functions
 */
 
inline StateType toStateType( const char* str ){
    int i;
    for (i=0; i<NUM_TYPES; i++)
    if (strcmp(str, stateTypeIdentifiers[i]) == 0)
        return (StateType) i;
    return TYPE_UNKNOWN;
}
 
inline Boolean isInitialExon(StateType type){
    return (type==initial0 || type==initial1 || type==initial2);
}
 
inline Boolean isInternalExon(StateType type){
    return (type==internal0 || type==internal1 || type==internal2);
}
 
inline Boolean isRTerminalExon(StateType type){
    return (type==rterminal0 || type==rterminal1 || type==rterminal2);
}
 
inline Boolean isRInternalExon(StateType type){
    return (type==rinternal0 || type==rinternal1 || type==rinternal2);
}
 
inline Boolean isFirstExon(StateType type){
    return (isInitialExon(type) || isRTerminalExon(type) || type==singleG || type==rsingleG);
}
 
inline Boolean isLastExon(StateType type){
    return (type==terminal || type==rinitial || type==singleG || type == rsingleG);
}
 
inline Boolean is5UTR(StateType type){
    return (utr5single <= type && type <= utr5term)
    || (rutr5single <= type && type <= rutr5term);
}
 
inline Boolean is5UTRIntron(StateType type){
    return (type == utr5intron || type == utr5intronvar 
        || type == rutr5intron || type == rutr5intronvar);
}
 
inline Boolean is5UTRExon(StateType type){
    return is5UTR(type) && !is5UTRIntron(type);
}
 
inline Boolean is3UTR(StateType type){
    return  (utr3single <= type && type <= utr3term)
    || (rutr3single <= type && type <= rutr3term);
}
 
inline Boolean is3UTRIntron(StateType type){
    return (type == utr3intron || type == utr3intronvar 
        || type == rutr3intron || type == rutr3intronvar);
}
 
inline Boolean is3UTRExon(StateType type){
    return is3UTR(type) && !is3UTRIntron(type);
}
 
inline Boolean isFirstUTRExon(StateType type){
    return (type == utr5single || type == utr5init || type == rutr3single || type == rutr3term);
}
 
inline Boolean isLastUTRExon(StateType type){
    return (type == utr3single || type == utr3term || type == rutr5single || type == rutr5init);
}
 
inline Boolean isLongUTRIntron(StateType type){
    return (type == utr5intronvar || type == utr3intronvar || type == rutr5intronvar || type == rutr3intronvar);
}
 
inline Boolean isCodingExon(StateType type) {
    return  (singleG <= type && type <= terminal) 
    || (rsingleG <= type && type <= rterminal2);
}
 
inline Boolean isNcIntron(StateType type) {
    return  (type == ncintron || type == ncintronvar || type == rncintron || type == rncintronvar);
}
 
inline Boolean isNcExon(StateType type) {
    return  (ncsingle <= type && !isNcIntron(type));
}
 
inline Boolean isNc(StateType type) {
    return  (ncsingle <= type);
}
 
inline Boolean isExon(StateType type) {
    return (isCodingExon(type) || is5UTRExon(type) || is3UTRExon(type) || isNcExon(type));
}
 
inline Boolean isCodingIntron(StateType type) {
    return  (lessD0 <= type && type <= longass2)
    || (rlessD0 <= type && type <= rlongass2);
}
 
inline Boolean isIntron(StateType type) {
    return (isCodingIntron(type) || is5UTRIntron(type) || is3UTRIntron(type) 
        || type == intron_type || type == rintron_type || isNcIntron(type));
}
 
inline Boolean isGeometricIntron(StateType type){
    return (type==geometric0 || type==geometric1 || type==geometric2);
}
 
inline Boolean isRGeometricIntron(StateType type){
    return (type==rgeometric0 || type==rgeometric1 || type==rgeometric2);
}
 
inline Boolean isLongAssIntron(StateType type){
    return (type==longass0 || type== longass1 || type==longass2);
}
 
inline Boolean isRLongDssIntron(StateType type){
    return (type==rlongdss0 || type== rlongdss1 || type==rlongdss2);
}
 
inline Boolean isLongDssIntron(StateType type){
    return (type==longdss0 || type== longdss1 || type==longdss2);
}
 
inline Boolean isOnFStrand(StateType type){
    return ((type >= singleG && type < rsingleG) || (type >= ncsingle && type < rncsingle));
}
 
inline StateType initialExon(int frame){
    frame = mod3(frame);
    if (frame == 0)
    return initial0;
    if (frame == 1)
    return initial1;
    if (frame == 2)
    return initial2;
    return TYPE_UNKNOWN;
}
 
inline StateType internalExon(int frame){
    frame = mod3(frame);
    if (frame == 0)
    return internal0;
    if (frame == 1)
    return internal1;
    if (frame == 2)
    return internal2;
    return TYPE_UNKNOWN;
}
 
inline StateType lessDIntron(int frame){
    frame = mod3(frame);
   if (frame == 0)
    return lessD0;
    if (frame == 1)
    return lessD1;
    if (frame == 2)
    return lessD2;
    return TYPE_UNKNOWN;
}
 
 
inline StateType geometricIntron(int frame){
    frame = mod3(frame);
    if (frame == 0)
    return geometric0;
    if (frame == 1)
    return geometric1;
    if (frame == 2)
    return geometric2;
    return TYPE_UNKNOWN;
}
 
inline StateType longdssIntron(int frame){
    frame = mod3(frame);
    if (frame == 0)
    return longdss0;
    if (frame == 1)
    return longdss1;
    if (frame == 2)
    return longdss2;
    return TYPE_UNKNOWN;
}
 
inline StateType longassIntron(int frame){
    frame = mod3(frame);
    if (frame == 0)
    return longass0;
    if (frame == 1)
    return longass1;
    if (frame == 2)
    return longass2;
    return TYPE_UNKNOWN;
}
 
inline StateType equalDIntron(int frame){
    frame = mod3(frame);
    if (frame == 0)
    return equalD0;
    if (frame == 1)
    return equalD1;
    if (frame == 2)
    return equalD2;
    return TYPE_UNKNOWN;
}
 
inline StateType rterminalExon(int frame){
    frame = mod3(frame);
    if (frame == 0)
    return rterminal0;
    if (frame == 1)
    return rterminal1;
    if (frame == 2)
    return rterminal2;
    return TYPE_UNKNOWN;
}
 
inline StateType rinternalExon(int frame){
    frame = mod3(frame);
    if (frame == 0)
    return rinternal0;
    if (frame == 1)
    return rinternal1;
    if (frame == 2)
    return rinternal2;
    return TYPE_UNKNOWN;
}
inline StateType rlongdssIntron(int frame){
    frame = mod3(frame);
    if (frame == 0)
    return rlongdss0;
    if (frame == 1)
    return rlongdss1;
    if (frame == 2)
    return rlongdss2;
    return TYPE_UNKNOWN;
}
 
inline StateType rlongassIntron(int frame){
    frame = mod3(frame);
    if (frame == 0)
    return rlongass0;
    if (frame == 1)
    return rlongass1;
    if (frame == 2)
    return rlongass2;
    return TYPE_UNKNOWN;
}
 
inline StateType rlessDIntron(int frame){
    frame = mod3(frame);
   if (frame == 0)
    return rlessD0;
    if (frame == 1)
    return rlessD1;
    if (frame == 2)
    return rlessD2;
    return TYPE_UNKNOWN;
}
 
inline StateType rgeometricIntron(int frame){
    frame = mod3(frame);
    if (frame == 0)
    return rgeometric0;
    if (frame == 1)
    return rgeometric1;
    if (frame == 2)
    return rgeometric2;
    return TYPE_UNKNOWN;
}
 
inline StateType requalDIntron(int frame){
    frame = mod3(frame);
    if (frame == 0)
    return requalD0;
    if (frame == 1)
    return requalD1;
    if (frame == 2)
    return requalD2;
    return TYPE_UNKNOWN;
}
 
int howOftenOccursIt(const char* haystack, const char* needle, const char* endhaystack=NULL);
bool containsJustNonNucs(const char *dna, int dnalen);
bool isNuc(const char *dna);
 
/*
 * expand the ~ to the $Home directory
 */
string expandHome(string filename);
 
/*
 * convert int to string
 */
inline string itoa(int n) {
    ostringstream strm;
    strm << n;
    return strm.str();
}
 
/*
 * convert double to string
 */
inline string ftoa(double n) {
    ostringstream strm;
    strm << n;
    return strm.str();
}
 
/*
 * copy string to newly allocated character array
 */
inline char* newstrcpy(const char* s, int len) {
    char* result = new char[len + 1];
    result[len] = '\0';
    return strncpy(result, s, len);
}
inline char* newstrcpy(const char* s) {
    int len = strlen(s);
    char* result = new char[len + 1];
    return strcpy(result, s);
}
 
inline char* newstrcpy(const string& s, int len) {
    char* result = new char[len + 1];
    s.copy(result, len);
    result[len] = '\0';
    return result;
}
inline char* newstrcpy(const string& s) {
    return newstrcpy(s, s.length());
}
 
inline void strip_newlines(char* p) {
    char* q = p;
    while (p != 0 && *p != '\0') {
        if (*p == '\n') {
            p++;
            *q = *p;
        } 
        else {
            *q++ = *p++;
        }
    }
    *q = '\0';
}
 
char *getRandomDNA(int len);
 
/* for 0 <= q <=1 get the q-th quantile of the values store in
 * vector v. 
 */
Double quantile(const vector<Double> &v, float q);
int quantile(const vector<int> &v, float q);
 
/* obtain a "hashtable" from a list of unique names, e.g.
 * [human, mouse, dog] becomes
 * human => 0, mouse => 1, dog => 2
 * Mario: This could become an unordered_map (expected constant time) once we commit to C++11,
 * rather than a logarithmic tree data structure
 */
map<string, size_t> *getMap (vector<string> names);
 
/* Structure to store and accumulate a vector of doubles.
 * Used for storing cumulative sums of log likelihoods for each rate matrix (e.g. for each omega)
 * (optional: cumulative sum of number of substitutions) 
 * we create an object of cumValues for each bit_vector and reading frame combination
 */
struct cumValues{
    vector<double> logliks;
    int numSubs;
    int id;
    int numSites; // number of codon alignment columns
    vector<string> rows; // for debugging, temporary
    
    cumValues(int i=0, int s=-1):
        numSubs(s), id(i), numSites(0) {};
    
    void addLogliks(vector<double>* ll){
        if (logliks.size() == 0){
            logliks.resize(ll->size(), 0.0);
        }
        for (unsigned u = 0; u < ll->size(); u++){
            logliks[u] += (*ll)[u];
        }
        numSites++;
    }
    void addNumSubs(int subs){
        numSubs += subs;
    }
    void addRows(vector<string> &cs) {
        if (rows.size() == 0)
            rows.assign(cs.size(), "");
        for (unsigned i=0; i < cs.size(); i++){
            rows[i].append(cs[i]);
        }
    }
};
 
/*
 * functions used in earlier versions
 
 *
 * sort an array increasingly p from index a to index b (included).
 */
// void QuickSort(Double *p, int a, int b);
 
/*
 * compute the unbiased empirical standard deviation of a distribution given by an array of doubles of size n
 * a distribution on the numbers {0,1, ..., n-1}
 */
// double variance(double *p, int n);
 
/*
 * compute the chi-square statistics for 
 * a multinomial distribution on the numbers {0,1, ..., r-1} given by p
 * reference: uniform distribution
 */
// double chiSquareUniform(int *p, int n);
 
// double chisquare(int a[4][4]);
// int branchPointPosition(char *seq, int asspos);
 
#endif   //  _TYPES_HH