Pyomo solving enhancements (#107)

* Ability to run compute-cn when building with HATCHET_BUILD_NOEXT=1
* Some enhancements for pyomo solving
* more atomic and faster unit tests

Author: vineetbansal

.github/workflows/main.yml

@@ -98,6 +98,14 @@ jobs:
           PLATFORM: ${{ matrix.python }}
           CXXFLAGS: -pthread
 
+      # We may end up with files that have a ':' in their filenames,
+      # due to the <chr>:start:end notation of chromosomes we process through SAMtools/BCFtools.
+      # upload-artifact is unable to handle these, so we replace colons with dashes
+      - name: Clean up artifact filenames
+        run: |
+          sudo apt install rename
+          find tests/out -name "*:*" -exec rename 's|:|-|g' {} \;
+ 
       - name: Save Pytest Output Data
         uses: actions/upload-artifact@v2
         with:

setup.py

@@ -76,7 +76,7 @@ def build_extension(self, ext):
 
 setup(
     name='hatchet',
-    version='0.4.8',
+    version='0.4.9',
     packages=['hatchet', 'hatchet.utils', 'hatchet.utils.solve', 'hatchet.bin', 'hatchet.data'],
     package_dir={'': 'src'},
     package_data={'hatchet': ['hatchet.ini'], 'hatchet.data': ['*']},

src/hatchet/__init__.py

@@ -1,4 +1,4 @@
-__version__ = '0.4.8'
+__version__ = '0.4.9'
 
 import os.path
 from importlib.resources import path

src/hatchet/utils/combine_counts.py

@@ -47,7 +47,7 @@ def main(args=None):
     nonzerobaf = (lambda rk : all(sample[4] + sample[5] > 0 for sample in rk))
     result = {key : result[key] for key in result if nonzerobaf(result[key])}
     for key in sorted(result, key=(lambda x : (sp.numericOrder(x[0]), int(x[1]), int(x[2])))):
-        for sample in result[key]:
+        for sample in sorted(result[key]):
             sys.stdout.write("{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n".format(key[0], key[1], key[2], sample[0], sample[1], sample[2], sample[3], sample[4], sample[5], sample[6]))
 
 

src/hatchet/utils/solve/__init__.py

@@ -64,7 +64,7 @@ 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()
+        ilp.create_model(pprint=True)
         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,

src/hatchet/utils/solve/cd.py

@@ -57,6 +57,10 @@ class CoordinateDescent:
     def __init__(self, f_a, f_b, n, mu, d, cn_max, cn, w, ampdel=True):
         # ilp attribute used here as a convenient storage container for properties
         self.ilp = ILPSubset(n=n, cn_max=cn_max, d=d, mu=mu, ampdel=ampdel, copy_numbers=cn, f_a=f_a, f_b=f_b, w=w)
+        # Building the model here is not strictly necessary, as, during execution,
+        #   self.carch and c.uarch will copy self.ilp and create+run those models.
+        # However, we do so here simply so we can print out some diagnostic information once for the user.
+        self.ilp.create_model(pprint=True)
         self.hcA, self.hcB = self.ilp.first_hot_start()
 
         self.seeds = None

src/hatchet/utils/solve/ilp_subset.py

@@ -1,3 +1,4 @@
+import textwrap
 import math
 import numpy as np
 import pandas as pd
@@ -55,6 +56,21 @@ def __copy__(self):
             f_a=self.f_a, f_b=self.f_b, w=self.w
         )
 
+    def __str__(self):
+        # Pyomo pprint gives us too much information - too unwieldy for large models
+        # This method is implemented to supply the bare-minimum but useful model information.
+        if self.model is None:
+            return ''
+        else:
+            return textwrap.dedent(f"""
+                # ------------------------------------------
+                #   Problem Information
+                # ------------------------------------------
+                #     Number of constraints: {self.model.nconstraints()}
+                #     Number of variables: {self.model.nvariables()}
+                # ------------------------------------------
+            """)
+
     @property
     def M(self):
         return math.floor(math.log2(self.cn_max)) + 1
@@ -88,7 +104,7 @@ def optimized_u(self):
         else:
             return self.u
 
-    def create_model(self):
+    def create_model(self, pprint=False):
 
         m, n, k = self.m, self.n, self.k
         f_a, f_b = self.f_a, self.f_b
@@ -370,6 +386,9 @@ def create_model(self):
         model.obj = pe.Objective(expr=obj, sense=pe.minimize)
         self.model = model
 
+        if pprint:
+            print(str(self))
+
     def build_symmetry_breaking(self, model):
         for i in range(1, self.n - 1):
             _sum1 = 0

src/hatchet/utils/solve/utils.py

@@ -38,14 +38,17 @@ def parse_clonal(clonal):
         cluster_id, cn_a, cn_b = [int(_c) for _c in c.split(':')]
 
         # The first two clonal clusters (used for scaling) MUST have different total copy numbers
-        if i==1:
+        if i == 1:
             _first_cn_a, _first_cn_b = list(copy_numbers.values())[0]
             if _first_cn_a + _first_cn_b == cn_a + cn_b:
                 raise ValueError('When >= 2 clonal copy numbers are given, the first two must be different in the two segmental clusters')
 
         cn_total = cn_a + cn_b
         if (cn_total == 2) and (n_clonal_parts > 1):
-            warnings.warn('Please specify a single cluster when CN_A+CN_B=2')
+            # warnings.warn('Please specify a single cluster when CN_A+CN_B=2')
+            # TODO: The C++ implementation generates a warning corresponding to the above
+            # This is suppressed by default, which is what we do here for now.
+            pass
         if cluster_id in copy_numbers:
             raise ValueError('Already encountered cluster_id =', str(cluster_id))
         copy_numbers[cluster_id] = cn_a, cn_b