Repository: ExcitedStates/qfit-3.0

Issue:
Install & Test open source optimization package 
To get qFit into SBGrid, we must find another solution to IBM's CPLEX (this is not open source). I suggest using https://github.com/cvxgrp/cvxpy. It seems to be the most robust of the packages, but I would love to get feedback from @blake-riley or @natechols, or anyone else. 

Discussion:

Patch:
diff --git a/environment.yml b/environment.yml
index 80a9d00..b9774c4 100644
--- a/environment.yml
+++ b/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
diff --git a/setup.py b/setup.py
index 5b10350..3f54dc0 100644
--- a/setup.py
+++ b/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={
diff --git a/src/qfit/qfit.py b/src/qfit/qfit.py
index 323cc7b..48c16be 100644
--- a/src/qfit/qfit.py
+++ b/src/qfit/qfit.py
@@ -14,7 +14,7 @@ from .clash import ClashDetector
 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 @@ logger = logging.getLogger(__name__)
 
 # 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 @@ class QFitOptions:
         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 @@ class _BaseQFit:
     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,7 +300,8 @@ class _BaseQFit:
 
         # 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:
@@ -311,8 +309,8 @@ class _BaseQFit:
             # 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 @@ class _BaseQFit:
                     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 @@ class _BaseQFit:
             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.
diff --git a/src/qfit/qfit_covalent_ligand.py b/src/qfit/qfit_covalent_ligand.py
index 7f4f200..b1ca589 100644
--- a/src/qfit/qfit_covalent_ligand.py
+++ b/src/qfit/qfit_covalent_ligand.py
@@ -14,6 +14,7 @@ import numpy as np
 
 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",
diff --git a/src/qfit/qfit_ligand.py b/src/qfit/qfit_ligand.py
index c34b2b7..23d028f 100644
--- a/src/qfit/qfit_ligand.py
+++ b/src/qfit/qfit_ligand.py
@@ -12,6 +12,8 @@ from string import ascii_uppercase
 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",
diff --git a/src/qfit/qfit_protein.py b/src/qfit/qfit_protein.py
index db96443..95b29b2 100644
--- a/src/qfit/qfit_protein.py
+++ b/src/qfit/qfit_protein.py
@@ -26,6 +26,7 @@ from .logtools import (
     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",
diff --git a/src/qfit/qfit_residue.py b/src/qfit/qfit_residue.py
index 4cfbf57..381e5d1 100644
--- a/src/qfit/qfit_residue.py
+++ b/src/qfit/qfit_residue.py
@@ -17,8 +17,10 @@ from .logtools import setup_logging, log_run_info
 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",
diff --git a/src/qfit/solvers.py b/src/qfit/solvers.py
index 2e053db..279bcc4 100644
--- a/src/qfit/solvers.py
+++ b/src/qfit/solvers.py
@@ -1,282 +1,789 @@
+from __future__ import annotations
+
+import inspect
 import logging
+import sys
+from abc import ABC, abstractmethod
+from types import ModuleType
+from typing import TYPE_CHECKING, Any, Optional, cast
 
 import numpy as np
-from scipy import sparse
+import scipy as sci
+import scipy.sparse  # pylint: disable=unused-import
+from numpy.typing import NDArray
 
+from .utils.optional_lazy_import import lazy_load_module_if_available
 
 logger = logging.getLogger(__name__)
 
-
-# Try load solver sets
-try:
-    import cvxopt
-    import cplex
-except ImportError:
-    CPLEX = False
+SolverError: tuple[type[Exception], ...] = (RuntimeError,)
+
+__all__ = [
+    "available_qp_solvers",
+    "available_miqp_solvers",
+    "get_qp_solver_class",
+    "get_miqp_solver_class",
+]
+
+
+###############################
+# Define solver "interfaces" / Abstractions
+#   All solvers in this module must subclass either QPSolver or MIQPSolver, and conform to the interface.
+#   The functions in the "Helper Methods" section depend on the subclass membership of these ABCs.
+###############################
+
+
+class GenericSolver(ABC):
+    # Class variables
+    driver_pkg_name: str
+    driver: Optional[ModuleType]
+
+    # Instance variables
+    target: NDArray[np.float_]  # 1D array of shape (n_voxels,)
+    models: NDArray[np.float_]  # 2D array of shape (n_models, n_voxels,)
+
+    weights: Optional[NDArray[np.float_]] = None
+    objective_value: Optional[float] = None
+
+
+class QPSolver(GenericSolver):
+    """Finds the combination of conformer-occupancies that minimizes difference density.
+
+    Problem statement
+    -----------------
+    We have observed density ρ^o from the user-provided map (target).
+    We also have a set of conformers, each with modelled/calculated density ρ^c_i.
+    We want find the vector of occupancies ω = <ω_0, ..., ω_n> that minimizes
+        the difference between the observed and modelled density --- that minimizes
+        a residual sum-of-squares function, rss(ω).
+    Mathematically, we wish to minimize:
+        min_ω rss(ω) = min_ω || ρ^c ω - ρ^o ||^2
+
+    Expanding & rearranging rss(ω):
+        rss(ω) = ( ρ^c ω - ρ^o ).T  ( ρ^c ω - ρ^o )
+                = ω.T ρ^c.T ρ^c ω - 2 ρ^o.T ρ^c ω + ρ^o.T ρ^o
+    We can rewrite this as
+        rss(ω) = ω.T P ω + 2 q.T ω + C
+    where
+        P =  ρ^c.T ρ^c
+        q = -ρ^c.T ρ^o
+        C =  ρ^o.T ρ^o
+
+    Noting that the canonical QP objective function is
+        g(x) = 1/2 x.T P x + q.T x
+    we can use a QP solver to find min_x g(x), which, by equivalence,
+        will provide the solution to min_ω rss(ω).
+
+    Solution constraints
+    --------------------
+    Furthermore, these occupancies are meaningful parameters, so we require
+    that their sum is within the unit interval:
+        Σ ω_i ≤ 1
+    and that each individual occupancy is a positive fractional number:
+        0 ≤ ω_i ≤ 1
+    """
+
+    @abstractmethod
+    def solve_qp(self) -> None:
+        ...
+
+
+class MIQPSolver(GenericSolver):
+    """Finds the combination of conformer-occupancies that minimizes difference density.
+
+    Problem statement
+    -----------------
+    We have observed density ρ^o from the user-provided map (target).
+    We also have a set of conformers, each with modelled/calculated density ρ^c_i.
+    We want find the vector of occupancies ω = <ω_0, ..., ω_n> that minimizes
+        the difference between the observed and modelled density --- that minimizes
+        a residual sum-of-squares function, rss(ω).
+    Mathematically, we wish to minimize:
+        min_ω rss(ω) = min_ω || ρ^c ω - ρ^o ||^2
+
+    Expanding & rearranging rss(ω):
+        rss(ω) = ( ρ^c ω - ρ^o ).T  ( ρ^c ω - ρ^o )
+                = ω.T ρ^c.T ρ^c ω - 2 ρ^o.T ρ^c ω + ρ^o.T ρ^o
+    We can rewrite this as
+        rss(ω) = ω.T P ω + 2 q.T ω + C
+    where
+        P =  ρ^c.T ρ^c
+        q = -ρ^c.T ρ^o
+        C =  ρ^o.T ρ^o
+
+    Noting that the canonical QP objective function is
+        g(x) = 1/2 x.T P x + q.T x
+    we can use a QP solver to find min_x g(x), which, by equivalence,
+        will provide the solution to min_ω rss(ω).
+
+    Solution constraints
+    --------------------
+    Furthermore, these occupancies are meaningful parameters, so we require
+    that their sum is within the unit interval:
+        Σ ω_i ≤ 1
+
+    We also want to have either:
+        (a) a set of conformers of known size (cardinality), or
+        (b) a set of conformers with _at least_ threshold occupancy, or else zero (threshold).
+    This can be achieved with a mixed-integer linear constraint:
+        z_i t_min ≤ ω_i ≤ z_i
+    where
+        z_i ∈ {0, 1}
+        t_min is the minimum-allowable threshold value for ω.
+    """
+
+    @abstractmethod
+    def solve_miqp(
+        self,
+        threshold: Optional[float] = None,
+        cardinality: Optional[int] = None,
+        exact: bool = False,
+    ) -> None:
+        ...
+
+
+###############################
+# Define solver implementations
+###############################
+
+
+class CVXOPTSolver(QPSolver):
+    driver_pkg_name = "cvxopt"
+    driver = lazy_load_module_if_available(driver_pkg_name)
+    if TYPE_CHECKING:
+        import cvxopt
+
+    def __init__(self, target: NDArray[np.float_], models: NDArray[np.float_]) -> None:
+        if TYPE_CHECKING:
+            self.driver = self.cvxopt
+            assert self.driver is not None
+
+        # Initialize variables
+        self.target = target
+        self.models = models
+
+        # self.quad_obj: self.driver.matrix
+        # self.lin_obj: self.driver.matrix
+        # self.le_constraints: self.driver.spmatrix
+        # self.le_bounds: self.driver.matrix
+
+        self.weights: Optional[NDArray[np.float_]] = None
+        self.objective_value: Optional[float] = None
+
+        self.solution: dict[str, Any]
+        self.nconformers = models.shape[0]
+
+        # Get the driver & set options
+        self.driver.solvers.options["show_progress"] = False
+        self.driver.solvers.options["abstol"] = 1e-8
+        self.driver.solvers.options["reltol"] = 1e-7
+        self.driver.solvers.options["feastol"] = 1e-8
+
+    def solve_qp(self) -> None:
+        if TYPE_CHECKING:
+            self.driver = self.cvxopt
+            assert self.driver is not None
+
+        # Set up the matrices and restraints
+        logger.debug(
+            "Building cvxopt matrix, size: (%i,%i)",
+            self.nconformers,
+            self.nconformers,
+        )
+
+        # minimize 0.5 x.T P x + q.T x
+        #   where P = self.quad_obj        =   ρ_model.T ρ_model
+        #         q = self.lin_obj         = - ρ_model.T ρ_obs
+        self.quad_obj = self.driver.matrix(np.inner(self.models, self.models), tc="d")
+        self.lin_obj = self.driver.matrix(-np.inner(self.models, self.target), tc="d")
+
+        # subject to:
+        #   G x ≤ h
+        #   where G = self.le_constraints
+        #         h = self.le_bounds
+        #
+        #   Each weight x falls in the closed interval [0..1], and the sum of all weights is <= 1.
+        #   This corresponds to (2 * nconformers + 1) constraints, imposed on (nconformers) variables:
+        #     the lower bound accounts for (nconformers) constraints,
+        #     the upper bound accounts for (nconformers) constraints,
+        #     the summation constraint accounts for 1.
+        #   We construct G with a sparse matrix.
+        #         constraint idx=[0..N)           constraint idx=[N..2N)                             constraint idx=2N
+        rowidxs = [*range(0, self.nconformers)] + [*range(self.nconformers, 2 * self.nconformers)] + (self.nconformers * [2 * self.nconformers])  # fmt:skip
+        #         -x_i ≤ 0                        x_i ≤ 1                                            Σ_i x_i ≤ 1
+        coeffs  = (self.nconformers * [-1.0])   + (self.nconformers * [1.0])                       + (self.nconformers * [1.0])  # fmt:skip
+        colidxs = [*range(0, self.nconformers)] + [*range(0, self.nconformers)]                    + [*range(0, self.nconformers)]  # fmt:skip
+        threshs = (self.nconformers * [0.0])    + (self.nconformers * [1.0])                       + [1.0]  # fmt:skip
+        self.le_constraints = self.driver.spmatrix(coeffs, rowidxs, colidxs, tc="d")
+        self.le_bounds = self.driver.matrix(threshs, tc="d")
+
+        # Solve
+        self.solution = self.driver.solvers.qp(
+            self.quad_obj, self.lin_obj, self.le_constraints, self.le_bounds
+        )
+
+        # Store the density residual and the weights
+        self.objective_value = (
+            2 * self.solution["primal objective"] +
+            np.inner(self.target, self.target)
+        )  # fmt:skip
+        self.weights = np.asarray(self.solution["x"]).ravel()
+
+
+class CPLEXSolver(QPSolver, MIQPSolver):
+    driver_pkg_name = "cplex"
+    driver = lazy_load_module_if_available(driver_pkg_name)
+    if TYPE_CHECKING:
+        import cplex
+
+    def __init__(
+        self, target: NDArray[np.float_], models: NDArray[np.float_], nthreads: int = 1
+    ) -> None:
+        if TYPE_CHECKING:
+            self.driver = self.cplex
+            assert self.driver is not None
+
+        # Initialize variables
+        self.target = target
+        self.models = models
+
+        self.quad_obj: list[object] = []
+        self.lin_obj: list[tuple[int, float]] = []
+
+        self.weights: Optional[NDArray[np.float_]] = None
+        self.objective_value: Optional[float] = None
+
+        self.nconformers = models.shape[0]
+        self.nthreads = nthreads
+
+        # Get the driver & append raisable Exceptions to SolverError class in module (global) scope
+        CplexSolverError: type[Exception] = self.driver.exceptions.CplexSolverError
+        global SolverError
+        SolverError += (CplexSolverError,)
+
+    def compute_quadratic_coeffs(self) -> None:
+        """Precompute the quadratic coefficients (P, q).
+
+        These values don't depend on threshold/cardinality.
+        Having these objectives pre-computed saves a few cycles when MIQP
+        is evaluated for multiple values of threshold/cardinality.
+        """
+        if TYPE_CHECKING:
+            self.driver = self.cplex
+            assert self.driver is not None
+
+        # minimize 0.5 x.T P x + q.T x
+        #   where P = self.quad_obj =   ρ_model.T ρ_model
+        #         q = self.lin_obj  = - ρ_model.T ρ_obs
+        quad_obj_coeffs = np.inner(self.models, self.models)
+        lin_obj_coeffs = -np.inner(self.models, self.target)
+
+        # We have to unpack the arrays for CPLEX into CSR format (even though they're dense)
+        self.quad_obj = []
+        for row in quad_obj_coeffs:
+            idxs, vals = zip(*enumerate(row.tolist()))
+            self.quad_obj.append(self.driver.SparsePair(ind=idxs, val=vals))
+
+        # CPLEX requires linear objectives as a list of (idx, val) tuples
+        self.lin_obj = list(enumerate(lin_obj_coeffs))
+
+    def solve_miqp(
+        self,
+        threshold: Optional[float] = None,
+        cardinality: Optional[int] = None,
+        exact: bool = False,
+    ) -> None:
+        if TYPE_CHECKING:
+            self.driver = self.cplex
+            assert self.driver is not None
+
+        if not (self.quad_obj and self.lin_obj):
+            self.compute_quadratic_coeffs()
+
+        # Create and configure the cplex object
+        miqp = self.driver.Cplex()
+        miqp.set_results_stream(None)
+        miqp.set_log_stream(None)
+        miqp.set_warning_stream(None)
+        miqp.set_error_stream(None)
+        if self.nthreads is not None:
+            miqp.parameters.threads.set(self.nthreads)
+
+        # Create variables and set linear constraints
+        # w_i ≤ 1
+        variable_names = [f"w{n}" for n in range(self.nconformers)]
+        upper_bounds = self.nconformers * [1.0]
+        miqp.variables.add(names=variable_names, ub=upper_bounds)
+
+        # Σ_i w_i ≤ 1
+        ind = [f"w{n}" for n in range(self.nconformers)]
+        val = self.nconformers * [1.0]
+        miqp.linear_constraints.add(
+            lin_expr=[self.driver.SparsePair(ind=ind, val=val)],
+            senses=["L"],
+            rhs=[1.0],
+        )
+
+        # Setup quadratic objective of the MIQP
+        miqp.objective.set_quadratic(self.quad_obj)
+        miqp.objective.set_linear(self.lin_obj)
+        miqp.objective.set_sense(miqp.objective.sense.minimize)
+
+        # If cardinality or threshold is specified, the problem is a MIQP,
+        # so we need to add binary integer variables z_i.
+        if cardinality not in (None, 0) or threshold not in (None, 0):
+            # z_i ∈ {0, 1}
+            integer_names = [f"z{n}" for n in range(self.nconformers)]
+            variable_types = self.nconformers * miqp.variables.type.binary
+            miqp.variables.add(names=integer_names, types=variable_types)
+
+            # Only count weights for which z_i is 1
+            # w_i - z_i ≤ 0
+            # (∵ z_i ∈ {0,1} and 0 ≤ w_i ≤ 1, this is only true when z_i = 1)
+            lin_expr = [
+                self.driver.SparsePair(ind=[f"w{n}", f"z{n}"], val=[1, -1])
+                for n in range(self.nconformers)
+            ]
+            senses = self.nconformers * ["L"]
+            rhs = self.nconformers * [0.0]
+            miqp.linear_constraints.add(lin_expr=lin_expr, senses=senses, rhs=rhs)
+
+        # Set the threshold constraint if applicable
+        if threshold not in (None, 0):
+            # tdmin z_i - w_i ≤ 0, i.e. w_i ≥ tdmin
+            lin_expr = [
+                self.driver.SparsePair(ind=[f"z{n}", f"w{n}"], val=[threshold, -1])
+                for n in range(self.nconformers)
+            ]
+            senses = self.nconformers * ["L"]
+            rhs = self.nconformers * [0.0]
+            miqp.linear_constraints.add(lin_expr=lin_expr, senses=senses, rhs=rhs)
+
+        # Set the cardinality constraint if applicable
+        if cardinality not in (None, 0):
+            # Σ z_i ≤ cardinality
+            cardinality = cast(int, cardinality)  # typing noop
+            ind = [f"z{n}" for n in range(self.nconformers)]
+            val = self.nconformers * [1.0]
+            lin_expr = [self.driver.SparsePair(ind=ind, val=val)]
+            if exact:
+                senses = ["E"]
+                rhs = [min(cardinality, self.nconformers)]
+            else:
+                senses = ["L"]
+                rhs = [cardinality]
+            miqp.linear_constraints.add(lin_expr=lin_expr, senses=senses, rhs=rhs)
+
+        # Solve
+        miqp.solve()
+
+        # Store the density residual and the weights
+        self.objective_value = (
+            2 * miqp.solution.get_objective_value()
+            + np.inner(self.target, self.target)
+        )  # fmt:skip
+        self.weights = np.array(miqp.solution.get_values()[: self.nconformers])
+
+        # Close the cplex object
+        miqp.end()
+
+    def solve_qp(self) -> None:
+        if TYPE_CHECKING:
+            self.driver = self.cplex
+            assert self.driver is not None
+
+        # We can re-use the MIQP code above, provided no cardinality or threshold.
+        self.solve_miqp(threshold=None, cardinality=None, exact=False)
+
+
+class OSQPSolver(QPSolver):
+    driver_pkg_name = "osqp"
+    driver = lazy_load_module_if_available(driver_pkg_name)
+    if TYPE_CHECKING:
+        import osqp
+
+    OSQP_SETTINGS = {
+        "eps_abs": 1e-06,
+        "eps_rel": 1e-06,
+        "eps_prim_inf": 1e-07,
+        "verbose": False,
+    }
+
+    def __init__(self, target: NDArray[np.float_], models: NDArray[np.float_]) -> None:
+        self.target = target
+        self.models = models
+
+        self.nconformers = models.shape[0]
+
+        self.quad_obj: sci.sparse.csc_matrix
+        self.lin_obj: NDArray[np.float_]
+        self.constraints: sci.sparse.csc_matrix
+        self.lower_bounds: NDArray[np.float_]
+        self.upper_bounds: NDArray[np.float_]
+
+        self.weights: Optional[NDArray[np.float_]] = None
+        self.objective_value: Optional[float] = None
+
+    def compute_quadratic_coeffs(self) -> None:
+        # minimize 0.5 x.T P x + q.T x
+        #   where P = self.quad_obj =   ρ_model.T ρ_model
+        #         q = self.lin_obj  = - ρ_model.T ρ_obs
+        quad_obj_coeffs = np.inner(self.models, self.models)
+        lin_obj_coeffs = -np.inner(self.models, self.target)
+
+        # OSQP requires quadratic objectives in a scipy CSC sparse matrix
+        self.quad_obj = sci.sparse.csc_matrix(quad_obj_coeffs)
+        self.lin_obj = lin_obj_coeffs
+
+    def compute_constraints(self) -> None:
+        # subject to:
+        #   l ≤ A x ≤ u
+        #   where l = self.lower_bounds
+        #         A = self.constraints
+        #         u = self.upper_bounds
+        #
+        #   Each weight x falls in the closed interval [0..1], and the sum of all weights is <= 1.
+        #   This corresponds to (1 * nconformers + 1) constraints, imposed on (nconformers) variables:
+        #     the lower and upper bounds account for (nconformers) constraints,
+        #     the summation constraint accounts for 1.
+        shape = (self.nconformers + 1, self.nconformers)
+        #   We construct A with a sparse matrix.
+        #         constraint idx=[0..N)           constraint idx=N
+        rowidxs = [*range(0, self.nconformers)] + (self.nconformers * [self.nconformers])  # fmt:skip
+        #         0 ≤ x_i ≤ 1                     0 ≤ Σ_i x_i ≤ 1
+        coeffs  = (self.nconformers * [1.0])    + (self.nconformers * [1.0])  # fmt:skip
+        colidxs = [*range(0, self.nconformers)] + [*range(0, self.nconformers)]  # fmt:skip
+        lowers  = (self.nconformers * [0.0])    + [0.0]  # fmt:skip
+        uppers  = (self.nconformers * [1.0])    + [1.0]  # fmt:skip
+
+        self.constraints = sci.sparse.csc_matrix(
+            (coeffs, (rowidxs, colidxs)),
+            shape=shape,
+        )
+        self.lower_bounds = np.array(lowers)
+        self.upper_bounds = np.array(uppers)
+
+    def solve_qp(self) -> None:
+        if TYPE_CHECKING:
+            self.driver = self.osqp
+            assert self.driver is not None
+
+        self.compute_quadratic_coeffs()
+        self.compute_constraints()
+
+        qp = self.driver.OSQP()
+        qp.setup(
+            P=self.quad_obj,
+            q=self.lin_obj,
+            A=self.constraints,
+            l=self.lower_bounds,
+            u=self.upper_bounds,
+            **self.OSQP_SETTINGS,
+        )
+        result = qp.solve()
+
+        self.weights = np.array(result.x).ravel()
+        self.objective_value = (
+            2 * result.info.obj_val
+            + np.inner(self.target, self.target)
+        )  # fmt:skip
+
+
+class MIOSQPSolver(MIQPSolver):
+    driver_pkg_name = "miosqp"
+    driver = lazy_load_module_if_available(driver_pkg_name)
+    if TYPE_CHECKING:
+        import miosqp
+
+    MIOSQP_SETTINGS = {
+        # integer feasibility tolerance
+        "eps_int_feas": 1e-06,
+        # maximum number of iterations
+        "max_iter_bb": 10000,
+        # tree exploration rule
+        #   [0] depth first
+        #   [1] two-phase: depth first until first incumbent and then best bound
+        "tree_explor_rule": 0,
+        # branching rule
+        #   [0] max fractional part
+        "branching_rule": 0,
+        "verbose": False,
+        "print_interval": 1,
+    }
+
+    OSQP_SETTINGS = {
+        "eps_abs": 1e-06,
+        "eps_rel": 1e-06,
+        "eps_prim_inf": 1e-07,
+        "verbose": False,
+    }
+
+    def __init__(self, target: NDArray[np.float_], models: NDArray[np.float_]) -> None:
+        self.target = target
+        self.models = models
+
+        self.nconformers = models.shape[0]
+
+        self.quad_obj: Optional[sci.sparse.csc_matrix] = None
+        self.lin_obj: Optional[NDArray[np.float_]] = None
+        self.constraints: sci.sparse.csc_matrix
+        self.lower_bounds: NDArray[np.float_]
+        self.upper_bounds: NDArray[np.float_]
+        self.binary_vars: NDArray[np.int_]
+        self.lower_ints: NDArray[np.int_]
+        self.upper_ints: NDArray[np.int_]
+
+        self.weights: Optional[NDArray[np.float_]] = None
+        self.objective_value: Optional[float] = None
+
+        # Append raisable Exceptions to SolverError class in module (global) scope
+        global SolverError
+        SolverError += (ValueError,)
+
+    def compute_quadratic_coeffs(self) -> None:
+        """Precompute the quadratic coefficients (P, q).
+
+        These values don't depend on threshold/cardinality.
+        Having these objectives pre-computed saves a few cycles when MIQP
+        is evaluated for multiple values of threshold/cardinality.
+        """
+        # Since this is an MIQP problem, the objective function f(x)
+        # takes the vector x = [w_0 .. w_i, z_0 .. z_i],
+        #   where w_i are the weights,
+        #         z_i are the integer selections.
+
+        # We wish to minimize 0.5 x.T P x + q.T x
+        #   where P = self.quad_obj =   ρ_model.T ρ_model
+        #         q = self.lin_obj  = - ρ_model.T ρ_obs
+        quad_obj_coeffs = np.inner(self.models, self.models)
+        lin_obj_coeffs = -np.inner(self.models, self.target)
+
+        # This is sufficient for the non-integer weights,
+        #   but we must extend our matrix so that
+        #     P is of shape (2*nconfs, 2*nconfs), and
+        #     q is of len (2*nconfs),
+        #   to match the dimensions of vector x.
+        P_shape = (2 * self.nconformers, 2 * self.nconformers)
+
+        # MIOSQP requires quadratic objectives in a scipy CSC sparse matrix
+        # Get an index into the dense _quad_obj_coeffs array,
+        #   and construct P with csc_array((data, (row_ind, col_ind)), [shape=(M, N)]).
+        rowidx, colidx = np.indices(quad_obj_coeffs.shape)
+        self.quad_obj = sci.sparse.csc_matrix(
+            (quad_obj_coeffs.ravel(), (rowidx.ravel(), colidx.ravel())),
+            shape=P_shape,
+        )
+        # Extend q to appropriate shape
+        self.lin_obj = np.append(lin_obj_coeffs, np.zeros((self.nconformers,)))
+
+    def compute_mixed_int_constraints(
+        self,
+        threshold: Optional[float],
+        cardinality: Optional[int],
+    ) -> None:
+        # subject to:
+        #   l ≤ A x ≤ u
+        #   x[i] ∈ Z, for i in i_idx
+        #   il[i] <= x[i] <= iu[i], for i in i_idx
+        #   where l = self.lower_bounds
+        #         A = self.constraints
+        #         u = self.upper_bounds
+        #         i_idx = self.binary_vars: a vector of indices of which variables are integer
+        #         il = self.lower_ints: the lower bounds on the integer variables
+        #         iu = self.upper_ints: the upper bounds on the integer variables.
+
+        # We will construct A with a sparse matrix.
+
+        # Presently, we have no constraints.
+        n_constraints = 0
+        # Our constraints will be applied over the variables [w_0 .. w_i, z_0 .. z_i],
+        #   where w_i are the weights,
+        #         z_i are the integer selections.
+        n_vars = 2 * self.nconformers
+
+        # Since cardinality or threshold will be specified, the problem is a MIQP,
+        # so we need to add binary integer variables z_i.
+        #     z_i ∈ {0, 1}
+        self.binary_vars = np.arange(self.nconformers, 2 * self.nconformers)
+        self.lower_ints = np.array(self.nconformers * [0])
+        self.upper_ints = np.array(self.nconformers * [1])
+
+        # fmt:off
+
+        #   Each weight w_i falls in the closed interval [0..1], and the sum of all weights is <= 1.
+        #   This corresponds to (1 * nconformers + 1) constraints, imposed on (nconformers) variables:
+        #     the lower and upper bounds account for (nconformers) constraints,
+        #     the summation constraint accounts for 1.
+        #        constraint idx=[0..N)           constraint idx=N
+        rowidx = [*range(0, self.nconformers)] + (self.nconformers * [self.nconformers])
+        #        0 ≤ w_i ≤ 1                     0 ≤ Σ_i w_i ≤ 1
+        coeffs = (self.nconformers * [1.0])    + (self.nconformers * [1.0])
+        colidx = [*range(0, self.nconformers)] + [*range(0, self.nconformers)]
+        lowers = (self.nconformers * [0.0])    + [0.0]
+        uppers = (self.nconformers * [1.0])    + [1.0]
+        n_constraints += self.nconformers + 1
+
+        # Introduce an implicit cardinality constraint
+        # Only count weights for which z_i is 1
+        #     0 <= z_i - w_i <= 1
+        #     (∵ z_i ∈ {0,1} and 0 ≤ w_i ≤ 1, this is only true when z_i = 1)
+        # constraint idx=[N+1..2N+1)
+        rowidx += (
+            [*range(n_constraints, n_constraints + self.nconformers)] +
+            [*range(n_constraints, n_constraints + self.nconformers)]
+        )
+        #    0 <= -w_i                           + z_i <= 1
+        coeffs += ((self.nconformers * [-1.0])   + (self.nconformers * [1.0]))
+        colidx += ([*range(0, self.nconformers)] + [*range(self.nconformers, 2 * self.nconformers)])
+        lowers += self.nconformers * [0.0]
+        uppers += self.nconformers * [1.0]
+        n_constraints += self.nconformers
+
+        # Set the threshold constraint if applicable
+        #     tdmin z_i - w_i ≤ 0, i.e. w_i ≥ tdmin
+        if threshold is not None:
+            # constraint idx=[2N+1..3N+1)
+            rowidx += (
+                [*range(n_constraints, n_constraints + self.nconformers)] +
+                [*range(n_constraints, n_constraints + self.nconformers)]
+            )
+            #    0 <= wi                            - t * zi <= 1
+            coeffs += (self.nconformers * [1.0])    + (self.nconformers * [-threshold])
+            colidx += [*range(0, self.nconformers)] + [*range(self.nconformers, 2 * self.nconformers)]
+            lowers += (self.nconformers * [0.0])
+            uppers += (self.nconformers * [1.0])
+            n_constraints += self.nconformers
+
+        # fmt:on
+
+        # Set the cardinality constraint if applicable
+        #     Σ z_i ≤ cardinality
+        if cardinality is not None:
+            # constraint idx=2N+2, if no threshold constraint; or
+            # constraint idx=3N+2, if a threshold constraint was applied
+            rowidx += self.nconformers * [n_constraints]
+            #    0 <= Σ z_i <= cardinality
+            coeffs += self.nconformers * [1]
+            colidx += [*range(self.nconformers, 2 * self.nconformers)]
+            lowers += [0]
+            uppers += [cardinality]
+            n_constraints += 1
+
+        # MIOSQP requires constraints as a sparse matrix
+        self.constraints = sci.sparse.csc_matrix(
+            (coeffs, (rowidx, colidx)),
+            shape=(n_constraints, n_vars),
+        )
+        self.lower_bounds = np.array(lowers)
+        self.upper_bounds = np.array(uppers)
+
+    def solve_miqp(
+        self,
+        threshold: Optional[float] = None,
+        cardinality: Optional[int] = None,
+        exact: bool = False,
+    ) -> None:
+        if TYPE_CHECKING:
+            self.driver = self.miosqp
+            assert self.driver is not None
+
+        if cardinality is threshold is None:
+            raise ValueError("Set either cardinality or threshold.")
+
+        if not (self.quad_obj and self.lin_obj):
+            self.compute_quadratic_coeffs()
+        self.compute_mixed_int_constraints(threshold, cardinality)
+
+        # Construct the MIOSQP solver & solve!
+        miqp = self.driver.MIOSQP()
+        miqp.setup(
+            P=self.quad_obj,
+            q=self.lin_obj,
+            A=self.constraints,
+            l=self.lower_bounds,
+            u=self.upper_bounds,
+            i_idx=self.binary_vars,
+            i_l=self.lower_ints,
+            i_u=self.upper_ints,
+            settings=self.MIOSQP_SETTINGS,
+            qp_settings=self.OSQP_SETTINGS,
+        )
+        result = miqp.solve()
+
+        # Destructure results
+        self.weights = np.array(result.x[0 : self.nconformers])
+        self.objective_value = (
+            2 * result.upper_glob
+            + np.inner(self.target, self.target)
+        )  # fmt:skip
+
+
+###############################
+# Helper methods
+###############################
+
+
+def _available_qp_solvers() -> dict[str, type]:
+    """List all available QP solver classes in this module."""
+    available_solvers = {}
+
+    # Get all classes defined in this module
+    #   use module.__dict__ because it preserves order
+    #     (unlike dir(module) or inspect.getmembers(module))
+    for name, obj in sys.modules[__name__].__dict__.items():
+        if inspect.isclass(obj) and obj.__module__ == __name__:
+            # Check the class implements QPSolver
+            if obj in QPSolver.__subclasses__():
+                # Check the driver module is loadable
+                if obj.driver is not None:
+                    available_solvers[name] = obj
+    return available_solvers
+
+
+def _available_miqp_solvers() -> dict[str, type]:
+    """List all available MIQP solver classes in this module."""
+    available_solvers = {}
+
+    # Get all classes defined in this module
+    #   use module.__dict__ because it preserves order
+    #     (unlike dir(module) or inspect.getmembers(module))
+    for name, obj in sys.modules[__name__].__dict__.items():
+        if inspect.isclass(obj) and obj.__module__ == __name__:
+            # Check the class implements MIQPSolver
+            if obj in MIQPSolver.__subclasses__():
+                # Check the driver module is loadable
+                if obj.driver is not None:
+                    available_solvers[name] = obj
+    return available_solvers
+
+
+available_qp_solvers = _available_qp_solvers()
+available_miqp_solvers = _available_miqp_solvers()
+if not available_qp_solvers:
     msg = (
-        "CPLEX or CVXOPT is not install. Please re-install these packages and qFit to run"
-          )
-    raise RuntimeError(msg)
-else:
-    CPLEX = True
-    SolverError = cplex.exceptions.CplexSolverError
-
-
-# Only these classes should be "public"
-__all__ = ["QPSolver", "MIQPSolver", "SolverError"]
-
-
-# Define the required functions for (MI)QP solver objects
-class _Base_QPSolver(object):
-    """Base class for Quadratic Programming solvers.
-
-    Declares required interface functions for child classes to overwrite."""
-
-    def initialize(self):
-        raise NotImplementedError
-
-    def solve(self):
-        raise NotImplementedError
-
-
-# If we can load the CPLEX solver set,
-# provide class definitions for a QP and MIQP solver.
-if CPLEX:
-    cvxopt.solvers.options["show_progress"] = False
-    cvxopt.solvers.options["abstol"] = 1e-8
-    cvxopt.solvers.options["reltol"] = 1e-7
-    cvxopt.solvers.options["feastol"] = 1e-8
+        "Could not find any QP solver engines.\n"
+        + "Please ensure that at least one of:\n  "
+        + str([solver.driver_pkg_name for solver in QPSolver.__subclasses__()])
+        + "\n"
+        + "is installed."
+    )
+    raise ImportError(msg)
+if not available_miqp_solvers:
+    msg = (
+        "Could not find any MIQP solver engines.\n"
+        + "Please ensure that at least one of:\n  "
+        + str([solver.driver_pkg_name for solver in MIQPSolver.__subclasses__()])
+        + "\n"
+        + "is installed."
+    )
+    raise ImportError(msg)
 
-    class CPLEX_QPSolver(_Base_QPSolver):
-        """Quadratic Programming solver based on CPLEX."""
 
-        def __init__(self, target, models):
-            self._target = target
-            self._models = models
-            self._nconformers = models.shape[0]
-            self.initialized = False
+def get_qp_solver_class(solver_type: str) -> type[QPSolver]:
+    """Return the class of the requested solver type, or raise a KeyError."""
+    return available_qp_solvers[solver_type]
 
-            self._solution = None
-            self.weights = None
 
-        def initialize(self):
-            # Set up the matrices and restraints
-            logger.debug(
-                f"Building cvxopt matrix, size: ({self._nconformers},{self._nconformers})"
-            )
-            self._quad_obj = cvxopt.matrix(np.inner(self._models, self._models), tc="d")
-            self._lin_obj = cvxopt.matrix(-np.inner(self._models, self._target), tc="d")
-
-            # lower-equal constraints.
-            # Each weight falls in the closed interval [0..1] and its sum is <= 1.
-            # There are 2 * nconformers bounds + 1 constraint
-            # Make a sparse matrix to represent this information.
-            self._le_constraints = cvxopt.spmatrix(
-                self._nconformers * [-1.0] + 2 * self._nconformers * [1.0],
-                list(range(2 * self._nconformers))
-                + self._nconformers * [2 * self._nconformers],
-                3 * list(range(self._nconformers)),
-            )
-            self._le_bounds = cvxopt.matrix(
-                self._nconformers * [0.0] + (self._nconformers + 1) * [1.0], tc="d"
-            )
-            self.initialized = True
-
-        def solve(self):
-            if not self.initialized:
-                self.initialize()
-
-            self._solution = cvxopt.solvers.qp(
-                self._quad_obj, self._lin_obj, self._le_constraints, self._le_bounds
-            )
-            self.obj_value = 2 * self._solution["primal objective"] + np.inner(
-                self._target, self._target
-            )
-            self.weights = np.asarray(self._solution["x"]).ravel()
-
-    class CPLEX_MIQPSolver(_Base_QPSolver):
-        """Mixed-Integer Quadratic Programming solver based on CPLEX."""
-
-        def __init__(self, target, models, threads=1):
-            self._target = target
-            self._models = models
-            self._nconformers = models.shape[0]
-            self.initialized = False
-            self.threads = threads
-
-        def initialize(self):
-            self._quad_obj = np.inner(self._models, self._models)
-            self._lin_obj = -np.inner(self._models, self._target)
-
-            self.initialized = True
-
-        def solve(self, cardinality=None, exact=False, threshold=None):
-            if not self.initialized:
-                self.initialize()
-
-            miqp = cplex.Cplex()
-            miqp.set_results_stream(None)
-            miqp.set_log_stream(None)
-            miqp.set_warning_stream(None)
-            miqp.set_error_stream(None)
-
-            # Set number of threads to use
-            if self.threads is not None:
-                miqp.parameters.threads.set(self.threads)
-
-            # Setup QP part of the MIQP
-            variable_names = [f"w{n}" for n in range(self._nconformers)]
-            upper_bounds = self._nconformers * [1.0]
-
-            miqp.variables.add(names=variable_names, ub=upper_bounds)
-            for i in range(self._nconformers):
-                for j in range(i, self._nconformers):
-                    miqp.objective.set_quadratic_coefficients(
-                        i, j, self._quad_obj[i, j]
-                    )
-                miqp.objective.set_linear(i, self._lin_obj[i])
-            ind = range(self._nconformers)
-            val = [1] * self._nconformers
-            lin_expr = [cplex.SparsePair(ind=ind, val=val)]
-            miqp.linear_constraints.add(
-                lin_expr=lin_expr,
-                rhs=[1],
-                senses=["L"],
-            )  # Sum of weights is <= 1
-
-            # If cardinality or threshold is specified the problem is a MIQP, else its
-            # a regular QP.
-            if cardinality not in (None, 0) or threshold not in (None, 0):
-                integer_names = [f"z{n}" for n in range(self._nconformers)]
-                variable_types = self._nconformers * miqp.variables.type.binary
-                miqp.variables.add(names=integer_names, types=variable_types)
-
-                # Only count weights for which zi is 1
-                for n in range(self._nconformers):
-                    w = f"w{n}"
-                    z = f"z{n}"
-                    miqp.linear_constraints.add(
-                        lin_expr=[cplex.SparsePair(ind=[w, z], val=[1, -1])],
-                        rhs=[0],
-                        senses="L",
-                    )
-                    # Set the threshold constraint
-                    if threshold not in (None, 0):
-                        miqp.linear_constraints.add(
-                            lin_expr=[
-                                cplex.SparsePair(ind=[z, w], val=[threshold, -1])
-                            ],
-                            rhs=[0],
-                            senses=["L"],
-                        )
-                # Set the cardinality constraint
-                if cardinality not in (None, 0):
-                    senses = "L"
-                    if exact:
-                        senses = "E"
-                        cardinality = min(cardinality, self._nconformers)
-                    lin_expr = [
-                        [
-                            list(range(self._nconformers, 2 * self._nconformers)),
-                            self._nconformers * [1],
-                        ]
-                    ]
-                    miqp.linear_constraints.add(
-                        lin_expr=lin_expr,
-                        rhs=[cardinality],
-                        senses=senses,
-                    )
-            miqp.solve()
-
-            self.obj_value = 2 * miqp.solution.get_objective_value() + np.inner(
-                self._target, self._target
-            )
-            self.weights = np.asarray(miqp.solution.get_values()[: self._nconformers])
-            miqp.end()
-            q = self._lin_obj.reshape(-1, 1)
-            P = self._quad_obj
-            w = self.weights.reshape(-1, 1)
-
-            # print("CPLEX MIQP")
-            # print('P:', P)
-            # print('q:', q)
-            # print('w:', w)
-
-            # obj = 0.5 * w.T @ P @ w + q.T @ w
-            # print("calculated myself OBJ:", obj)
-            # print('from solver OBJ:', self.obj_value)
-            # print("TOTAL:", self.obj_value * 2 + np.inner(self._target, self._target))
-
-
-# Create a "pseudo-class" to abstract the choice of solver.
-class QPSolver(object):
-    def __new__(cls, *args, **kwargs):
-        """Return a constructor for the appropriate solver."""
-
-        # Default to using cplex if not specified.
-        use_cplex = kwargs.pop("use_cplex", True)
-
-        # Walk through solver decision tree.
-        if use_cplex:
-            if CPLEX:
-                return CPLEX_QPSolver(*args, **kwargs)
-            else:
-                raise ImportError(
-                    "qFit could not load modules for Quadratic Programming solver.\n"
-                    "Please install cvxopt & CPLEX."
-                )
-        else:
-            if CPLEX:
-                print(
-                    "WARNING: A different solver was requested, but only CPLEX solvers found.\n"
-                    "         Using CPLEX solver as fallback."
-                )
-                return CPLEX_QPSolver(*args, **kwargs)
-            else:
-                raise ImportError(
-                    "qFit could not load modules for Quadratic Programming solver.\n"
-                    "Please install cvxopt & CPLEX."
-                )
-
-    def __init__(self, *args, **kwargs):
-        # This pseudo-class does not generate class instances.
-        # As there are no instances of this class,
-        #   it does not need an instance initialiser.
-        raise NotImplementedError
-
-
-# Create a "pseudo-class" to abstract the choice of solver.
-class MIQPSolver(object):
-    def __new__(cls, *args, **kwargs):
-        """Construct an initialised instance of the appropriate solver."""
-
-        # Default to using cplex if not specified.
-        use_cplex = kwargs.pop("use_cplex", True)
-
-        # Walk through solver decision tree.
-        if use_cplex:
-            if CPLEX:
-                return CPLEX_MIQPSolver(*args, **kwargs)
-            else:
-                raise ImportError(
-                    "qFit could not load modules for Quadratic Programming solver.\n"
-                    "Please install cvxopt & CPLEX."
-                )
-        else:
-            if CPLEX:
-                print(
-                    "WARNING: A different solver was requested, but only CPLEX solvers found.\n"
-                    "         Using CPLEX solver as fallback."
-                )
-                return CPLEX_MIQPSolver(*args, **kwargs)
-            else:
-                raise ImportError(
-                    "qFit could not load modules for Quadratic Programming solver.\n"
-                    "Please install cvxopt & CPLEX."
-                )
-
-    def __init__(self, *args, **kwargs):
-        # This pseudo-class does not generate class instances.
-        # As there are no instances of this class,
-        #   it does not need an instance initialiser.
-        raise NotImplementedError
+def get_miqp_solver_class(solver_type: str) -> type[MIQPSolver]:
+    """Return the class of the requested solver type, or raise a KeyError."""
+    return available_miqp_solvers[solver_type]
diff --git a/src/qfit/utils/optional_lazy_import.py b/src/qfit/utils/optional_lazy_import.py
new file mode 100644
index 0000000..a0c0bd5
--- /dev/null
+++ b/src/qfit/utils/optional_lazy_import.py
@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+
+import importlib.util
+import sys
+from types import ModuleType
+from typing import Optional
+
+
+def lazy_load_module_if_available(name: str) -> Optional[ModuleType]:
+    """Import a module if available, otherwise return None."""
+    # If it's already loaded, no need to re-import
+    if name in sys.modules:
+        return sys.modules[name]
+
+    # Can we find the module?
+    spec = importlib.util.find_spec(name)
+    if spec is None:
+        return None
+    if not hasattr(spec, "loader") or spec.loader is None:
+        return None
+
+    # Lazy load, and put in sys.modules
+    module = importlib.util.module_from_spec(spec)
+    loader = importlib.util.LazyLoader(spec.loader)
+    loader.exec_module(module)
+    sys.modules[name] = module
+    return module


Test patch:
diff --git a/tests/test_qfit_ligand.py b/tests/test_qfit_ligand.py
index 83689a9..c54c0cb 100644
--- a/tests/test_qfit_ligand.py
+++ b/tests/test_qfit_ligand.py
@@ -10,6 +10,10 @@ from qfit.logtools import (
     setup_logging,
     log_run_info,
 )
+from qfit.solvers import (
+    available_qp_solvers,
+    available_miqp_solvers,
+)
 
 
 logger = logging.getLogger(__name__)
@@ -40,6 +44,8 @@ class TestQFitLigand:
         options = QFitOptions()
         options.apply_command_args(args)
         options.debug = True  # For debugging in tests
+        options.qp_solver = next(iter(available_qp_solvers.keys()))
+        options.miqp_solver = next(iter(available_miqp_solvers.keys()))
 
         # Setup logger
         setup_logging(options=options, filename="qfit_ligand.log")
diff --git a/tests/test_solvers.py b/tests/test_solvers.py
index edaefa5..9829f81 100644
--- a/tests/test_solvers.py
+++ b/tests/test_solvers.py
@@ -1,29 +1,158 @@
+import inspect
+
+import numpy as np
 import pytest
+from numpy.typing import NDArray
+
+import qfit.solvers
+from qfit.solvers import (
+    available_miqp_solvers,
+    available_qp_solvers,
+    get_miqp_solver_class,
+    get_qp_solver_class,
+)
+
+
+def test_missing_solvers() -> None:
+    with pytest.raises(KeyError):
+        get_qp_solver_class("NotASolver")
+    with pytest.raises(KeyError):
+        get_miqp_solver_class("NotASolver")
+
+
+def test_get_qp_solver() -> None:
+    qp_solver_class = get_qp_solver_class(next(iter(available_qp_solvers.keys())))
+    assert inspect.isclass(qp_solver_class)
+    assert issubclass(qp_solver_class, qfit.solvers.QPSolver)
+
+
+def test_get_miqp_solver() -> None:
+    miqp_solver_class = get_qp_solver_class(next(iter(available_miqp_solvers.keys())))
+    assert inspect.isclass(miqp_solver_class)
+    assert issubclass(miqp_solver_class, qfit.solvers.MIQPSolver)
+
+
+@pytest.mark.parametrize("solver_class", available_qp_solvers.values())
+class TestQPSolver:
+    target = np.array([2.0, 3.0, 7.0])
+    models = np.array([[6.0, 0.0, 0.0], [0.0, 9.0, 0.0], [0.0, 0.0, 21.0]])
+
+    def test_qp_solver(self, solver_class: type[qfit.solvers.QPSolver]) -> None:
+        solver = solver_class(self.target, self.models)
+        solver.solve_qp()
+
+        assert np.allclose(solver.weights, [1 / 3, 1 / 3, 1 / 3], atol=1e-3)
+        assert np.isclose(solver.objective_value, 0.0, atol=1e-6)
+
+
+@pytest.mark.parametrize("solver_class", available_miqp_solvers.values())
+class TestMIQPSolver:
+    target = np.array([2.0, 3.0, 7.0])
+    models = np.array([[6.0, 0.0, 0.0], [0.0, 9.0, 0.0], [0.0, 0.0, 21.0]])
+
+    def expected_objective(self, expected_weights: NDArray[np.float_]) -> float:
+        return np.sum(np.square(np.inner(self.models, expected_weights) - self.target))
+
+    def test_miqp_solver_with_threshold(
+        self, solver_class: type[qfit.solvers.MIQPSolver]
+    ) -> None:
+        solver = solver_class(self.target, self.models)
+        solver.solve_miqp(threshold=0.4)
+
+        expected_weights = np.array([0.0, 0.4, 0.4])
+        expected_objective = self.expected_objective(expected_weights)
+
+        assert np.allclose(solver.weights, expected_weights, atol=1e-3)
+        assert np.isclose(solver.objective_value, expected_objective, atol=1e-6)
+
+    def test_miqp_solver_with_cardinality_3(
+        self, solver_class: type[qfit.solvers.MIQPSolver]
+    ) -> None:
+        solver = solver_class(self.target, self.models)
+        solver.solve_miqp(cardinality=3)
+
+        expected_weights = np.array([1 / 3, 1 / 3, 1 / 3])
+        expected_objective = 0.0
+        assert np.allclose(solver.weights, expected_weights, atol=1e-3)
+        assert np.isclose(solver.objective_value, expected_objective, atol=1e-6)
+
+    def test_miqp_solver_with_cardinality_2(
+        self, solver_class: type[qfit.solvers.MIQPSolver]
+    ) -> None:
+        solver = solver_class(self.target, self.models)
+        solver.solve_miqp(cardinality=2)
+
+        expected_weights = np.array([0.0, 1 / 3, 1 / 3])
+        expected_objective = self.expected_objective(expected_weights)
+        assert np.allclose(solver.weights, expected_weights, atol=1e-3)
+        assert np.isclose(solver.objective_value, expected_objective, atol=1e-6)
+
+    def test_miqp_solver_with_cardinality_1(
+        self, solver_class: type[qfit.solvers.MIQPSolver]
+    ) -> None:
+        solver = solver_class(self.target, self.models)
+        solver.solve_miqp(cardinality=1)
+
+        expected_weights = np.array([0.0, 0.0, 1 / 3])
+        expected_objective = self.expected_objective(expected_weights)
+        assert np.allclose(solver.weights, expected_weights, atol=1e-3)
+        assert np.isclose(solver.objective_value, expected_objective, atol=1e-6)
+
+    def test_miqp_solver_with_threshold_and_cardinality_1(
+        self, solver_class: type[qfit.solvers.MIQPSolver]
+    ) -> None:
+        solver = solver_class(self.target, self.models)
+        solver.solve_miqp(threshold=0.4, cardinality=1)
+
+        expected_weights = np.array([0.0, 0.0, 0.4])
+        expected_objective = self.expected_objective(expected_weights)
+        assert np.allclose(solver.weights, expected_weights, atol=1e-3)
+        assert np.isclose(solver.objective_value, expected_objective, atol=1e-6)
+
+
+@pytest.mark.parametrize("solver_class", available_miqp_solvers.values())
+class TestMIQPSolverReuse:
+    target = np.array([2.0, 3.0, 7.0])
+    models = np.array([[6.0, 0.0, 0.0], [0.0, 9.0, 0.0], [0.0, 0.0, 21.0]])
+
+    @pytest.fixture
+    def solver(
+        self, solver_class: type[qfit.solvers.MIQPSolver]
+    ) -> qfit.solvers.MIQPSolver:
+        """Instantiate a solver within the scope of this class, to be re-used."""
+        return solver_class(self.target, self.models)
+
+    def expected_objective(self, expected_weights: NDArray[np.float_]) -> float:
+        return np.sum(np.square(np.inner(self.models, expected_weights) - self.target))
 
-from qfit.solvers import *
+    def test_miqp_solver_with_threshold_and_cardinality_1(
+        self, solver_class: type[qfit.solvers.MIQPSolver]
+    ) -> None:
+        solver = solver_class(self.target, self.models)
+        solver.solve_miqp(threshold=0.4, cardinality=1)
 
+        expected_weights = np.array([0.0, 0.0, 0.4])
+        expected_objective = self.expected_objective(expected_weights)
+        assert np.allclose(solver.weights, expected_weights, atol=1e-3)
+        assert np.isclose(solver.objective_value, expected_objective, atol=1e-6)
 
-def setup_module(module):
-    import qfit.solvers
+    def test_miqp_solver_with_cardinality_1(
+        self, solver: type[qfit.solvers.MIQPSolver]
+    ) -> None:
+        solver.solve_miqp(cardinality=1)
 
-    # Manually set CPLEX to False
-    # in this module ONLY
-    qfit.solvers.CPLEX = False
+        expected_weights = np.array([0.0, 0.0, 1 / 3])
+        expected_objective = self.expected_objective(expected_weights)
+        assert np.allclose(solver.weights, expected_weights, atol=1e-3)
+        assert np.isclose(solver.objective_value, expected_objective, atol=1e-6)
 
+    def test_miqp_solver_with_threshold(
+        self, solver: type[qfit.solvers.MIQPSolver]
+    ) -> None:
+        solver.solve_miqp(threshold=0.4)
 
-class TestForcedNoSolvers:
-    # QPSolver.__new__ should raise ImportError
-    # when CPLEX flags are false
-    def test_QPSolvers_fail(self):
-        with pytest.raises(ImportError):
-            QPSolver(None, None, use_cplex=False)
-        with pytest.raises(ImportError):
-            QPSolver(None, None, use_cplex=True)
+        expected_weights = np.array([0.0, 0.4, 0.4])
+        expected_objective = self.expected_objective(expected_weights)
 
-    # MIQPSolver.__new__ should raise ImportError
-    # when CPLEX flags are false
-    def test_MIQPSolvers_fail(self):
-        with pytest.raises(ImportError):
-            MIQPSolver(None, None, use_cplex=False)
-        with pytest.raises(ImportError):
-            MIQPSolver(None, None, use_cplex=True)
+        assert np.allclose(solver.weights, expected_weights, atol=1e-3)
+        assert np.isclose(solver.objective_value, expected_objective, atol=1e-6)