Passing timelimit specified on command line to pyomo solver (#97)

* Passing timelimit specified on command line to pyomo solver; bug fix with cd; check-solver in docs

* Added an entry for the check-solver command and an associated doc page

Author: vineetbansal

docs/source/README.md

@@ -137,6 +137,8 @@ section of your `hatchet.ini` (if using the [hatchet run](doc_runhatchet.html) c
 solver. Make sure the relevant solver binaries are in your `$PATH`, otherwise Pyomo will not be able to find them
 correctly.
 
+One HATCHet command that is very useful to sanity-check your solver is `hatchet check-solver`. This command runs the HATCHet `compute_cn` step on a small set of pre-packaged data files and completes fairly quickly (a few seconds for the Gurobi optimizer, but up to a few minutes for glpk). Running this command will ensure that you have your solver settings (including licenses) set up correctly, so you should always run this command first before trying out the `compute_cn` step on your large data files.
+
 ### Required data
 <a name="requireddata"></a>
 

docs/source/doc_check_solver.md

@@ -0,0 +1,19 @@
+# check-solver
+
+This command of HATCHet runs the `compute_cn` step on a set of small data files (.bbc/.seg) pre-packaged with HATCHet, and is a quick way to verify if your solver is working correctly.
+If you are unable to run this command, it likely indicates a licensing issue with default (Gurobi) solver. To use alternative solvers, see the
+[Using a different Pyomo-supported solver](README.html#usingasolver_other) section of the README for more details).
+
+## Input
+
+This command takes no inputs.
+
+## Output
+
+This command produces debugging output from the solver. Look for the message towards the end of the output:
+
+```
+# Your current solver ... seems to be working correctly
+```
+
+to verify if your selected solver works correctly.

docs/source/doc_compute_cn.md

@@ -1,4 +1,4 @@
-# hatchet
+# compute-cn
 
 This step computes allele-specific fractional copy numbers, solves a constrained distance-based simultaneous factorization to compute allele and clone-specific copy numbers and clone proportions, and deploys a model-selection criterion select the number of clone by explicitly considering the trade-off between subclonal copy-number aberrations and whole-genome duplication.
 The step offers some parameters to control each of these features.

docs/source/doc_fullpipeline.md

@@ -64,6 +64,7 @@ hatchet <step>
 | (5)   | [*plot-bins*](doc_plot_bins.html)           | BBot     | This step produces informative plots concerning the computed RDRs, BAFs, and clusters. The information produced by this step are important to validate the compute clusters of genomic regions. |
 | (6)   | [*compute-cn*](doc_compute_cn.html)         | hatchet  | This step computes allele-specific fractional copy numbers, solves a constrained distance-based simultaneous factorization to compute allele and clone-specific copy numbers and clone proportions, and deploys a model-selection criterion select the number of clone by explicitly considering the trade-off between subclonal copy-number aberrations and whole-genome duplication. |
 | (7)   | [*plot-cn*](doc_plot_cn.html)               | BBeval   | This step analyzes the inferred copy-number states and clone proportions and produces informative plots jointly considering all samples from the same patient. In addition, this step can also combine results obtained for different patients and perform integrative analysis. |
+| (8)   | [*check-solver*](doc_check_solver.html)     |          | This step runs the `compute_cn` step on a set of small data files pre-packaged with HATCHet, and is a quick way to verify if your solver is working correctly. |
 
 ## Recommendations and quality control
 <a name="recommendations"></a>

src/hatchet/bin/HATCHet.py

@@ -634,7 +634,8 @@ def execute_python(solver, args, n, outprefix):
         n_seed=args['p'],
         n_worker=args['j'],
         random_seed=args['r'],
-        max_iters=args['f']
+        max_iters=args['f'],
+        timelimit=args['s']
     )
 
     segmentation(cA, cB, u, cluster_ids, sample_ids, bbc_file=args['bbc'], bbc_out_file=bbc_out_file,

src/hatchet/utils/check_solver.py

@@ -27,7 +27,7 @@ def main(args=None):
                 '-p', '5',
                 '-v', '3',
                 '-u', '0.03',
-                '--mode', '1',
+                '--mode', '0',
                 '-r', '6700',
                 '-j', '1',
                 '-eD', '6',

src/hatchet/utils/solve/__init__.py

@@ -19,7 +19,7 @@ def solver_available(solver=None):
 
 
 def solve(clonal, seg_file, n, solver='gurobi', solve_mode='cd', d=-1, cn_max=-1, mu=0.01, diploid_threshold=0.1,
-          ampdel=True, n_seed=400, n_worker=8, random_seed=None, max_iters=None):
+          ampdel=True, n_seed=400, n_worker=8, random_seed=None, max_iters=None, timelimit=None):
 
     assert solve_mode in ('ilp', 'cd', 'both'), 'Unrecognized solve_mode'
     assert solver_available(solver), f'Solver {solver} not available or not licensed'
@@ -65,17 +65,17 @@ def solve(clonal, seg_file, n, solver='gurobi', solve_mode='cd', d=-1, cn_max=-1
     if solve_mode == 'ilp':
         ilp = ILPSubset(n, cn_max, d=d, mu=mu, ampdel=ampdel, copy_numbers=copy_numbers, f_a=f_a, f_b=f_b, w=weights)
         ilp.create_model()
-        return ilp.run(solver_type=solver)
+        return ilp.run(solver_type=solver, timelimit=timelimit)
     elif solve_mode == 'cd':
         cd = CoordinateDescent(f_a=f_a, f_b=f_b, n=n, mu=mu, d=d, cn_max=cn_max, w=weights, ampdel=ampdel,
                                cn=copy_numbers)
-        return cd.run(solver_type=solver, max_iters=max_iters, n_seed=n_seed, j=n_worker, random_seed=random_seed)
+        return cd.run(solver_type=solver, max_iters=max_iters, n_seed=n_seed, j=n_worker, random_seed=random_seed, timelimit=timelimit)
     else:
         cd = CoordinateDescent(f_a=f_a, f_b=f_b, n=n, mu=mu, d=d, cn_max=cn_max, w=weights, ampdel=ampdel,
                                cn=copy_numbers)
-        _, cA, cB, _ = cd.run(solver_type=solver, max_iters=max_iters, n_seed=n_seed, j=n_worker, random_seed=random_seed)
+        _, cA, cB, _, _, _ = cd.run(solver_type=solver, max_iters=max_iters, n_seed=n_seed, j=n_worker, random_seed=random_seed, timelimit=timelimit)
 
         ilp = ILPSubset(n, cn_max, d=d, mu=mu, ampdel=ampdel, copy_numbers=copy_numbers, f_a=f_a, f_b=f_b, w=weights)
         ilp.create_model()
         ilp.hot_start(cA, cB)
-        return ilp.run(solver_type=solver)
+        return ilp.run(solver_type=solver, timelimit=timelimit)

src/hatchet/utils/solve/cd.py

@@ -11,7 +11,7 @@ def __init__(self, ilp, solver):
         self.ilp = ilp
         self.solver_type = solver
 
-    def run(self, cA, cB, u, max_iters, max_convergence_iters, tol=0.001):
+    def run(self, cA, cB, u, max_iters, max_convergence_iters, tol=0.001, timelimit=None):
         _iters = _convergence_iters = 0
         _u = u
         _cA, _cB = cA, cB  # first hot-start values
@@ -22,15 +22,15 @@ def run(self, cA, cB, u, max_iters, max_convergence_iters, tol=0.001):
             carch.fix_u(_u)
             carch.create_model()
             carch.hot_start(_cA, _cB)
-            carch_results = carch.run(self.solver_type)
+            carch_results = carch.run(self.solver_type, timelimit=timelimit)
             if carch_results is None:
                 return None
             _obj_c, _cA, _cB, _, _, _ = carch_results
 
             uarch = copy(self.ilp)
             uarch.fix_c(_cA, _cB)
             uarch.create_model()
-            uarch_results = uarch.run(self.solver_type)
+            uarch_results = uarch.run(self.solver_type, timelimit=timelimit)
             if uarch_results is None:
                 return None
             _obj_u, _, _, _u, _, _ = uarch_results
@@ -47,9 +47,9 @@ def run(self, cA, cB, u, max_iters, max_convergence_iters, tol=0.001):
 
 
 # Top-level 'work' function that can be pickled for multiprocessing
-def _work(cd, u, solver_type, max_iters, max_convergence_iters):
+def _work(cd, u, solver_type, max_iters, max_convergence_iters, timelimit):
     worker = Worker(cd.ilp, solver_type)
-    return worker.run(cd.hcA, cd.hcB, u, max_iters=max_iters, max_convergence_iters=max_convergence_iters)
+    return worker.run(cd.hcA, cd.hcB, u, max_iters=max_iters, max_convergence_iters=max_convergence_iters, timelimit=timelimit)
 
 
 class CoordinateDescent:
@@ -61,15 +61,15 @@ def __init__(self, f_a, f_b, n, mu, d, cn_max, cn, w, ampdel=True):
 
         self.seeds = None
 
-    def run(self, solver_type='gurobi', max_iters=10, max_convergence_iters=2, n_seed=400, j=8, random_seed=None):
+    def run(self, solver_type='gurobi', max_iters=10, max_convergence_iters=2, n_seed=400, j=8, random_seed=None, timelimit=None):
         with Random(random_seed):
             seeds = [self.ilp.build_random_u() for _ in range(n_seed)]
 
         result = {}  # obj. value => (cA, cB, u) mapping
         to_do = []
         with ProcessPoolExecutor(max_workers=min(j, n_seed)) as executor:
             for u in seeds:
-                future = executor.submit(_work, self, u, solver_type, max_iters, max_convergence_iters)
+                future = executor.submit(_work, self, u, solver_type, max_iters, max_convergence_iters, timelimit)
                 to_do.append(future)
 
             for future in as_completed(to_do):

src/hatchet/utils/solve/ilp_subset.py

@@ -538,8 +538,9 @@ def run(self, solver_type='gurobi', timelimit=None, write_path=None):
         else:
             solver = pe.SolverFactory(solver_type)
 
-        timelimit = int(timelimit or config.compute_cn.timelimit or 100)
-        kwargs = {'timelimit': timelimit, 'report_timing': False}
+        kwargs = {'report_timing': False}
+        if timelimit is not None:
+            kwargs['timelimit'] = int(timelimit)
         if solver.warm_start_capable():
             kwargs['warmstart'] = self.warmstart