<commit_before># -*- coding: utf-8 -*-

# Copyright 2017 Novo Nordisk Foundation Center for Biosustainability,
# Technical University of Denmark.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""Supporting functions for stoichiometric consistency checks."""

from __future__ import absolute_import, division

import logging
import multiprocessing
from operator import attrgetter

import numpy as np
from cobra import Configuration, Reaction
from cobra.exceptions import Infeasible
from cobra.flux_analysis import flux_variability_analysis
from optlang.interface import INFEASIBLE, OPTIMAL
from optlang.symbolics import Zero

import memote.support.consistency_helpers as con_helpers
import memote.support.helpers as helpers


LOGGER = logging.getLogger(__name__)

# The following dictionary is based on the list of energy metabolites chosen
# as part of the following publication:
# <NAME>., <NAME>., <NAME>., <NAME>., & Lercher,
# <NAME>. (2017). Erroneous energy-generating cycles in published genome scale
# metabolic networks: Identification and removal. PLoS Computational
# Biology, 13(4), 1–14. http://doi.org/10.1371/journal.pcbi.1005494
ENERGY_COUPLES = {
    "MNXM3": "MNXM7",
    "MNXM63": "MNXM220",
    "MNXM51": "MNXM30",
    "MNXM121": "MNXM17",
    "MNXM423": "MNXM495",
    "MNXM6": "MNXM5",
    "MNXM10": "MNXM8",
    "MNXM38": "MNXM33",
    "MNXM208": "MNXM119",
    "MNXM191": "MNXM232",
    "MNXM223": "MNXM509",
    "MNXM7517": "MNXM12235",
    "MNXM12233": "MNXM12236",
    "MNXM558": "MNXM2178",
    "MNXM21": "MNXM12",
    "MNXM89557": "MNXM20",
}
TOLERANCE_THRESHOLD = 1e-07

cobra_configuration = Configuration()


def check_stoichiometric_consistency(model):
    """
    Verify the consistency of the model's stoichiometry.

    Parameters
    ----------
    model : cobra.Model
        The metabolic model under investigation.

    Notes
    -----
    See [1]_ section 3.1 for a complete description of the algorithm.

    .. [1] <NAME>., <NAME>, and <NAME>.
           "Detection of Stoichiometric Inconsistencies in Biomolecular
           Models."
           Bioinformatics 24, no. 19 (2008): 2245.

    """
    problem = model.problem
    # The transpose of the stoichiometric matrix N.T in the paper.
    stoich_trans = problem.Model()
    # We clone the existing configuration in order to apply non-default
    # settings, for example, for solver verbosity or timeout.
    stoich_trans.configuration = problem.Configuration.clone(
        config=model.solver.configuration, problem=stoich_trans
    )
    internal_rxns = con_helpers.get_internals(model)
    metabolites = set(met for rxn in internal_rxns for met in rxn.metabolites)
    LOGGER.info("model '%s' has %d internal reactions", model.id, len(internal_rxns))
    LOGGER.info("model '%s' has %d internal metabolites", model.id, len(metabolites))

    stoich_trans.add([problem.Variable(m.id, lb=1) for m in metabolites])
    stoich_trans.update()
    con_helpers.add_reaction_constraints(
        stoich_trans, internal_rxns, problem.Constraint
    )
    # The objective is to minimize the metabolite mass vector.
    stoich_trans.objective = problem.Objective(Zero, direction="min", sloppy=True)
    stoich_trans.objective.set_linear_coefficients(
        {var: 1.0 for var in stoich_trans.variables}
    )
    status = stoich_trans.optimize()
    if status == OPTIMAL:
        return True
    elif status == INFEASIBLE:
        return False
    else:
        raise RuntimeError(
            "Could not determine stoichiometric consistencty."
            " Solver status is '{}'"
            " (only optimal or infeasible expected).".format(status)
        )


def find_unconserved_metabolites(model):
    """
    Detect unconserved metabolites.

    Parameters
    ----------
    model : cobra.Model
        The metabolic model under investigation.

    Notes
    -----
    See [1]_ section 3.2 for a complete description of the algorithm.


    .. [1] <NAME>., <NAME>, and <NAME>.
           "Detection of Stoichiometric Inconsistencies in Biomolecular
           Models."
           Bioinformatics 24, no. 19 (2008): 2245.

    """
    problem = model.problem
    stoich_trans = problem.Model()
    # We clone the existing configuration in order to apply non-default
    # settings, for example, for solver verbosity or timeout.
    stoich_trans.configuration = problem.Configuration.clone(
        config=model.solver.configuration, problem=stoich_trans
    )
    internal_rxns = con_helpers.get_internals(model)
    metabolites = set(met for rxn in internal_rxns for met in rxn.metabolites)
    # The binary variables k[i] in the paper.
    k_vars = list()
    for met in metabolites:
        # The element m[i] of the mass vector.
        m_var = problem.Variable(met.id)
        k_var = problem.Variable("k_{}".format(met.id), type="binary")
        k_vars.append(k_var)
        stoich_trans.add([m_var, k_var])
        # This constraint is equivalent to 0 <= k[i] <= m[i].
        stoich_trans.add(
            problem.Constraint(k_var - m_var, ub=0, name="switch_{}".format(met.id))
        )
    stoich_trans.update()
    con_helpers.add_reaction_constraints(
        stoich_trans, internal_rxns, problem.Constraint
    )
    # The objective is to maximize the binary indicators k[i], subject to the
    # above inequality constraints.
    stoich_trans.objective = problem.Objective(Zero, sloppy=True, direction="max")
    stoich_trans.objective.set_linear_coefficients({var: 1.0 for var in k_vars})
    status = stoich_trans.optimize()
    if status == OPTIMAL:
        # TODO: See if that could be a Boolean test `bool(var.primal)`.
        return set(
            [
                model.metabolites.get_by_id(var.name[2:])
                for var in k_vars
                if var.primal < 0.8
            ]
        )
    else:
        raise RuntimeError(
            "Could not compute list of unconserved metabolites."
            " Solver status is '{}' (only optimal expected).".format(status)
        )


def find_inconsistent_min_stoichiometry(model, atol=1e-13, max_mets_computed=100):<commit_msg>fix: lower the soft cap to ten<commit_after># -*- coding: utf-8 -*-

# Copyright 2017 Novo Nordisk Foundation Center for Biosustainability,
# Technical University of Denmark.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""Supporting functions for stoichiometric consistency checks."""

from __future__ import absolute_import, division

import logging
import multiprocessing
from operator import attrgetter

import numpy as np
from cobra import Configuration, Reaction
from cobra.exceptions import Infeasible
from cobra.flux_analysis import flux_variability_analysis
from optlang.interface import INFEASIBLE, OPTIMAL
from optlang.symbolics import Zero

import memote.support.consistency_helpers as con_helpers
import memote.support.helpers as helpers


LOGGER = logging.getLogger(__name__)

# The following dictionary is based on the list of energy metabolites chosen
# as part of the following publication:
# <NAME>., <NAME>., <NAME>., <NAME>., & Lercher,
# <NAME>. (2017). Erroneous energy-generating cycles in published genome scale
# metabolic networks: Identification and removal. PLoS Computational
# Biology, 13(4), 1–14. http://doi.org/10.1371/journal.pcbi.1005494
ENERGY_COUPLES = {
    "MNXM3": "MNXM7",
    "MNXM63": "MNXM220",
    "MNXM51": "MNXM30",
    "MNXM121": "MNXM17",
    "MNXM423": "MNXM495",
    "MNXM6": "MNXM5",
    "MNXM10": "MNXM8",
    "MNXM38": "MNXM33",
    "MNXM208": "MNXM119",
    "MNXM191": "MNXM232",
    "MNXM223": "MNXM509",
    "MNXM7517": "MNXM12235",
    "MNXM12233": "MNXM12236",
    "MNXM558": "MNXM2178",
    "MNXM21": "MNXM12",
    "MNXM89557": "MNXM20",
}
TOLERANCE_THRESHOLD = 1e-07

cobra_configuration = Configuration()


def check_stoichiometric_consistency(model):
    """
    Verify the consistency of the model's stoichiometry.

    Parameters
    ----------
    model : cobra.Model
        The metabolic model under investigation.

    Notes
    -----
    See [1]_ section 3.1 for a complete description of the algorithm.

    .. [1] <NAME>., <NAME>, and <NAME>.
           "Detection of Stoichiometric Inconsistencies in Biomolecular
           Models."
           Bioinformatics 24, no. 19 (2008): 2245.

    """
    problem = model.problem
    # The transpose of the stoichiometric matrix N.T in the paper.
    stoich_trans = problem.Model()
    # We clone the existing configuration in order to apply non-default
    # settings, for example, for solver verbosity or timeout.
    stoich_trans.configuration = problem.Configuration.clone(
        config=model.solver.configuration, problem=stoich_trans
    )
    internal_rxns = con_helpers.get_internals(model)
    metabolites = set(met for rxn in internal_rxns for met in rxn.metabolites)
    LOGGER.info("model '%s' has %d internal reactions", model.id, len(internal_rxns))
    LOGGER.info("model '%s' has %d internal metabolites", model.id, len(metabolites))

    stoich_trans.add([problem.Variable(m.id, lb=1) for m in metabolites])
    stoich_trans.update()
    con_helpers.add_reaction_constraints(
        stoich_trans, internal_rxns, problem.Constraint
    )
    # The objective is to minimize the metabolite mass vector.
    stoich_trans.objective = problem.Objective(Zero, direction="min", sloppy=True)
    stoich_trans.objective.set_linear_coefficients(
        {var: 1.0 for var in stoich_trans.variables}
    )
    status = stoich_trans.optimize()
    if status == OPTIMAL:
        return True
    elif status == INFEASIBLE:
        return False
    else:
        raise RuntimeError(
            "Could not determine stoichiometric consistencty."
            " Solver status is '{}'"
            " (only optimal or infeasible expected).".format(status)
        )


def find_unconserved_metabolites(model):
    """
    Detect unconserved metabolites.

    Parameters
    ----------
    model : cobra.Model
        The metabolic model under investigation.

    Notes
    -----
    See [1]_ section 3.2 for a complete description of the algorithm.


    .. [1] <NAME>., <NAME>, and <NAME>.
           "Detection of Stoichiometric Inconsistencies in Biomolecular
           Models."
           Bioinformatics 24, no. 19 (2008): 2245.

    """
    problem = model.problem
    stoich_trans = problem.Model()
    # We clone the existing configuration in order to apply non-default
    # settings, for example, for solver verbosity or timeout.
    stoich_trans.configuration = problem.Configuration.clone(
        config=model.solver.configuration, problem=stoich_trans
    )
    internal_rxns = con_helpers.get_internals(model)
    metabolites = set(met for rxn in internal_rxns for met in rxn.metabolites)
    # The binary variables k[i] in the paper.
    k_vars = list()
    for met in metabolites:
        # The element m[i] of the mass vector.
        m_var = problem.Variable(met.id)
        k_var = problem.Variable("k_{}".format(met.id), type="binary")
        k_vars.append(k_var)
        stoich_trans.add([m_var, k_var])
        # This constraint is equivalent to 0 <= k[i] <= m[i].
        stoich_trans.add(
            problem.Constraint(k_var - m_var, ub=0, name="switch_{}".format(met.id))
        )
    stoich_trans.update()
    con_helpers.add_reaction_constraints(
        stoich_trans, internal_rxns, problem.Constraint
    )
    # The objective is to maximize the binary indicators k[i], subject to the
    # above inequality constraints.
    stoich_trans.objective = problem.Objective(Zero, sloppy=True, direction="max")
    stoich_trans.objective.set_linear_coefficients({var: 1.0 for var in k_vars})
    status = stoich_trans.optimize()
    if status == OPTIMAL:
        # TODO: See if that could be a Boolean test `bool(var.primal)`.
        return set(
            [
                model.metabolites.get_by_id(var.name[2:])
                for var in k_vars
                if var.primal < 0.8
            ]
        )
    else:
        raise RuntimeError(
            "Could not compute list of unconserved metabolites."
            " Solver status is '{}' (only optimal expected).".format(status)
        )


def find_inconsistent_min_stoichiometry(model, atol=1e-13, max_mets_computed=10):
{'fraction_non_alphanumeric': 0.10037764350453172, 'fraction_numerical': 0.025453172205438065, 'mean_word_length': 3.077918078226055, 'pattern_counts': {'":': 32, '<': 31, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5196202', 'n_tokens_mistral': 4401, 'n_tokens_neox': 3931, 'n_words': 1195}