Merge pull request #1248 from veg/develop

Possible fix for an initial value error in Optimize
veg · Nov 17, 2020 · 41c8a6e · 41c8a6e
2 parents 3ae1dbe + ec43aa2
commit 41c8a6e
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 24 deletions.
diff --git a/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf b/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf
@@ -685,17 +685,27 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod
         for (vi, vv; in; v) {
             values[pnames[vi]] = +vv;
         }
+        if (^"fel.srv"){
+         values[^"fel.alpha.scaler"] = Random (0.75, 1.25);
+        } 
         start.grid + values;
     }
 
-
+    /*
+    GetString (p, ^lf, -1);
+    console.log (p);
+    console.log (values);
+    assert (0);
+    */
+
     Optimize (results, ^lf, {
                 "OPTIMIZATION_METHOD" : "nedler-mead", 
                 "OPTIMIZATION_PRECISION" : 1e-5,
                 "OPTIMIZATION_START_GRID" : start.grid    
             }
            );
-
+
+
     /**
         infer and report ancestral substitutions
     */
@@ -738,20 +748,17 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod
     ref_parameter = (^"fel.scaler_parameter_names")[(^"fel.branches.testable")["VALUEINDEXORDER"][0]];
 
     ^ref_parameter = sum /denominator;
+
 
     if (testable == 1) {
         parameters.SetConstraint ((^"fel.scaler_parameter_names")[^"terms.tree_attributes.background"],ref_parameter, "");
     } else {
-        utility.ForEach (^"fel.branches.testable", "_gname_",
-        '
-            //console.log ("REF " + `&ref_parameter`);
+        for (_gname_; in; ^"fel.branches.testable") {
             _pname_ =  (^"fel.scaler_parameter_names")[_gname_];
-            if (_pname_ != `&ref_parameter`) {
-            	//console.log (_pname_ + "=>" + `&ref_parameter`);
-                parameters.SetConstraint (_pname_,`&ref_parameter`, "");
+            if (_pname_ != ref_parameter) {
+                ^_pname_ := ^ref_parameter;
             }
-        '
-        );
+        }
     }
 
     Optimize (results, ^lf, {"OPTIMIZATION_METHOD" : "nedler-mead"});

diff --git a/src/core/likefunc.cpp b/src/core/likefunc.cpp
@@ -256,8 +256,14 @@ void         BenchmarkThreads (_LikelihoodFunction* lf) {
     lf->SetThreadCount (1);
 #endif
     TimeDifference timer;
-    lf->SetIthIndependent (alterIndex,lf->GetIthIndependent(alterIndex));
-    lf->Compute           ();
+#ifdef  __HYPHYMPI__
+    if (hy_mpi_node_rank == 0)
+#endif
+
+    {
+        lf->SetIthIndependent (alterIndex,lf->GetIthIndependent(alterIndex));
+        lf->Compute           ();
+    }
 
 
 
@@ -276,6 +282,7 @@ hyFloat            tdiff = timer.TimeSinceStart();
             lf->SetThreadCount              (k);
             TimeDifference timer;
             lf->SetIthIndependent           (alterIndex,lf->GetIthIndependent(alterIndex));
+
             lf->Compute                     ();
             tdiff = timer.TimeSinceStart();
             if (tdiff < minDiff) {
@@ -4063,15 +4070,8 @@ _Matrix*        _LikelihoodFunction::Optimize (_AssociativeList const * options)
     SetupLFCaches       ();
     SetupCategoryCaches ();
     computationalResults.Clear();
-    _variables_changed_during_last_compute = new _SimpleList ();
-    variables_changed_during_last_compute = new _AVLList (_variables_changed_during_last_compute);
 
-    for (unsigned long i=0UL; i<indexInd.lLength; i++) {
-        if (GetIthIndependentVar(i)->varFlags & HY_VARIABLE_CHANGED) {
-            variables_changed_during_last_compute->InsertNumber (GetIthIndependentVar(i)->get_index());
-        }
-    }
-
+
 
 #ifdef __HYPHYMPI__
     if (hy_mpi_node_rank == 0) {
@@ -4104,14 +4104,14 @@ _Matrix*        _LikelihoodFunction::Optimize (_AssociativeList const * options)
 #ifdef __HYPHYMPI__
     if (hy_mpi_node_rank == 0) {
 #endif
-
     BenchmarkThreads (this);
-
-
 #ifdef __HYPHYMPI__
     }
 #endif
 
+    _variables_changed_during_last_compute = new _SimpleList ();
+    variables_changed_during_last_compute = new _AVLList (_variables_changed_during_last_compute);
+
 
     bool            skipCG                  = ! CheckEqual (get_optimization_setting (kSkipConjugateGradient, 0.0), 0.0),
                     keepStartingPoint       = ! CheckEqual (get_optimization_setting (kUseLastResults, 0.0), 0.0),
@@ -4133,7 +4133,7 @@ _Matrix*        _LikelihoodFunction::Optimize (_AssociativeList const * options)
     if (keepStartingPoint) {
         indexInd.Each ([this, &untag] (long v, unsigned long i) -> void {
             _Variable *iv = GetIthIndependentVar (i);
-            if (iv->HasBeenInitialized()) {
+            if (!iv->IsGlobal() && iv->HasBeenInitialized()) {
                 if (!iv->IsModified()) {
                     iv->MarkModified();
                     untag << i;
@@ -4178,6 +4178,13 @@ _Matrix*        _LikelihoodFunction::Optimize (_AssociativeList const * options)
             }
         }
     }
+
+    for (unsigned long i=0UL; i<indexInd.lLength; i++) {
+        if (GetIthIndependentVar(i)->varFlags & HY_VARIABLE_CHANGED) {
+            variables_changed_during_last_compute->InsertNumber (GetIthIndependentVar(i)->get_index());
+        }
+    }
+
 
 #if !defined __UNIX__ || defined __HEADLESS__
     SetStatusBarValue (5,maxSoFar,(likeFuncEvalCallCount-evalsIn)/timer.TimeSinceStart());