Merge pull request #294 from ronquist/mpi-checkpoint-fix

Better checking of clock tree consistency
NBISweden · Feb 2, 2024 · e7b8431 · e7b8431
2 parents 3903b3b + 9c3f865
commit e7b8431
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 30 deletions.
diff --git a/src/mcmc.c b/src/mcmc.c
@@ -2584,6 +2584,10 @@ int DoMcmc (void)
             goto errorExit;
         }
 
+    /* Check again that the model is consistent after potentially resetting starting values (co-dependencies etc) */
+    if (CheckModel() == ERROR)
+        goto errorExit;
+
 #ifndef NDEBUG
     /* checking tree consistency for debug purposes */
     for (i=0; i<numParams; i++)

diff --git a/src/model.c b/src/model.c
@@ -1906,7 +1906,7 @@ void CheckCharCodingType (Matrix *m, CharInfo *ci)
 -------------------------------------------------------------*/
 int CheckModel (void)
 {
-    int         i, j, k, answer;
+    int         ch, i, j, k, answer;
     Tree        *t = NULL;
     TreeNode    *p;
     MrBFlt      treeAge, clockRate;
@@ -1948,44 +1948,48 @@ int CheckModel (void)
     /* 
      * Check that the clock rate is consistent with the tree age prior. We cannot check this earlier
      * because the clock rate and tree age are set separately, and we do not know in which order they are set.
+     * We need to check all chains because start values are set separately for each chain.
      */
     for (i=0; i<numTrees; i++)
         {
-        t = GetTreeFromIndex(i,0,0);
-        if (t->isClock == YES && t->isCalibrated == YES)
+        for (ch=0; ch<numGlobalChains; ch++)
             {
-            clockRate = *GetParamVals(modelSettings[t->relParts[0]].clockRate, 0, 0);
-            treeAge = t->root->left->nodeDepth / clockRate;
-            if (!AreDoublesEqual(treeAge, t->root->left->age, 0.000001))
-                {
-                MrBayesPrint("%s   ERROR: The tree age setting is inconsistent with the specified tree age prior.\n", spacer);
-                return (ERROR);
-                }
-            if (modelParams[t->relParts[0]].treeAgePr.prior == fixed)
+            t = GetTreeFromIndex(i,ch,0);
+            if (t->isClock == YES && t->isCalibrated == YES)
                 {
-                if (!AreDoublesEqual(treeAge, modelParams[t->relParts[0]].treeAgePr.priorParams[0], 0.000001))
+                clockRate = *GetParamVals(modelSettings[t->relParts[0]].clockRate, ch, 0);
+                treeAge = t->root->left->nodeDepth / clockRate;
+                if (!AreDoublesEqual(treeAge, t->root->left->age, 0.000001))
                     {
-                    MrBayesPrint("%s   ERROR: The clock rate is inconsistent with the specified tree age prior.\n", spacer);
+                    MrBayesPrint("%s   ERROR: The tree age setting is inconsistent with the specified tree age prior.\n", spacer);
                     return (ERROR);
                     }
-                }
-            else if (modelParams[t->relParts[0]].treeAgePr.prior == uniform)
-                {
-                if (treeAge < modelParams[t->relParts[0]].treeAgePr.priorParams[0] || treeAge > modelParams[t->relParts[0]].treeAgePr.priorParams[1])
-                    {    
-                    MrBayesPrint("%s   ERROR: The clock rate is inconsistent with the specified tree age prior.\n", spacer);
-                    return (ERROR);
+                if (modelParams[t->relParts[0]].treeAgePr.prior == fixed)
+                    {
+                    if (!AreDoublesEqual(treeAge, modelParams[t->relParts[0]].treeAgePr.priorParams[0], 0.000001))
+                        {
+                        MrBayesPrint("%s   ERROR: The clock rate is inconsistent with the specified tree age prior.\n", spacer);
+                        return (ERROR);
+                        }
                     }
-                }
-            else if (modelParams[t->relParts[0]].treeAgePr.prior == offsetExponential ||
-                     modelParams[t->relParts[0]].treeAgePr.prior == offsetGamma ||
-                     modelParams[t->relParts[0]].treeAgePr.prior == truncatedNormal ||
-                     modelParams[t->relParts[0]].treeAgePr.prior == offsetLogNormal)
-                {
-                if (treeAge < modelParams[t->relParts[0]].treeAgePr.priorParams[0])
-                    {    
-                    MrBayesPrint("%s   ERROR: The clock rate is inconsistent with the specified tree age prior.\n", spacer);
-                    return (ERROR);
+                else if (modelParams[t->relParts[0]].treeAgePr.prior == uniform)
+                    {
+                    if (treeAge < modelParams[t->relParts[0]].treeAgePr.priorParams[0] || treeAge > modelParams[t->relParts[0]].treeAgePr.priorParams[1])
+                        {    
+                        MrBayesPrint("%s   ERROR: The clock rate is inconsistent with the specified tree age prior.\n", spacer);
+                        return (ERROR);
+                        }
+                    }
+                else if (modelParams[t->relParts[0]].treeAgePr.prior == offsetExponential ||
+                         modelParams[t->relParts[0]].treeAgePr.prior == offsetGamma ||
+                         modelParams[t->relParts[0]].treeAgePr.prior == truncatedNormal ||
+                        modelParams[t->relParts[0]].treeAgePr.prior == offsetLogNormal)
+                    {
+                    if (treeAge < modelParams[t->relParts[0]].treeAgePr.priorParams[0])
+                        {    
+                        MrBayesPrint("%s   ERROR: The clock rate is inconsistent with the specified tree age prior.\n", spacer);
+                        return (ERROR);
+                        }
                     }
                 }
             }