src/utils.h

   1 #ifndef __UTILS_H__
   2 #define __UTILS_H__
   3
   4 #undef complex
   5 struct complex
   6 {
   7     MrBFlt re;
   8     MrBFlt im;
   9 };
  10 typedef struct complex complex;
  11
  12 typedef struct
  13     {
  14     MrBFlt     mean;
  15     MrBFlt     median;
  16     MrBFlt     lower;
  17     MrBFlt     upper;
  18     MrBFlt     var;
  19     MrBFlt     PSRF;
  20     MrBFlt     avrESS;
  21     MrBFlt     minESS;
  22     }
  23     Stat;
  24
  25 int      AddBitfield (BitsLong ***list, int listLen, int *set, int setLen);
  26 #if defined (SSE_ENABLED)
  27 void    *AlignedMalloc (size_t size, size_t alignment);
  28 void     AlignedSafeFree (void **ptr);
  29 #endif
  30 int      AreBitfieldsEqual (BitsLong *p, BitsLong *q, int length);
  31 int      Bit (int n, BitsLong *p);
  32 void     ClearBit (int i, BitsLong *bits);
  33 void     ClearBits (BitsLong *bits, int nLongs);
  34 void     CopyBits (BitsLong *dest, BitsLong *source, int nLongs);
  35 int      CopyResults (FILE *toFile, char *fromFileName, int lastGen);
  36 int      CopyProcessSsFile (FILE *toFile, char *fromFileName, int lastStep, MrBFlt *marginalLnLSS, MrBFlt *splitfreqSS);
  37 int      CopyTreeResults (FILE *toFile, char *fromFileName, int lastGen, int *treeNum);
  38 int      FirstTaxonInPartition (BitsLong *partition, int length);
  39 long     FirstTree (FILE *fp, char *lineBuf, int longestLine);
  40 int      Flip01 (int x);
  41 void     FlipBits (BitsLong *partition, int length, BitsLong *mask);
  42 void     FlipOneBit (int n, BitsLong *p);
  43 int      FromGrowthFxnToIndex (int *growthFxn);
  44 void     FromIndexToGrowthFxn (int index, int *growthFxn);
  45 void     GetIntSummary (int **vals, int nRows, int *rowCount, Stat *theStats, int HPD);
  46 int      GetKFromGrowthFxn (int *growthFxn);
  47 void     GetSummary (MrBFlt **vals, int nRows, int *rowCount, Stat *theStats, int HPD);
  48 int      HarmonicArithmeticMeanOnLogs (MrBFlt *vals, int nVals, MrBFlt *mean, MrBFlt *harm_mean);
  49 int      IsBitSet (int i, BitsLong *bits);
  50 int      IsConsistentWith (const char *token, const char *expected);
  51 int      IsPartNested (BitsLong *smaller, BitsLong *larger, int length);
  52 int      IsPartCompatible (BitsLong *smaller, BitsLong *larger, int length);
  53 int      IsSectionEmpty (BitsLong *bitField1, BitsLong *bitField2, int length);
  54 int      IsUnionEqThird (BitsLong *bitField1, BitsLong *bitField2, BitsLong *bitField3, int length);
  55 long     LastBlock (FILE *fp, char *lineBuf, int longestLine);
  56 int      LineTermType (FILE *fp);
  57 int      LongestLine (FILE *fp);
  58 void     LowerUpperMedian (MrBFlt *vals, int nVals, MrBFlt *lower, MrBFlt *upper, MrBFlt *median);
  59 void     LowerUpperMedianHPD (MrBFlt *vals, int nVals, MrBFlt *lower, MrBFlt *upper, MrBFlt *median);
  60 void     MeanVariance (MrBFlt *vals, int nVals, MrBFlt *mean, MrBFlt *var);
  61 void     MeanVarianceLog (MrBFlt *vals, int nVals, MrBFlt *mean, MrBFlt *var, MrBFlt *varEst);
  62 int      NextTaxonInPartition (int currentTaxon, BitsLong *partition, int length);
  63 int      NBits (int x);
  64 int      NumBits (BitsLong *x, int len);
  65 char    *MbPrintNum (MrBFlt num);
  66 void     MrBayesPrint (char *format, ...);
  67 void     MrBayesPrintf (FILE *f, char *format, ...);
  68 FILE    *OpenBinaryFileR (char *name);
  69 FILE    *OpenTextFileA (char *name);
  70 FILE    *OpenTextFileR (char *name);
  71 FILE    *OpenTextFileRQuait (char *name);
  72 FILE    *OpenTextFileW (char *name);
  73 MrBFlt   PotentialScaleReduction (MrBFlt **vals, int nRows, int *count);
  74 void     EstimatedSampleSize (MrBFlt **vals, int nRuns, int *count, MrBFlt *returnESS);
  75 void    *SafeCalloc (size_t n, size_t s);
  76 int      SafeFclose (FILE **fp);
  77 void     SafeFree (void **ptr);
  78 void    *SafeMalloc (size_t s);
  79 void    *SafeRealloc (void *ptr, size_t s);
  80 char    *SafeStrcat (char **target, const char *source);
  81 char    *SafeStrcpy (char **target, const char *source);
  82 void     SetBit (int i, BitsLong *bits);
  83 void     SortInts (int *item, int *assoc, int count, int descendingOrder);
  84 void     SortInts2 (int *item, int *assoc, int left, int right, int descendingOrder);
  85 void     SortMrBFlt (MrBFlt *item, int left, int right);
  86 int      StrCmpCaseInsensitive (char *s, char *t);
  87 void     StripComments (char *s);
  88 FILE    *TestOpenTextFileR (char *name);
  89 void     UpdateGrowthFxn (int *growthFxn);
  90 int      UpperTriangIndex (int i, int j, int size);
  91 int      WantTo (const char *msg);
  92
  93 /* tree utility functions */
  94 int       AddToTreeList (TreeList *treeList, Tree *tree);
  95 Tree     *AllocateTree (int numTaxa);
  96 Tree     *AllocateFixedTree (int numTaxa, int isRooted);
  97 int       AllocateTreePartitions (Tree *t);
  98 PolyTree *AllocatePolyTree (int numTaxa);
  99 int       AllocatePolyTreePartitions (PolyTree *pt);
 100 int       AllocatePolyTreeRelClockParams (PolyTree *pt, int nBSets, int nESets);
 101 int       AreTopologiesSame (Tree *t1, Tree *t2);
 102 int       AreTreesSame (Tree *t1, Tree *t2);
 103 int       BuildConstraintTree (Tree *t, PolyTree *pt, char **localTaxonNames);
 104 int       BuildRandomRTopology (Tree *t, RandLong *seed);
 105 int       BuildRandomUTopology (Tree *t, RandLong *seed);
 106 int       CheckConstraints (Tree *t);
 107 int       CheckSetConstraints (Tree *t);
 108 void      ColorClusters (TreeNode *p, int *index);
 109 void      CopySubtreeToTree (Tree *subtree, Tree *t);
 110 int       CopyToPolyTreeFromPolyTree (PolyTree *to, PolyTree *from);
 111 int       CopyToSpeciesTreeFromPolyTree (Tree *to, PolyTree *from);
 112 int       CopyToTreeFromPolyTree (Tree *to, PolyTree *from);
 113 void      CopyPolyNodes (PolyNode *p, PolyNode *q, int nLongsNeeded);
 114 int       CopyToTreeFromTree (Tree *to, Tree *from);
 115 void      CopyTreeNodes (TreeNode *p, TreeNode *q, int nLongsNeeded);
 116 void      CopyTreeToSubtree (Tree *t, Tree *subtree);
 117 int       Deroot (PolyTree *pt);
 118 void      EraseTreeList (TreeList *treeList);
 119 void      findAllowedClockrate (Tree *t, MrBFlt *minClockRate, MrBFlt *maxClockRate);
 120 void      FreePolyTree (PolyTree *pt);
 121 void      FreePolyTreePartitions (PolyTree *pt);
 122 void      FreePolyTreePopSizeParams (PolyTree *pt);
 123 void      FreePolyTreeRelClockParams (PolyTree *pt);
 124 void      FreeTree (Tree *t);
 125 void      FreeTreePartitions (Tree *pt);
 126 void      GetDatedNodeDepths (TreeNode *p, MrBFlt *nodeDepths);
 127 void      GetDatedNodes (TreeNode *p, TreeNode **datedNodes);
 128 void      GetDownPass (Tree *t);
 129 void      GetNodeDownPass (Tree *t, TreeNode *p, int *i, int *j);
 130 void      GetPolyAges (PolyTree *t);
 131 void      GetPolyDepths (PolyTree *t);
 132 void      GetPolyDownPass (PolyTree *t);
 133 void      GetPolyNodeDownPass (PolyTree *t, PolyNode *p, int *i, int *j);
 134 int       GetRandomEmbeddedSubtree (Tree *t, int nTerminals, RandLong *seed, int *nEmbeddedTrees);
 135 int       GetFromTreeList (TreeList *treeList, Tree *tree);
 136 int       InitBrlens (Tree *t, MrBFlt v);
 137 int       InitCalibratedBrlens (Tree *t, MrBFlt minLength, RandLong *seed);
 138 int       InitClockBrlens (Tree *t);
 139 int       IsCalibratedClockSatisfied (Tree *t,MrBFlt *minClockRate,MrBFlt *maxClockRate , MrBFlt tol);
 140 int       IsClockSatisfied (Tree *t, MrBFlt tol);
 141 int       IsTreeConsistent (Param *param, int chain, int state);
 142 int       LabelTree (Tree *t, char **taxonNames);
 143 void      Mark (TreeNode *p);
 144 void      MarkDistance (TreeNode *p, int YESorNO, int dist, int *n);
 145 void      MarkUnconstrained (TreeNode *p);
 146 int       MoveCalculationRoot (Tree *t, int outgroup);
 147 int       MovePolyCalculationRoot (PolyTree *t, int outgroup);
 148 int       NumConstrainedTips (TreeNode *p);
 149 int       NumDatedTips (TreeNode *p);
 150 void      OrderTips (PolyTree *t);
 151 void      PrintNewick (char **s, int *len, Tree *t);
 152 void      PrintNodes (Tree *t);
 153 void      PrintPolyNodes (PolyTree *pt);
 154 void      PrintTranslateBlock (FILE *fp, Tree *t);
 155 int       PrunePolyTree (PolyTree *pt);
 156 int       RandPerturb (Tree *t, int nPert, RandLong *seed);
 157 int       RandResolve (Tree *tt, PolyTree *t, RandLong *seed, int destinationIsRooted);
 158 int       ResetBrlensFromTree (Tree *tree, Tree *vTree);
 159 void      ResetIntNodeIndices (PolyTree *t);
 160 void      ResetPolyTree (PolyTree *t);
 161 void      ResetPolyTreePartitions (PolyTree *pt);
 162 void      ResetPolyTreeRelClockParams (PolyTree *pt);
 163 int       ResetRootHeight (Tree *t, MrBFlt rootHeight);
 164 void      ResetTipIndices (PolyTree *pt);
 165 int       ResetTopology (Tree *t, char *s);
 166 int       ResetTopologyFromTree (Tree *tree, Tree *top);
 167 int       ResetTopologyFromPolyTree (Tree *tree, PolyTree *top);
 168 void      ResetTreePartitions (Tree *t);
 169 int       RetrieveRTopology (Tree *t, int *order);
 170 int       RetrieveRTree (Tree *t, int *order, MrBFlt *brlens);
 171 int       RetrieveRTreeWithIndices (Tree *t, int *order, MrBFlt *brlens);
 172 int       RetrieveUTopology (Tree *t, int *order);
 173 int       RetrieveUTree (Tree *t, int *order, MrBFlt *brlens);
 174 void      SetDatedNodeAges (Param* param, int chain, int state);
 175 void      SetNodeDepths (Tree *t);
 176 int       SetTreeNodeAges (Param *param, int chain, int state);
 177 int       ShowPolyNodes (PolyTree *pt);
 178 int       ShowTree (Tree *t);
 179 int       StoreRPolyTopology (PolyTree *t, int *order);
 180 int       StoreRPolyTree (PolyTree *t, int *order, MrBFlt *brlens);
 181 int       StoreRTopology (Tree *t, int *order);
 182 int       StoreRTree (Tree *t, int *order, MrBFlt *brlens);
 183 int       StoreRTreeWithIndices (Tree *t, int *order, MrBFlt *brlens);
 184 int       StoreUPolyTopology (PolyTree *t, int *order);
 185 int       StoreUPolyTree (PolyTree *t, int *order, MrBFlt *brlens);
 186 int       StoreUTopology (Tree *t, int *order);
 187 int       StoreUTree (Tree *t, int *order, MrBFlt *brlens);
 188 MrBFlt    TreeLen (Tree *t);
 189 void      Unmark (TreeNode *p);
 190 void      UpdateTreeWithClockrate (Tree *t, MrBFlt clockRate);
 191 void      WriteEventTree (TreeNode *p, int chain, Param *param);
 192 void      WriteEventTreeToPrintString (TreeNode *p, int chain, Param *param, int printAll);
 193 void      WriteNoEvtTreeToPrintString (TreeNode *p, int chain, Param *param, int showBrlens, int isRooted);
 194 void      WriteEvolTree (TreeNode *p, int chain, Param *param);
 195 void      WriteTopologyToFile (FILE *fp, TreeNode *p, int isRooted);
 196
 197 /* math utility functions */
 198 complex **AllocateSquareComplexMatrix (int dim);
 199 MrBFlt  **AllocateSquareDoubleMatrix (int dim);
 200 int     **AllocateSquareIntegerMatrix (int dim);
 201 int       AutodGamma (MrBFlt *M, MrBFlt rho, int K);
 202 void      BetaBreaks (MrBFlt alpha, MrBFlt beta, MrBFlt *values, int K);
 203 MrBFlt    BetaQuantile (MrBFlt alpha, MrBFlt beta, MrBFlt x);
 204 void      CalcCijk (int dim, MrBFlt *c_ijk, MrBFlt **u, MrBFlt **v);
 205 void      CopyComplexMatrices (int dim, complex **from, complex **to);
 206 void      CopyDoubleMatrices (int dim, MrBFlt **from, MrBFlt **to);
 207 void      DirichletRandomVariable (MrBFlt *alp, MrBFlt *z, int n, RandLong *seed);
 208 int       DiscreteGamma (MrBFlt *rK, MrBFlt alfa, MrBFlt beta, int K, int median);
 209 void      FreeSquareComplexMatrix (complex **m);
 210 void      FreeSquareDoubleMatrix (MrBFlt **m);
 211 void      FreeSquareIntegerMatrix (int **m);
 212 int       GetEigens (int dim, MrBFlt **q, MrBFlt *eigenValues, MrBFlt *eigvalsImag, MrBFlt **eigvecs, MrBFlt **inverseEigvecs, complex **Ceigvecs, complex **CinverseEigvecs);
 213 MrBFlt    LnFactorial (int value);
 214 MrBFlt    LnGamma (MrBFlt alp);
 215 MrBFlt    LnPriorProbExponential (MrBFlt val, MrBFlt *params);
 216 MrBFlt    LnPriorProbExponential_Param_Mean (MrBFlt val, MrBFlt *params);
 217 MrBFlt    LnPriorProbFix (MrBFlt val, MrBFlt *params);
 218 MrBFlt    LnPriorProbGamma (MrBFlt val, MrBFlt *params);
 219 MrBFlt    LnPriorProbGamma_Param_Mean_Sd (MrBFlt val, MrBFlt *params);
 220 MrBFlt    LnPriorProbLognormal (MrBFlt val, MrBFlt *params);
 221 MrBFlt    LnPriorProbLognormal_Param_Mean_Sd (MrBFlt val, MrBFlt *params);
 222 MrBFlt    LnPriorProbNormal (MrBFlt val, MrBFlt *params);
 223 MrBFlt    LnPriorProbOffsetExponential (MrBFlt val, MrBFlt *params);
 224 MrBFlt    LnPriorProbOffsetExponential_Param_Offset_Mean (MrBFlt val, MrBFlt *params);
 225 MrBFlt    LnPriorProbOffsetGamma (MrBFlt val, MrBFlt *params);
 226 MrBFlt    LnPriorProbOffsetGamma_Param_Offset_Mean_Sd (MrBFlt val, MrBFlt *params);
 227 MrBFlt    LnPriorProbOffsetLognormal (MrBFlt val, MrBFlt *params);
 228 MrBFlt    LnPriorProbOffsetLognormal_Param_Offset_Mean_Sd (MrBFlt val, MrBFlt *params);
 229 MrBFlt    LnPriorProbTruncatedNormal (MrBFlt val, MrBFlt *params);
 230 MrBFlt    LnPriorProbTruncatedNormal_Param_Trunc_Mean_Sd (MrBFlt val, MrBFlt *params);
 231 MrBFlt    LnPriorProbUniform (MrBFlt val, MrBFlt *params);
 232 MrBFlt    LnProbRatioExponential (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 233 MrBFlt    LnProbRatioExponential_Param_Mean (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 234 MrBFlt    LnProbRatioFix (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 235 MrBFlt    LnProbRatioGamma (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 236 MrBFlt    LnProbRatioGamma_Param_Mean_Sd (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 237 MrBFlt    LnProbRatioLognormal (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 238 MrBFlt    LnProbRatioLognormal_Param_Mean_Sd (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 239 MrBFlt    LnProbRatioNormal (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 240 MrBFlt    LnProbRatioOffsetExponential (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 241 MrBFlt    LnProbRatioOffsetExponential_Param_Offset_Mean (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 242 MrBFlt    LnProbRatioOffsetGamma (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 243 MrBFlt    LnProbRatioOffsetGamma_Param_Offset_Mean_Sd (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 244 MrBFlt    LnProbRatioOffsetLognormal (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 245 MrBFlt    LnProbRatioOffsetLognormal_Param_Offset_Mean_Sd (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 246 MrBFlt    LnProbRatioTruncatedNormal (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 247 MrBFlt    LnProbRatioTruncatedNormal_Param_Trunc_Mean_Sd (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 248 MrBFlt    LnProbRatioUniform (MrBFlt newX, MrBFlt oldX, MrBFlt *params);
 249 MrBFlt    LnProbGamma (MrBFlt alpha, MrBFlt beta, MrBFlt x);
 250 MrBFlt    LnProbTruncGamma (MrBFlt alpha, MrBFlt beta, MrBFlt x, MrBFlt min, MrBFlt max);
 251 MrBFlt    LnProbLogNormal (MrBFlt exp, MrBFlt sd, MrBFlt x);
 252 MrBFlt    LnRatioLogNormal (MrBFlt exp, MrBFlt sd, MrBFlt xNew, MrBFlt xOld);
 253 MrBFlt    LnProbTK02LogNormal (MrBFlt mean, MrBFlt var, MrBFlt x);
 254 MrBFlt    LnRatioTK02LogNormal (MrBFlt exp, MrBFlt sd, MrBFlt xNew, MrBFlt xOld);
 255 MrBFlt    LogNormalRandomVariable (MrBFlt mean, MrBFlt var, RandLong *seed);
 256 MrBFlt    MaximumValue (MrBFlt x, MrBFlt y);
 257 MrBFlt    MinimumValue (MrBFlt x, MrBFlt y);
 258 void      MultiplyMatrices (int dim, MrBFlt **a, MrBFlt **b, MrBFlt **result);
 259 int       MultiplyMatrixNTimes (int dim, MrBFlt **Mat, int power, MrBFlt **Result);
 260 MrBFlt    PointNormal (MrBFlt prob);
 261 MrBFlt    PsiGammaLnProb (MrBFlt alpha, MrBFlt value);
 262 MrBFlt    PsiGammaLnRatio (MrBFlt alpha, MrBFlt numerator, MrBFlt denominator);
 263 MrBFlt    PsiGammaRandomVariable (MrBFlt alpha, RandLong *seed);
 264 MrBFlt    QuantileGamma (MrBFlt x, MrBFlt alfa, MrBFlt beta);
 265 MrBFlt    RandomNumber (RandLong *seed);
 266 MrBFlt    QuantileLogNormal (MrBFlt prob, MrBFlt mu, MrBFlt sigma);
 267 int       DiscreteLogNormal (MrBFlt *rK, MrBFlt sigma, int K, int median);
 268 MrBFlt    LogNormalPoint (MrBFlt x, MrBFlt mu, MrBFlt sigma);
 269
 270 #endif  /* __UTILS_H__ */