Re-implement open source solvers

environment.yml

@@ -11,4 +11,7 @@ dependencies:
   - pyparsing>=2.2.0
   - cvxopt
   - ibmdecisionoptimization::cplex
+  - osqp
+  - pip:
+    - "git+https://github.com/osqp/miosqp.git@ac672338b0593d865dd15b7a76434f25e24244a9#egg=miosqp"
   - tqdm

setup.py

@@ -44,6 +44,10 @@ def main():
         ext_modules=ext_modules,
         setup_requires=setup_requires,
         install_requires=install_requires,
+        extra_dependencies={
+            "cplex": ["cvxopt", "cplex"],
+            "osqp": ["osqp", "miosqp @ git+https://github.com/osqp/miosqp.git@ac672338b0593d865dd15b7a76434f25e24244a9#egg=miosqp"],
+        },
         zip_safe=False,
         python_requires=">=3.9",
         entry_points={

src/qfit/qfit.py

@@ -14,7 +14,7 @@
 from .samplers import ChiRotator, CBAngleRotator, BondRotator
 from .samplers import CovalentBondRotator, GlobalRotator
 from .samplers import RotationSets, Translator
-from .solvers import QPSolver, MIQPSolver, SolverError
+from .solvers import SolverError, get_qp_solver_class, get_miqp_solver_class
 from .structure import Structure, _Segment, calc_rmsd
 from .structure.residue import residue_type
 from .structure.ligand import BondOrder
@@ -30,12 +30,7 @@
 
 # Create a namedtuple 'class' (struct) which carries info about an MIQP solution
 MIQPSolutionStats = namedtuple(
-    "MIQPSolutionStats", ["threshold", "BIC", "rss", "objective", "weights"]
-)
-
-# Create a namedtuple 'class' (struct) which carries info about an MIQP solution
-MIQPSolutionStats = namedtuple(
-    "MIQPSolutionStats", ["threshold", "BIC", "rss", "objective", "weights"]
+    "MIQPSolutionStats", ["threshold", "BIC", "rss", "objective_value", "weights"]
 )
 
 
@@ -78,7 +73,8 @@ def __init__(self):
         self.rmsd_cutoff = 0.01
 
         # MIQP options
-        self.cplex = True
+        self.qp_solver = None
+        self.miqp_solver = None
         self.cardinality = 5
         self.threshold = 0.20
         self.bic_threshold = True
@@ -277,14 +273,15 @@ def _convert(self):
     def _solve_qp(self):
         # Create and run solver
         logger.info("Solving QP")
-        solver = QPSolver(self._target, self._models, use_cplex=self.options.cplex)
-        solver.solve()
+        qp_solver_class = get_qp_solver_class(self.options.qp_solver)
+        solver = qp_solver_class(self._target, self._models)
+        solver.solve_qp()
 
         # Update occupancies from solver weights
         self._occupancies = solver.weights
 
         # Return solver's objective value (|ρ_obs - Σ(ω ρ_calc)|)
-        return solver.obj_value
+        return solver.objective_value
 
     def _solve_miqp(
         self,
@@ -303,16 +300,17 @@ def _solve_miqp(
 
         # Create solver
         logger.info("Solving MIQP")
-        solver = MIQPSolver(self._target, self._models, use_cplex=self.options.cplex)
+        miqp_solver_class = get_miqp_solver_class(self.options.miqp_solver)
+        solver = miqp_solver_class(self._target, self._models)
 
         # Threshold selection by BIC:
         if do_BIC_selection:
             # Iteratively test decreasing values of the threshold parameter tdmin (threshold)
             # to determine if the better fit (RSS) justifies the use of a more complex model (k)
             miqp_solutions = []
             for threshold in loop_range:
-                solver.solve(cardinality=None, threshold=threshold)
-                rss = solver.obj_value * self._voxel_volume
+                solver.solve_miqp(cardinality=None, threshold=threshold)
+                rss = solver.objective_value * self._voxel_volume
                 n = len(self._target)
 
                 natoms = self._coor_set[0].shape[0]
@@ -328,7 +326,7 @@ def _solve_miqp(
                     threshold=threshold,
                     BIC=BIC,
                     rss=rss,
-                    objective=solver.obj_value.copy(),
+                    objective_value=solver.objective_value.copy(),
                     weights=solver.weights.copy(),
                 )
                 miqp_solutions.append(solution)
@@ -338,17 +336,17 @@ def _solve_miqp(
             self._occupancies = miqp_solution_lowest_bic.weights
 
             # Return solver's objective value (|ρ_obs - Σ(ω ρ_calc)|)
-            return miqp_solution_lowest_bic.objective
+            return miqp_solution_lowest_bic.objective_value
 
         else:
             # Run solver with specified parameters
-            solver.solve(cardinality=cardinality, threshold=threshold)
+            solver.solve_miqp(cardinality=cardinality, threshold=threshold)
 
             # Update occupancies from solver weights
             self._occupancies = solver.weights
 
             # Return solver's objective value (|ρ_obs - Σ(ω ρ_calc)|)
-            return solver.obj_value
+            return solver.objective_value
 
     def sample_b(self):
         """Create copies of conformers that vary in B-factor.

src/qfit/qfit_covalent_ligand.py

@@ -14,6 +14,7 @@
 
 from . import MapScaler, Structure, XMap, Covalent_Ligand
 from . import QFitCovalentLigand, QFitOptions
+from .solvers import available_qp_solvers, available_miqp_solvers
 
 os.environ["OMP_NUM_THREADS"] = "1"
 
@@ -290,6 +291,22 @@ def parse_args():
         help="Do not use BIC to select the most parsimonious MIQP threshold",
     )
 
+    # Solver options
+    p.add_argument(
+        "--qp-solver",
+        dest="qp_solver_choice",
+        choices=available_qp_solvers.keys(),
+        default=next(iter(available_qp_solvers.keys())),
+        help="Select the QP solver",
+    )
+    p.add_argument(
+        "--miqp-solver",
+        dest="miqp_solver_choice",
+        choices=available_miqp_solvers.keys(),
+        default=next(iter(available_miqp_solvers.keys())),
+        help="Select the MIQP solver",
+    )
+
     # Output options
     p.add_argument(
         "-d",

src/qfit/qfit_ligand.py

@@ -12,6 +12,8 @@
 from . import MapScaler, Structure, XMap, _Ligand
 from .qfit import QFitLigand, QFitOptions
 from .logtools import setup_logging, log_run_info
+from .solvers import available_qp_solvers, available_miqp_solvers
+
 
 logger = logging.getLogger(__name__)
 os.environ["OMP_NUM_THREADS"] = "1"
@@ -249,6 +251,22 @@ def build_argparser():
         help="Use BIC to select the most parsimonious MIQP threshold",
     )
 
+    # Solver options
+    p.add_argument(
+        "--qp-solver",
+        dest="qp_solver_choice",
+        choices=available_qp_solvers.keys(),
+        default=next(iter(available_qp_solvers.keys())),
+        help="Select the QP solver",
+    )
+    p.add_argument(
+        "--miqp-solver",
+        dest="miqp_solver_choice",
+        choices=available_miqp_solvers.keys(),
+        default=next(iter(available_miqp_solvers.keys())),
+        help="Select the MIQP solver",
+    )
+
     # Output options
     p.add_argument(
         "-d",

src/qfit/qfit_protein.py

@@ -26,6 +26,7 @@
     QueueListener,
 )
 from . import MapScaler, Structure, XMap
+from .solvers import available_qp_solvers, available_miqp_solvers
 from .structure.rotamers import ROTAMERS
 
 
@@ -325,6 +326,22 @@ def build_argparser():
         help="Number of processors to use",
     )
 
+    # Solver options
+    p.add_argument(
+        "--qp-solver",
+        dest="qp_solver",
+        choices=available_qp_solvers.keys(),
+        default=next(iter(available_qp_solvers.keys())),
+        help="Select the QP solver",
+    )
+    p.add_argument(
+        "--miqp-solver",
+        dest="miqp_solver",
+        choices=available_miqp_solvers.keys(),
+        default=next(iter(available_miqp_solvers.keys())),
+        help="Select the MIQP solver",
+    )
+
     # qFit Segment options
     p.add_argument(
         "--only-segment",

src/qfit/qfit_residue.py

@@ -17,8 +17,10 @@
 from . import MapScaler, Structure, XMap
 from .qfit import QFitOptions
 from . import QFitRotamericResidue
+from .solvers import available_qp_solvers, available_miqp_solvers
 from .structure import residue_type
 
+
 logger = logging.getLogger(__name__)
 os.environ["OMP_NUM_THREADS"] = "1"
 
@@ -321,6 +323,22 @@ def build_argparser():
         help="Use BIC to select the most parsimonious MIQP threshold",
     )
 
+    # Solver options
+    p.add_argument(
+        "--qp-solver",
+        dest="qp_solver_choice",
+        choices=available_qp_solvers.keys(),
+        default=next(iter(available_qp_solvers.keys())),
+        help="Select the QP solver",
+    )
+    p.add_argument(
+        "--miqp-solver",
+        dest="miqp_solver_choice",
+        choices=available_miqp_solvers.keys(),
+        default=next(iter(available_miqp_solvers.keys())),
+        help="Select the MIQP solver",
+    )
+
     # Output options
     p.add_argument(
         "-d",