Merge pull request #21 from raphael-group/gc-fix-2

Adding functional tests and updating flag behavior

Examples.md

@@ -1,5 +1,5 @@
 # Prerequisites
-To try the following examples, first do the followings:
+To try the following examples, first do the following:
 1. Download the data from [examples.zip](https://github.com/raphael-group/laml/tree/master/examples.zip)
 2. Unzip the downloaded file. After unzipping, you should see a folder named ``examples``
 3. Change directory to ``examples``

README.md

@@ -70,7 +70,7 @@ There are four output files:
 4. `LAML_output.log`: The LAML logfile.
 
 ## Examples
-We provide two examples for the two common use cases of LAML. See [Examples.md](https://github.com/raphael-group/laml/tree/master/Examples.md) for more details.
+We provide two examples for two common LAML use cases. See [Examples.md](https://github.com/raphael-group/laml/tree/master/Examples.md) for more details.
 
 ## Advanced I/O options
 LAML has the following additional options for I/O
@@ -91,7 +91,7 @@ Note: LAML also accepts a character matrix that contains negative integers and/o
 
 ### The mutation priors
 While not strictly required, **mutation priors** can have a large effect on the outputs. If no mutation priors are provided, LAML uses uniform priors
-by default. However, if possible we highly recommend specifying mutation prior using `-p`. We accept the following two formats for mutation priors 
+by default. However, if possible we highly recommend specifying mutation prior using `-p`. We accept the following two formats for mutation priors:
 
 **Recommended** A file containing the prior matrix, a [comma-separated values (CSV) file](https://en.wikipedia.org/wiki/Comma-separated_values), with three columns: site index, character state, and probability. The site index and character states must be integers, and the probability must be a float. We do *not* expect the unmutated state to appear in the alphabet. See an example input prior file in
 [examples/example1/priors.csv](https://github.com/raphael-group/LAML/tree/laml/examples/example1/priors.csv).

extra_laml_tests.py

@@ -1,15 +1,16 @@
-from unit_tests.unit_tests_MLSolver import *
-from unit_tests.unit_tests_Simulator import *
-#from unit_tests.unit_tests_SpaLinSolver import *
-#from unit_tests.unit_tests_TopoSearch import *
-#from unit_tests.unit_tests_TopoSearchParallel import *
-#from unit_tests.unit_tests_EMSolver import *
-#from unit_tests.utils import *
+#! /usr/bin/env python
+##from laml_unit_tests.unit_tests_MLSolver import *
+##from laml_unit_tests.unit_tests_Simulator import *
+from laml_unit_tests.unit_tests_TopoSearch import *
+from laml_unit_tests.unit_tests_TopoSearchParallel import *
+from laml_unit_tests.unit_tests_io import *
 import sys
 import os
 
-def main():
-    sys.path.append(os.path.join(os.path.dirname(__file__), 'helpers'))
-
 if __name__ == '__main__':
+
+    sys.path.append(os.path.dirname(__file__))
+    sys.path.append(os.path.join(os.path.dirname(__file__), 'helpers'))
+    print("Running extra tests for LAML...")
     unittest.main()
+

laml_libs/EM_solver.py

@@ -401,35 +401,36 @@ def __optimize_nu__(d): # d is a vector of all branch lengths
 
         nIters = 1
         nu_star = self.params.nu
-        if (nu_star == 0):
-            d_star,status_d = __compute_brlen__()
-            status_nu = "optimal"
-        else:
-            for r in range(nIters):
+        # [deprecated]: will always scale to timescale! cannot use closed-form solution
+        #if (nu_star == 0): 
+        #    d_star,status_d = __compute_brlen__()
+        #    status_nu = "optimal"
+        #else:
+        for r in range(nIters):
+            if verbose > 0:
+                print("Optimizing branch lengths. Current phi: " + str(phi_star) + ". Current nu:" + str(nu_star))
+            try:
+                d_star,status_d = __optimize_brlen__(nu_star,verbose=False)
+            except:
+                d_star = d_ini
+                status_d = "failure"
+            #d_star = np.array([max(x,self.dmin) for x in d_star])     
+            if status_d == "infeasible": # should only happen with local EM 
+                return False,"d_infeasible"
+            if not optimize_nu:
+                if verbose > 0:
+                    print("Fixing nu to " + str(self.params.nu))
+                nu_star = self.params.nu
+                status_nu = "optimal"    
+            else:    
                 if verbose > 0:
-                    print("Optimizing branch lengths. Current phi: " + str(phi_star) + ". Current nu:" + str(nu_star))
+                    print("Optimizing nu")
                 try:
-                    d_star,status_d = __optimize_brlen__(nu_star,verbose=False)
+                    nu_star,status_nu = __optimize_nu__(d_star)                 
                 except:
-                    d_star = d_ini
-                    status_d = "failure"
-                #d_star = np.array([max(x,self.dmin) for x in d_star])     
-                if status_d == "infeasible": # should only happen with local EM 
-                    return False,"d_infeasible"
-                if not optimize_nu:
-                    if verbose > 0:
-                        print("Fixing nu to " + str(self.params.nu))
-                    nu_star = self.params.nu
-                    status_nu = "optimal"    
-                else:    
-                    if verbose > 0:
-                        print("Optimizing nu")
-                    try:
-                        nu_star,status_nu = __optimize_nu__(d_star)                 
-                    except:
-                        status_nu = "failure"
-                    #if status_nu != "optimal":
-                    #    return False,status_nu
+                    status_nu = "failure"
+                #if status_nu != "optimal":
+                #    return False,status_nu
         # place the optimal value back to params
         self.params.phi = phi_star
         self.params.nu = nu_star

laml_unit_tests/test_data/test_inputs/example1/README

@@ -0,0 +1 @@
+This is s50d50p01_sub250_r01.