miom.py

class BaseModel(ABC):
    """Base class for building LP/MIP optimization models
    using a high-level API for metabolic problems.

    It implements the chainable methods to set-up a LP/MIP 
    metabolic network problem. Two implementations are available: 
    PicosModel and PythonMipModel. The PicosModel uses the 
    [PICOS library](https://picos-api.gitlab.io/picos/) 
    as a backend to interact with different solvers. The PythonMipModel uses 
    the [Python-MIP](https://www.python-mip.com/) library to solve the
    model using the CBC or GUROBI solvers.

    !!! note
        Do not try to instantiate this class directly. Use the [load()][miom.miom.load] 
        function instead. The method automatically selects the right implementation 
        depending on the solver.

    """
    def __init__(self, previous_step_model=None, miom_network=None, solver_name=None):
        self.previous_step_model = previous_step_model
        self.network = miom_network
        self.variables = None
        self.objective = None
        self._last_start_time = None
        self.last_solver_time = None
        if previous_step_model is not None:
            self._options = previous_step_model._options
            if miom_network is None:
                self.network = previous_step_model.network
            if solver_name is not None:
                self._options["solver"] = solver_name
        else:
            # Default recommended options
            self._options = {
                'int_tol': 1e-8,
                'feas_tol': 1e-8,
                'opt_tol': 1e-5,
                'verbosity': 0,
                'solver': solver_name
            }
        self.problem = self.initialize_problem()
        self.setup(**self._options)

    @abstractmethod
    def initialize_problem(self):
        pass

    @abstractmethod
    def get_solver_status(self):
        pass

    @property
    def status(self):
        return self.get_solver_status()


    @_composable
    def setup(self, **kwargs):
        """Provide the options for the solver.

        Attributes:
            int_tol (float): Integrality tolerance for integer variables.
                Defaults to 1e-8.
            feas_tol (float): Feasibility tolerance. Defaults to 1e-8.
            opt_tol (float): Relative MIP gap tolerance for MIP problems.
                Defaults to 1e-5.
            verbosity (int): Values above 0 force the backends to be verbose.
                Use a value of 1 to show useful information about the search.
                Defaults to 0.

        Returns:
            BaseModel: the configured instance of the BaseModel
        """
        if self.problem is None:
            warnings.warn("Problem cannot be configured since it was not initialized")
            return False
        options = dict()
        int_tol = kwargs["int_tol"] if "int_tol" in kwargs else None
        feas_tol = kwargs["feas_tol"] if "feas_tol" in kwargs else None
        opt_tol = kwargs["opt_tol"] if "opt_tol" in kwargs else None
        verbosity = kwargs["verbosity"] if "verbosity" in kwargs else None
        solver = kwargs["solver"] if "solver" in kwargs else None
        if int_tol is not None:
            options["int_tol"] = int_tol
        if feas_tol is not None:
            options["feas_tol"] = feas_tol
        if opt_tol is not None:
            options["opt_tol"] = opt_tol
        if verbosity is not None:
            options["verbosity"] = verbosity
        if solver is not None:
            options["solver"] = solver
        self._options.update(options)
        # Calculate a min eps value to avoid numerical issues
        self._options["_min_eps"] = np.sqrt(2*self._options["int_tol"])
        return self._options

    @_composable
    def steady_state(self):
        """Add the required constraints for finding steady-state fluxes

        The method adds the $S * V = 0$ set of constraints, where $S$
        is the stoichiometric matrix and $V$ the flux variables.

        Returns:
            BaseModel: instance of BaseModel with the modifications applied.
        """
        pass

    @_composable
    def keep(self, reactions):
        """Force the inclusion of a list of reactions in the solution.

        Reactions have to be associated with positive weights in order to
        keep them in the final solution. Note that once keep() is called,
        the weights associated to the reactions selected to be kept will
        not be taken into account, as they will be forced to be kept in
        the solution.

        Args:
            reactions (list): List of reaction names, a binary vector
                indicating the reactions to keep, or a list of indexes
                with the reactions to keep.

        Returns:
            BaseModel: instance of BaseModel with the modifications applied.
        """
        if self.variables.indicators is None:
            raise ValueError("""No indicator variables for reactions, 
                             transform it to a subset selection problem calling
                             subset_selection first, providing positive weights for the
                             reactions you want to keep.""")
        if reactions is None or len(reactions) == 0:
            return False
        if isinstance(reactions, str):
            reactions = [reactions]
        # Check if it's a binary vector
        if len(reactions) == self.network.num_reactions and max(reactions)==1:
            reactions = np.where(reactions==1)[0]
        # Check if the reactions have an indicator variable
        reactions = set(self.network.find_reaction(rxn)[0] for rxn in reactions)
        available = set(rxn.index for rxn in self.variables.assigned_reactions if rxn.cost > 0)
        diff = reactions - available
        if len(diff) != 0:
            raise ValueError(f"""Only reactions associated with positive weights
                             can be forced to be selected. The following reaction 
                             indexes have no indicator variables or are associated with 
                             negative weights: {diff}.""")
        valid_rxn_idx = reactions & available
        # Get the indexes of the indicator variables
        idxs = [i for i, r in enumerate(self.variables.assigned_reactions)
                if r.index in valid_rxn_idx]
        return dict(idxs=idxs, valid_rxn_idx=valid_rxn_idx)

    @_composable
    def subset_selection(self, rxn_weights, direction='max', eps=1e-2, strengthen=True):
        """Transform the current model into a subset selection problem.

        The subset selection problem consists of selecting the positive weighted
        reactions and remove the negative weighted reactions, subject to steady
        state constraints and/or additional constraints on fluxes, and maximizing
        the weighted sum of the (absolute) weights for the successfully selected reactions
        (with positive weights) and the successfully removed reactions (with negative
        weights). Selected reactions are forced to have an absolute flux value greater 
        or equal to the threshold `eps` (1e-2 by default). Removed reactions should have a
        flux equal to 0.

        Each reaction is associated with a weight (positive or negative) provided
        in the parameter `rxn_weights`, and the objective is to select the reactions 
        that optimizes the following expression:

        $$
        f(x) = \sum_i^n |w_i| * x_i
        $$

        where $x_i$ are the indicator variables for the reactions $i$ and $w_i$ are
        the weights for the reactions associated to the indicator variable. Indicator
        variables are automatically created for each reaction associated to a non-zero
        weight. Two (mutually exclusive) indicator variables are used for positive weighted 
        reactions that are reversible to indicate whether there is positive or negative flux 
        through the reaction. A single indicator variable is created for positive weighted
        non-reversible reactions, to indicate if the reaction is selected (has a non-zero 
        flux greater or equal to `eps`) in which case the indicator variable is 1, 
        or 0 otherwise. 

        A single binary indicator variable is also created for negative weighted reactions, 
        indicating whether the reaction was not selected (i.e, has 0 flux, in which case the
        indicator variable is 1) or not (in which case the indicator variable is 0).

        Args:
            rxn_weights (list): List of weights for each reaction. If a single value
                is provided, it is assumed to be the weight for all reactions.
            eps (float, optional): Min absolute flux value for weighted reactions
                to consider them active or inactive. Defaults to 1e-2.

        Returns:
            BaseModel: instance of BaseModel with the modifications applied.
        """
        # Calculate min valid EPS based on integrality tolerance
        min_eps = self._options["_min_eps"]
        if eps < min_eps:
            warnings.warn(f"Minimum epsilon value below min. allowed value, changed to {min_eps}.")
        eps = max(eps, min_eps)
        if not isinstance(rxn_weights, Iterable):
            rxn_weights = [rxn_weights] * self.network.num_reactions
        rxnw = _weighted_rxns(self.network.R, rxn_weights)
        if self.variables.indicators is None:
            self.variables._assigned_reactions = rxnw
            return dict(eps=eps, direction=direction, strengthen=strengthen)
        else:
            raise ValueError("The current model is already a subset selection problem.")

    @_composable
    def set_flux_bounds(self, rxn_id, min_flux=None, max_flux=None):
        """Change the flux bounds of a reaction.

        Args:
            rxn_id (str/int): name or id of the reaction to change
            min_flux (float, optional): Min flux value. Defaults to None.
            max_flux (float, optional): Max flux value. Defaults to None.

        Returns:
            BaseModel: instance of BaseModel with the modifications applied.
        """
        i, _ = self.network.find_reaction(rxn_id)
        return i

    @_composable
    def add_constraints(self, constraints):
        """Add a list of constraint to the model

        Args:
            constraints (list): List of expressions with the
                constraints.

        Returns:
            BaseModel: instance of BaseModel with the modifications applied.
        """
        for c in constraints:
            self.add_constraint(c)
        return len(constraints) > 0

    @_composable
    def exclude(self, indicator_values=None):
        """Exclude a solution from the set of feasible solutions.

        If the problem is a subset selection problem, it adds a new constraint
        to exclude the given values (or the current values of the indicator variables)
        from the set of feasible solutions. This is useful to force the solver to find
        alternative solutions in a manual fashion. https://doi.org/10.1371/journal.pcbi.1008730


        Args:
            indicator_values (list, optional): List of values for each indicator variable. Defaults to None.

        Returns:
            BaseModel: instance of BaseModel with the modifications applied.
        """
        if self.variables.indicators is None:
            raise ValueError("""The optimization model does not contain indicator variables.
            Make sure that the problem is a subset selection problem.""")
        if indicator_values is None:
            indicator_values = np.array(self.variables.indicator_values)
        return dict(indicator_values=indicator_values)

    @_composable
    def add_constraint(self, constraint):
        """Add a specific constraint to the model.
        The constraint should use existing variables already included in the model.

        Args:
            constraint: affine expression using model's variables.
        """
        pass

    @_composable
    def set_objective(self, cost_vector, variables, direction='max'):
        """Set the optmization objective of the model.

        Args:
            cost_vector (Iterable): List with the cost for each variable
            variables (Iterable): Variables used for the objective function
            direction (str, optional): Optimization direction (min or max). Defaults to 'max'.
        """
        if self.objective is not None:
            warnings.warn("Previous objective changed")
        self.objective = (cost_vector, variables, direction)

    def set_rxn_objective(self, rxn, direction='max'):
        """Set a flux objective

        Maximize or minimize the flux through a given reaction.

        Args:
            rxn (str): Name of the reaction to use for the optimization
            direction (str, optional): Minimization or maximization of 
                the flux ('min' or 'max'). Defaults to 'max'.

        Returns:
            BaseModel: instance of BaseModel with the modifications applied.
        """
        i = self.network.get_reaction_id(rxn)
        cost = np.zeros(self.network.num_reactions)
        cost[i] = 1
        self.set_objective(cost, self.variables.fluxvars, direction=direction)
        return self


    def set_fluxes_for(self, reactions, tolerance=1e-6):
        warnings.warn("This method was renamed to fix_fluxes. It will dissappear in the v0.9.0", DeprecationWarning)
        return self.fix_fluxes(reactions, tolerance)


    def fix_fluxes(self, reactions, tolerance=1e-6):
        """Force the flux of certain reactions to match current values.

        After calling `.solve()` for a flux optimization problem (e.g. FBA), this
        method adds a new constraint to force the flux of the given reactions to
        match the current flux values found after the optimization.

        This is interesting for example to implement methods like sparse-FBA, where
        the optimization is no longer the flux but the number of active reactions,
        and a new constraint is required to preserve optimality of fluxes.

        Args:
            reactions (list): reaction or list of reactions
            tolerance (float, optional): Tolerance for the flux values
                (a solution is valid if the flux is within optimal value +/- tol. 
                Defaults to 1e-6.

        Returns:
            BaseModel: instance of BaseModel with the modifications applied.
        """
        if isinstance(reactions, str):
            reactions = [reactions]

        for rid in reactions:
            if isinstance(rid, np.ndarray):
                # TODO: Test this path
                rid = rid['id']
            idx, rxn = self.network.find_reaction(rid)
            lb = max(rxn['lb'], self.variables.flux_values[idx] - tolerance)
            ub = min(rxn['ub'], self.variables.flux_values[idx] + tolerance)
            self.add_constraint(self.variables.fluxvars[idx] >= lb)
            self.add_constraint(self.variables.fluxvars[idx] <= ub)
        return self

    @_composable
    def reset(self):
        """Resets the original problem (removes all modifications)

        Returns:
            BaseModel: instance of BaseModel with the modifications applied.
        """
        if self.problem is None:
            warnings.warn("Problem is not initialized, nothing to reset")
            return False
        else:
            return True

    def obtain_subnetwork(
            self,
            extraction_mode=ExtractionMode.ABSOLUTE_FLUX_VALUE,
            comparator=Comparator.GREATER_OR_EQUAL,
            value=1e-8
    ):
        """Same as [select_subnetwork][miom.miom.BaseModel] but returns the network instead.

        See [select_subnetwork][miom.miom.BaseModel] for a detailed description of the method.

        Returns:
            MiomNetwork: A MiomNetwork with the selected subnetwork.
        """
        # If indicators are present and assigned,
        # take the subset of the network for which
        # the indicators of positive weighted reactions
        # are equal to 1
        if extraction_mode == ExtractionMode.ABSOLUTE_FLUX_VALUE:
            variables = self.variables.flux_values
            if variables is None:
                raise ValueError("The model does not contain flux variables. "
                                 "You need to call first steady_state() to add "
                                 "the flux variables")
        elif extraction_mode == ExtractionMode.INDICATOR_VALUE:
            variables = self.variables.indicator_values
            if variables is None:
                raise ValueError("The model does not contain indicator variables. "
                                 "You need to transform it to a subset selection problem "
                                 "by invoking subset_selection() first.")
        elif extraction_mode == ExtractionMode.REACTION_ACTIVITY:
            variables = self.variables.reaction_activity
            if variables is None:
                raise ValueError("The model does not contain reaction activity variables. "
                                 "You need to transform it to a subset selection problem "
                                 "by invoking subset_selection() first.")
        else:
            raise ValueError("Invalid extraction mode")

        if comparator == Comparator.GREATER_OR_EQUAL:
            selected = np.where(np.abs(variables) >= value)[0]
        elif comparator == Comparator.LESS_OR_EQUAL:
            selected = np.where(np.abs(variables) <= value)[0]
        else:
            raise ValueError("Invalid comparison")
        # TODO: Assigned reactions work only for indicator variables
        if extraction_mode == ExtractionMode.INDICATOR_VALUE:
            rxns = [self.variables.assigned_reactions[x] for x in selected]
            selected_idx = [rxn.index for rxn in rxns]
        elif extraction_mode == ExtractionMode.ABSOLUTE_FLUX_VALUE:
            selected_idx = selected
        else:
            raise NotImplementedError("Only indicator variables and absolute flux values are supported")
        return self.network.subnet(selected_idx)

    @_composable
    def select_subnetwork(
            self,
            mode=ExtractionMode.ABSOLUTE_FLUX_VALUE,
            comparator=Comparator.GREATER_OR_EQUAL,
            value=1e-8
    ):
        """Select a subnetwork and create a new BaseModel to operate on it.

        The new instance of the BaseModel is a new problem instance with no constraints.
        If the idea is to perform FBA simulations on this new subnetwork, remember to add
        the new constraints, especially the `steady_state`.

        Args:
            mode (ExtractionMode, optional): Method used to extract the subnetwork 
                (based on flux values or using the indicator values). 
                Defaults to ExtractionMode.ABSOLUTE_FLUX_VALUE.
            comparator (Comparator, optional): Comparator for the selected mode. 
                Defaults to Comparator.GREATER_OR_EQUAL.
            value (float, optional): Value threshold for the mode and comparator selected. 
                Defaults to 1e-8.

        Returns:
            [type]: [description]
        """
        return self.obtain_subnetwork(extraction_mode=mode,
                                      comparator=comparator,
                                      value=value)

    def get_values(self):
        """Get the values for the variables

        Returns:
            tuple: (V, X) where V are the flux values and X are the indicator values
                (if the problem is a MIP problem, for example if `subset_selection` was
                called)
        """
        return self.variables.values()


    def get_fluxes(self, reactions=None):
        """Get the flux values.

        Args:
            reactions (list, optional): Reaction or subset of reactions
                For which to obtain the flux values. Defaults to None.

        Returns:
            list: List with the flux values for all or the selected reactions.
        """
        if isinstance(reactions, str):
            return self.variables.flux_values[self.network.get_reaction_id(reactions)]
        if isinstance(reactions, Iterable):
            return {
                r['id']: self.variables.flux_values[self.network.get_reaction_id(r['id'])]
                for r in reactions
            }
        if reactions is None:
            return self.variables.flux_values
        else:
            raise ValueError("reactions should be an iterable of strings or a single string")

    @_composable
    def solve(self, verbosity=None, max_seconds=None):
        """Solve the current model and assign the values to the variables of the model.

        Args:
            verbosity (int, optional): Level of verbosity for the solver. 
                Values above 0 will force the backend to show output information of the search. Defaults to None.
            max_seconds (int, optional): Max time in seconds for the search. Defaults to None.
        """
        self._last_start_time = perf_counter()

    @_composable
    def copy(self):
        pass

class Comparator(str, Enum):
    """Comparator enum to use with [select_subnetwork()][miom.miom.BaseModel]

    Attributes:
        GREATER_OR_EQUAL (str): Select those variables 
            with a value greater or equal than the one provided.
        LESS_OR_EQUAL (str): Select those variables 
            with a value less or equal than the one provided.
    """
    GREATER_OR_EQUAL = "geq",
    LESS_OR_EQUAL = "leq"

class ExtractionMode(str, Enum):
    """Method to select the subnetwork to be extracted.

    This mode works with the method [select_subnetwork][miom.miom.BaseModel] 
    only after a subset selection problem was succesfully solved.
    If the model is configured as a subset selection problem, the user can
    extract a subnetwork after the method solve() was called. The selection
    of the subnetwork can be done using the indicator variables or the flux
    variables.

    For more information, please read the documentation of the method
    [subset_selection][miom.miom.BaseModel]

    Attributes:
        ABSOLUTE_FLUX_VALUE (str): Use a decision criterion based on the
            value of the fluxes. For example, by selecting all the reactions
            that have an absolute flux value above a certain threshold.

        INDICATOR_VALUE (str): Use the binary indicator variables to select
            the subnetwork. Binary indicator variables are created for a
            subset selection problem (after calling [subset_selection][miom.miom.BaseModel]).
            Two indicators are created for each positive weighted and reversible reaction,
            to indicate if there is a non-zero positive flux or negative flux respectively.
            A single binary indicator variable is created for each negative weighted reaction.
            In this case, an indicator value of 1 indicates that the reaction was succesfully
            removed, and 0 otherwise. You can use the value of the indicator variables to
            select a subnetwork after solving. For example, if all the reactions have a
            negative weight (since the goal is, for example, to minimize the number of 
            reactions subject to some other constraints) and you want to select only the
            reactions that were not removed after solving, you can use the indicator
            variables with a value of 0.

            Usage example:
            ```python
                m.
                steady_state().
                add_constraints(...).
                # Convert to a subset selection, where each reaction
                # has a weight of -1.
                subset_selection(-1).
                # Solve the optimization problem
                .solve()
                # Get the subnetwork selecting the reactions
                # with an indicator value below 0.5. Since variables
                # are binary, it basically selects all the reactions
                # associated with a binary indicator value of 0.
                # This corresponds to the reactions that were not
                # removed after solving the problem (since all the
                # reactions have negative weights, and their binary
                # indicator variables are 1 if they were successfully
                # removed, and 0 if they could not be removed).
                .select_subnetwork(
                    mode=ExtractionMode.INDICATOR_VALUE,
                    comparator=Comparator.LESS_OR_EQUAL,
                    value=0.5
                ).network
            ```
    """
    ABSOLUTE_FLUX_VALUE = "flux_value",
    INDICATOR_VALUE = "indicator_value",
    REACTION_ACTIVITY = "reaction_activity"

class Solvers(str, Enum):
    """Solvers supported by the miom module.

    Please refer to https://picos-api.gitlab.io/picos/introduction.html to see
    the list of supported solvers using the PICOS backend. The Python-MIP backend
    only supports the GUROBI and CBC solvers.

    Note that in some cases, the Python-MIP backend (GUROBI/CBC) might be faster
    setting up the problem than the PICOS backend since it has less overhead.

    Attributes:
        GUROBI_PYMIP (str): Recommended solver for most problems (LP, MIP).
            It uses the Python-MIP backend. Note that GUROBI is a commercial
            software which require a license. Free academic licenses are 
            also available at https://gurobi.com/free/.
        GUROBI (str): For using GUROBI with the PICOS backend instead of Python-MIP.
        COIN_OR_CBC (str): Free LP/MIP solver with good performance, provided with Python-MIP.
        CPLEX (str): Commercial solver with a performance similar to GUROBI. Only
            supported by the PICOS backend.
        SCIP (str): Free academic solver, supported by the PICOS backend.
        GLPK (str): Open source solver, supported by the PICOS backend.
        MOSEK (str): Commercial solver, supported by the PICOS backend.
    """
    GUROBI_PYMIP = "gurobi_pymip",
    GUROBI = "gurobi",
    COIN_OR_CBC = "cbc",
    CPLEX = "cplex",
    GLPK = "glpk",
    SCIP = "scip",
    CVXOPT = "cvxopt",
    MOSEK = "mosek"

class Status(str, Enum):
    """Subset of PICOS solver status codes.

    This subset of states offers compatibility between Python-MIP and PICOS.
    Status codes can be obtained after solving an optimization problem with
    the property `status`. For example:

    ```python
    >>> import miom
    >>> miom.load('@iMM1865').steady_state().set_rxn_objective('BIOMASS_reaction').solve().status
    {'status': 'optimal', 'objective_value': 798.811, 'elapsed_seconds': 0.479}
    ```
    """
    OPTIMAL = "optimal",
    FEASIBLE = "feasible"
    INFEASIBLE = "infeasible",
    UNBOUNDED = "unbounded",
    ILLPOSED = "illposed",
    EMPTY = "empty",
    PREMATURE = "premature",
    UNKNOWN = "unknown"