// This file is part of BULL, a program for phylogenetic simulations // most of the code was written by Mark T. Holder. // This program is for internal use by the lab of Dr. Tandy Warnow only. // Do not redistribute the code. It is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. // // Some of the code is from publically available source by Paul Lewis, Ziheng Yang, // John Huelsenbeck, David Swofford , and others (as noted in the code). // In fact the main structure of the program was created by modifying Paul Lewis' // basiccmdline.cpp from his NCL // // This code was used in Mark's dissertation, some changes were made in order to // get it to compile on gcc. It is possible that this porting introduced bugs (very // little debugging has been done on UNIX platforms). I would suggest checking // the simulator by generating data on trees with short branches, etc. // This file is part of BULL, a program for phylogenetic simulations // most of the code was written by Mark T. Holder. // This program is for internal use by the lab of Dr. Tandy Warnow only. // Do not redistribute the code. It is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. // // Some of the code is from publically available source by Paul Lewis, Ziheng Yang, // John Huelsenbeck, David Swofford , and others (as noted in the code). // In fact the main structure of the program was created by modifying Paul Lewis' // basiccmdline.cpp from his NCL // // This code was used in Mark's dissertation, some changes were made in order to // get it to compile on gcc. It is possible that this porting introduced bugs (very // little debugging has been done on UNIX platforms). I would suggest checking // the simulator by generating data on trees with short branches, etc. #ifndef SSRFSUBMODCODON #define SSRFSUBMODCODON #include #include #if defined HAVE_CONFIG_H && HAVE_CONFIG_H # include #endif #define MINIMUMAAFREQABOVEZERO 1.0e-12 #include "model.hpp" #include "ssrf_codon_sub_mod.hpp" namespace bull { class SSRFCodonSubMod : public ModelWEig { double CalculateBeta(void); FreqParamGroup *baseFreqs,*aaFreqs; double **codFreqs,SSRkContrib,overflowMultiplier; Parameter *sharedBrLenInterpreter; void CalculateCodonFreqs(); //global codon code: AAA=0 ... TTT=63 global amino acid code A=0 ... Y=19 *=20 int *codonLocToGlob;//codonTrans[i] is the global codon code where i is the local codon code int *aminoAcidLocToGlob;//aminoAcidTLocToGlob[i] is the global amino acid code where i is the local amino acid code int codonGlobToLoc[64];//thisCodeByStd[i] is the local codon code where i is the global code int *codLocToAminoAcidLoc;//translationToLocal[i] is the local amino acid code where i is the local codon code // in the cpp file there is a global int MitotcodNum[i] that is the global amino acid code when i is the global codon code int nPossAAs; double **mutMat; public: static bool mutParamsDirty; static int geneticCode; static double defFreqA,defFreqC,defFreqG,defrAC,defrAG,defrAT,defrCG,defrCT,defrGT, defA,defC,defD,defE,defF,defG,defH,defI,defK,defL,defM,defN,defP,defQ,defR,defS,defT, defV,defW,defY,defStop,defRate; int nCodonsInProtein;//total # of codons, in the models that share this mutational set of parameters #ifdef ALLOWMULTIHITS enum params {blenMult=0 ,pMultHit,freqA , freqC , freqG , freqT , rAC , rAG , rAT , rCG , rCT , rGT ,fAAs}; #else enum params {blenMult=0 ,freqA , freqC , freqG , freqT , rAC , rAG , rAT , rCG , rCT , rGT ,fAAs}; #endif #ifdef ELIMINATE_ALL_ZEROS int *maxcodonLocToGlob; int *maxcodLocToAminoAcidLoc; int *maxaminoAcidLocToGlob; int origCodonGlobToLoc[64]; int maxNStates,maxNPossAAs; Parameter **OrigParamArray; #endif SSRFCodonSubMod(int ncods,int naas,bool *codP,bool *aaP,Parameter **gtrParams,FreqParamGroup *nucleoFreqs,double *ssrfParams,double *vecPreAlloc,int codnu,double **mm); virtual ~SSRFCodonSubMod(); void CalculateQ(); void DeleteSharedMemory(); //written to allow odd coding of SSRFCodonSubModel which overrides void EncodeACharacter(short *dest,short *inp,int datatype,bool keepGap) ; void FreqParamChangesShouldSumToOne(FreqParamGroup *p); int GetEncodingType(); FreqParamGroup *GetFreqParamGroup(int i); int GetNStates(); double **GetStateFreqs(); void InitializeParameters(); int NumShortsPerCharacter(); int NumStatesInLastShort(); void ParameterHasChanged(FreqParamGroup *p) ; void ParameterHasChanged(Parameter *p); void RecalculateMutMatrix(); static void set_default(double *allParams); void SharedConstruction(bool *codP,bool *aaP, FreqParamGroup *nucleoFreqs,double *ssrfParams,int codnu); inline void UpdatePmat(double b); inline void UpdatePmat(double b,int onlycol); inline void UpdatePMatrix(double **pm,double b);//for externally supplied Pmatrix inline void UpdatePMatrix(double **pm,double b,int onlycol); inline void UpdatePRow(double **pm,double b,int onlycol); inline void UpdatePmatWithOutSharedMatrix(double **pmats,double b); inline void UpdatePmatWithOutSharedMatrix(double **pmats,double b,int onlycol); //AlertSharedMemory() is called because the EigInvEigMult matrix is shared, if the eigen vectors don't need to be recalculated //(the mod parameter haven't changed) this will need to be called whenever the model being scored is changed called by Tree::PrepareLikeAttrAndAllNodes(int partnum,int modn) void AlertSharedMemory(); double GetMultiplier() const; void SetMultiplier(double x); void AddAminoAcidFreqs(std::ofstream &dest); #ifdef ELIMINATE_ALL_ZEROS void ResizeModel(); bool NeedToExpandPossibleAA(); int GetMaxNStates(); int *GetOrigLocToGlob(); int *GetCurrentGlobToLoc(); int *GetCurrentLocToGlob(); #endif void SetBlenMultBasedOnSBLI(void); void Summarize(std::ofstream &outpf); void Fill4TaxonSpectrum(double *spec,double bOne,double bTwo,double bInt,double bThr,double bFou); void Fill4TaxonSpectrum(double *specOne,double *specTwo, double *specThr,double bOne,double bTwo,double bInt,double bThr,double bFou);//for testing site specific rates double GetDeviationFromEquilibriumFreq(double branchLen); double GetProportionOnSubOptimalHill(); int GetIndexOfMutationalNeighbor(int inIndex,int numOfMutNeighbor); }; void CreateSetOfSSRFCodonSubMods(int nm,double *gtrp/* at least 8 doubles long */,double **setofSSRFparams/* full size=nm*21 */,SSRFCodonSubMod **arrayOfMods,double *multipliers,double treesc,int gcode); void ModifyBlenMultSoBranchesAreScaled(SSRFCodonSubMod **modArr,int ns); int GetNumberOfPossibleAminoAcids(bool *codOfObsAA,bool *possAAs); int GetNumberOfPossibleCodons(bool *possAAs); void SetCodeRelatedGlobals(int gcode); void GetNeutralAminoAcidFreq(double *aafreq,double *baseFreqParams); //bool CheckNotZero(double x); int GetGlobalIndexOfMutationalNeighbor(int inIndex,int numOfMutNeighbor); } //namespace bull { #endif