# Copyright 2021 <NAME>

# 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.

"""Read and write SBML directly with `geckopy.Model`s.

Proteins are Species, members of Group Protein, with:
    - initialAmount: Concentration.
    - boundary_condition: True

The MW are initialAssingments.
"""

import datetime
import logging
import numbers
import re
import warnings
from collections import defaultdict
from math import isnan
from pathlib import Path

import libsbml
from cobra import Configuration
from cobra.core.gene import parse_gpr
from cobra.io.sbml import (
    BOUND_MINUS_INF,
    BOUND_PLUS_INF,
    F_GENE,
    F_GENE_REV,
    F_GROUP,
    F_GROUP_REV,
    F_REACTION,
    F_REACTION_REV,
    F_REPLACE,
    F_SPECIE,
    F_SPECIE_REV,
    LONG_SHORT_DIRECTION,
    LOWER_BOUND_ID,
    SBO_DEFAULT_FLUX_BOUND,
    SBO_EXCHANGE_REACTION,
    SBO_FBA_FRAMEWORK,
    SHORT_LONG_DIRECTION,
    UNITS_FLUX,
    UPPER_BOUND_ID,
    ZERO_BOUND_ID,
    CobraSBMLError,
    Gene,
    Group,
    Metabolite,
    _check,
    _check_required,
    _create_bound,
    _create_parameter,
    _get_doc_from_filename,
    _parse_annotations,
    _parse_notes_dict,
    _sbase_annotations,
    _sbase_notes_dict,
    linear_reaction_coefficients,
)
from cobra.util.solver import set_objective

from geckopy.model import Model
from geckopy.protein import Protein
from geckopy.reaction import Reaction

from .standard import group_proteins


__all__ = ["read_sbml_ec_model", "write_sbml_ec_model"]

LOGGER = logging.getLogger(__name__)
config = Configuration()

PROT_PATTERN = re.compile(
    r"prot_[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}"
)
PROT_EX_PATTERN = re.compile(
    r"prot_[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}_exchange"
)


def read_sbml_ec_model(
    filename: str,
    number: float = float,
    # re.Pattern does not exist in py36 so this type hint cannot be added now
    f_replace=F_REPLACE,
    set_missing_bounds: bool = False,
    hardcoded_rev_reactions: bool = True,
    **kwargs,
) -> Model:
    """Create `geckopy.Model` from SBMLDocument.

    Parameters
    ----------
    filename: str
    number: data type of stoichiometry: {float, int}
        In which data type should the stoichiometry be parsed.
    f_replace : dict of replacement functions for id replacement
    set_missing_bounds : flag to set missing bounds
    hardcoded_rev_reactions: bool
        if reversible reaction to account for proteins being consumed on both
        directions are written explicitly for in the SBML

    Returns
    -------
    cobra.core.Model

    """
    try:
        fsanitized = str(filename) if isinstance(filename, Path) else filename
        doc = _get_doc_from_filename(fsanitized)
    except IOError as e:
        raise e
    except Exception as original_error:
        raise CobraSBMLError(
            "Something went wrong reading the SBML model. Most likely the SBML"
            " model is not valid. Please check that your model is valid using "
            "the `cobra.io.sbml.validate_sbml_model` function or via the "
            "online validator at http://sbml.org/validator .\n"
            "\t`(model, errors) = validate_sbml_model(filename)`"
            "\nIf the model is valid and cannot be read please open an issue "
            f"at https://github.com/opencobra/cobrapy/issues: {original_error}"
        )

    if f_replace is None:
        f_replace = {}

    # SBML model
    model: libsbml.Model = doc.getModel()
    if model is None:
        raise CobraSBMLError("No SBML model detected in file.")
    model_fbc: libsbml.FbcModelPlugin = model.getPlugin("fbc")

    if not model_fbc:
        LOGGER.warning("Model does not contain SBML fbc package information.")
    else:
        if not model_fbc.isSetStrict():
            LOGGER.warning('Loading SBML model without fbc:strict="true"')

        # fbc-v1 (legacy)
        doc_fbc = doc.getPlugin("fbc")  # type: libsbml.FbcSBMLDocumentPlugin
        fbc_version = doc_fbc.getPackageVersion()

        if fbc_version == 1:
            LOGGER.warning(
                "Loading SBML with fbc-v1 (models should be encoded" " using fbc-v2)"
            )
            conversion_properties = libsbml.ConversionProperties()
            conversion_properties.addOption(
                "convert fbc v1 to fbc v2", True, "Convert FBC-v1 model to FBC-v2"
            )
            result = doc.convert(conversion_properties)
            if result != libsbml.LIBSBML_OPERATION_SUCCESS:
                raise Exception("Conversion of SBML fbc v1 to fbc v2 failed")

    # Model
    model_id = model.getIdAttribute()
    if not libsbml.SyntaxChecker.isValidSBMLSId(model_id):
        LOGGER.error("'%s' is not a valid SBML 'SId'." % model_id)
    geckopy_model = Model(model_id, hardcoded_rev_reactions=hardcoded_rev_reactions)
    geckopy_model.name = model.getName()

    # meta information
    meta = {
        "model.id": model_id,
        "level": model.getLevel(),
        "version": model.getVersion(),
        "packages": [],
    }
    # History
    creators = []
    created = None
    if model.isSetModelHistory():
        history = model.getModelHistory()  # type: libsbml.ModelHistory

        if history.isSetCreatedDate():
            created = history.getCreatedDate()

        for c in history.getListCreators():  # type: libsbml.ModelCreator
            creators.append(
                {
                    "familyName": c.getFamilyName() if c.isSetFamilyName() else None,
                    "givenName": c.getGivenName() if c.isSetGivenName() else None,
                    "organisation": c.getOrganisation()
                    if c.isSetOrganisation()
                    else None,
                    "email": c.getEmail() if c.isSetEmail() else None,
                }
            )

    meta["creators"] = creators
    meta["created"] = created
    meta["notes"] = _parse_notes_dict(doc)
    meta["annotation"] = _parse_annotations(doc)

    info = "<{}> SBML L{}V{}".format(model_id, model.getLevel(), model.getVersion())
    packages = {}
    for k in range(doc.getNumPlugins()):
        plugin = doc.getPlugin(k)  # type:libsbml.SBasePlugin
        key, value = plugin.getPackageName(), plugin.getPackageVersion()
        packages[key] = value
        info += ", {}-v{}".format(key, value)
        if key not in ["fbc", "groups", "l3v2extendedmath"]:
            LOGGER.warning(
                "SBML package '%s' not supported by cobrapy, "
                "information is not parsed",
                key,
            )
    meta["info"] = info
    meta["packages"] = packages
    geckopy_model._sbml = meta

    # notes and annotations
    geckopy_model.notes = _parse_notes_dict(model)
    geckopy_model.annotation = _parse_annotations(model)

    # Compartments
    # FIXME: update with new compartments
    compartments = {}
    for (
        compartment
    ) in model.getListOfCompartments():  # noqa: E501 type: libsbml.Compartment
        cid = _check_required(compartment, compartment.getIdAttribute(), "id")
        compartments[cid] = compartment.getName()
    geckopy_model.compartments = compartments

    # Species
    metabolites = []
    # proteins that rely on the naming convention "prot_UNIPROT_ID" will be
    # catched here. Those who are annotated by groups membership will be parsed
    # after the groups are processed.
    proteins = []
    boundary_metabolites = []
    if model.getNumSpecies() == 0:
        LOGGER.warning("No metabolites in model")

    for specie in model.getListOfSpecies():  # type: libsbml.Species
        sid = _check_required(specie, specie.getIdAttribute(), "id")
        if f_replace and F_SPECIE in f_replace:
            sid = f_replace[F_SPECIE](sid)

        met = Metabolite(sid)
        met.name = specie.getName()
        met.notes = _parse_notes_dict(specie)
        met.annotation = _parse_annotations(specie)
        met.compartment = specie.getCompartment()
        initial_amount = specie.getInitialAmount()

        specie_fbc = specie.getPlugin("fbc")  # type: libsbml.FbcSpeciesPlugin
        if specie_fbc:
            met.charge = specie_fbc.getCharge()
            met.formula = specie_fbc.getChemicalFormula()
        else:
            if specie.isSetCharge():
                LOGGER.warning(
                    "Use of the species charge attribute is "
                    "discouraged, use fbc:charge "
                    "instead: %s",
                    specie,
                )
                met.charge = specie.getCharge()
            else:
                if "CHARGE" in met.notes:
                    LOGGER.warning(
                        "Use of CHARGE in the notes element is "
                        "discouraged, use fbc:charge "
                        "instead: %s",
                        specie,
                    )
                    try:
                        met.charge = int(met.notes["CHARGE"])
                    except ValueError:
                        # handle nan, na, NA, ...
                        pass

            if "FORMULA" in met.notes:
                LOGGER.warning(
                    "Use of FORMULA in the notes element is "
                    "discouraged, use fbc:chemicalFormula "
                    "instead: %s",
                    specie,
                )
                met.formula = met.notes["FORMULA"]

        # Detect boundary metabolites
        if specie.getBoundaryCondition() is True:
            boundary_metabolites.append(met)

        if not PROT_PATTERN.match(met.id):
            metabolites.append(met)
        else:
            proteins.append(Protein(met, concentration=initial_amount))

    geckopy_model.add_metabolites(metabolites)
    geckopy_model.add_proteins(proteins)

    # Add exchange reactions for boundary metabolites
    ex_reactions = []
    for met in boundary_metabolites:
        ex_rid = "EX_{}".format(met.id)
        ex_reaction = Reaction(ex_rid)
        ex_reaction.name = ex_rid
        ex_reaction.annotation = {"sbo": SBO_EXCHANGE_REACTION}
        ex_reaction.lower_bound = config.lower_bound
        ex_reaction.upper_bound = config.upper_bound
        LOGGER.warning(
            "Adding exchange reaction %s with default bounds "
            "for boundary metabolite: %s." % (ex_reaction.id, met.id)
        )
        # species is reactant
        ex_reaction.add_metabolites({met: -1})
        ex_reactions.append(ex_reaction)
    geckopy_model.add_reactions(ex_reactions)

    # Genes
    if model_fbc:
        for (
            gp
        ) in model_fbc.getListOfGeneProducts():  # noqa: E501 type: libsbml.GeneProduct
            gid = _check_required(gp, gp.getIdAttribute(), "id")
            if f_replace and F_GENE in f_replace:
                gid = f_replace[F_GENE](gid)
            cobra_gene = Gene(gid)
            cobra_gene.name = gp.getName()
            if cobra_gene.name is None:
                cobra_gene.name = gid
            cobra_gene.annotation = _parse_annotations(gp)
            cobra_gene.notes = _parse_notes_dict(gp)

            geckopy_model.genes.append(cobra_gene)
    else:
        for (
            cobra_reaction
        ) in model.getListOfReactions():  # noqa: E501 type: libsbml.Reaction
            # fallback to notes information
            notes = _parse_notes_dict(cobra_reaction)
            if "GENE ASSOCIATION" in notes:
                gpr = notes["GENE ASSOCIATION"]
            elif "GENE_ASSOCIATION" in notes:
                gpr = notes["GENE_ASSOCIATION"]
            else:
                gpr = ""

            if len(gpr) > 0:
                gpr = gpr.replace("(", ";")
                gpr = gpr.replace(")", ";")
                gpr = gpr.replace("or", ";")
                gpr = gpr.replace("and", ";")
                # Interaction of the above replacements can lead to multiple
                # ;, which results in empty gids
                gids = [t.strip() for t in gpr.split(";")]
                gids = set(gids).difference({""})

                # create missing genes
                for gid in gids:
                    if f_replace and F_GENE in f_replace:
                        gid = f_replace[F_GENE](gid)

                    if gid not in geckopy_model.genes:
                        cobra_gene = Gene(gid)
                        cobra_gene.name = gid
                        geckopy_model.genes.append(cobra_gene)

    # GPR rules
    def process_association(ass):
        """Recursively convert gpr association to a gpr string.

        Defined as inline functions to not pass the replacement dict around.
        """
        if ass.isFbcOr():
            return " ".join(
                [
                    "(",
                    " or ".join(
                        process_association(c) for c in ass.getListOfAssociations()
                    ),
                    ")",
                ]
            )
        elif ass.isFbcAnd():
            return " ".join(
                [
                    "(",
                    " and ".join(
                        process_association(c) for c in ass.getListOfAssociations()
                    ),
                    ")",
                ]
            )
        elif ass.isGeneProductRef():
            gid = ass.getGeneProduct()
            if f_replace and F_GENE in f_replace:
                return f_replace[F_GENE](gid)
            else:
                return gid

    # Reactions
    missing_bounds = False
    reactions = []
    if model.getNumReactions() == 0:
        LOGGER.warning("No reactions in model")

    for reaction in model.getListOfReactions():  # type: libsbml.Reaction
        rid = _check_required(reaction, reaction.getIdAttribute(), "id")
        # proteins are parsed based on Species, prot exchanges are ignored
        if PROT_EX_PATTERN.search(rid):
            continue
        if f_replace and F_REACTION in f_replace:
            rid = f_replace[F_REACTION](rid)
        cobra_reaction = Reaction(rid)
        cobra_reaction.name = reaction.getName()
        cobra_reaction.annotation = _parse_annotations(reaction)
        cobra_reaction.notes = _parse_notes_dict(reaction)

        # set bounds
        p_ub, p_lb = None, None
        r_fbc = reaction.getPlugin("fbc")  # type: libsbml.FbcReactionPlugin
        if r_fbc:
            # bounds in fbc
            lb_id = r_fbc.getLowerFluxBound()
            if lb_id:
                p_lb = model.getParameter(lb_id)  # type: libsbml.Parameter
                if p_lb and p_lb.getConstant() and (p_lb.getValue() is not None):
                    cobra_reaction.lower_bound = p_lb.getValue()
                else:
                    raise CobraSBMLError(
                        "No constant bound '%s' for " "reaction: %s" % (p_lb, reaction)
                    )

            ub_id = r_fbc.getUpperFluxBound()
            if ub_id:
                p_ub = model.getParameter(ub_id)  # type: libsbml.Parameter
                if p_ub and p_ub.getConstant() and (p_ub.getValue() is not None):
                    cobra_reaction.upper_bound = p_ub.getValue()
                else:
                    raise CobraSBMLError(
                        "No constant bound '%s' for " "reaction: %s" % (p_ub, reaction)
                    )

        elif reaction.isSetKineticLaw():
            # some legacy models encode bounds in kinetic laws
            klaw = reaction.getKineticLaw()  # type: libsbml.KineticLaw
            p_lb = klaw.getParameter(
                "LOWER_BOUND"
            )  # noqa: E501 type: libsbml.LocalParameter
            if p_lb:
                cobra_reaction.lower_bound = p_lb.getValue()
            p_ub = klaw.getParameter(
                "UPPER_BOUND"
            )  # noqa: E501 type: libsbml.LocalParameter
            if p_ub:
                cobra_reaction.upper_bound = p_ub.getValue()

            if p_ub is not None or p_lb is not None:
                LOGGER.warning(
                    "Encoding LOWER_BOUND and UPPER_BOUND in "
                    "KineticLaw is discouraged, "
                    "use fbc:fluxBounds instead: %s",
                    reaction,
                )

        if p_lb is None:
            missing_bounds = True
            lower_bound = config.lower_bound
            cobra_reaction.lower_bound = lower_bound
            LOGGER.warning(
                "Missing lower flux bound set to '%s' for " " reaction: '%s'",
                lower_bound,
                reaction,
            )

        if p_ub is None:
            missing_bounds = True
            upper_bound = config.upper_bound
            cobra_reaction.upper_bound = upper_bound
            LOGGER.warning(
                "Missing upper flux bound set to '%s' for " " reaction: '%s'",
                upper_bound,
                reaction,
            )

        # add reaction
        reactions.append(cobra_reaction)

        # parse equation
        stoichiometry = defaultdict(lambda: 0)
        for (
            sref
        ) in reaction.getListOfReactants():  # noqa: E501 type: libsbml.SpeciesReference
            sid = _check_required(sref, sref.getSpecies(), "species")

            if f_replace and F_SPECIE in f_replace:
                sid = f_replace[F_SPECIE](sid)
            stoichiometry[sid] -= number(
                _check_required(sref, sref.getStoichiometry(), "stoichiometry")
            )

        for (
            sref
        ) in reaction.getListOfProducts():  # noqa: E501 type: libsbml.SpeciesReference
            sid = _check_required(sref, sref.getSpecies(), "species")

            if f_replace and F_SPECIE in f_replace:
                sid = f_replace[F_SPECIE](sid)
            stoichiometry[sid] += number(
                _check_required(sref, sref.getStoichiometry(), "stoichiometry")
            )

        # convert to metabolite objects
        object_stoichiometry = {}
        for met_id in stoichiometry:
            target_set = (
                geckopy_model.proteins
                if met_id in geckopy_model.proteins
                else geckopy_model.metabolites
            )
            metabolite = target_set.get_by_id(met_id)
            object_stoichiometry[metabolite] = stoichiometry[met_id]
        cobra_reaction.add_metabolites(object_stoichiometry)

        # GPR
        if r_fbc:
            gpr = ""
            gpa = (
                r_fbc.getGeneProductAssociation()
            )  # noqa: E501 type: libsbml.GeneProductAssociation
            if gpa is not None:
                association = (
                    gpa.getAssociation()
                )  # noqa: E501 type: libsbml.FbcAssociation
                gpr = process_association(association)
        else:
            # fallback to notes information
            notes = cobra_reaction.notes
            if "GENE ASSOCIATION" in notes:
                gpr = notes["GENE ASSOCIATION"]
            elif "GENE_ASSOCIATION" in notes:
                gpr = notes["GENE_ASSOCIATION"]
            else:
                gpr = ""

            if len(gpr) > 0:
                LOGGER.warning(
                    "Use of GENE ASSOCIATION or GENE_ASSOCIATION "
                    "in the notes element is discouraged, use "
                    "fbc:gpr instead: %s",
                    reaction,
                )
                if f_replace and F_GENE in f_replace:
                    gpr = " ".join(f_replace[F_GENE](t) for t in gpr.split(" "))

        # remove outside parenthesis, if any
        if gpr.startswith("(") and gpr.endswith(")"):
            try:
                parse_gpr(gpr[1:-1].strip())
                gpr = gpr[1:-1].strip()
            except (SyntaxError, TypeError) as e:
                LOGGER.warning(
                    f"Removing parenthesis from gpr {gpr} leads to "
                    f"an error, so keeping parenthesis, error: {e}",
                )

        cobra_reaction.gene_reaction_rule = gpr

    geckopy_model.add_reactions(reactions)

    # Objective
    obj_direction = "max"
    coefficients = {}
    if model_fbc:
        obj_list = (
            model_fbc.getListOfObjectives()
        )  # noqa: E501 type: libsbml.ListOfObjectives
        if obj_list is None:
            LOGGER.warning("listOfObjectives element not found")
        elif obj_list.size() == 0:
            LOGGER.warning("No objective in listOfObjectives")
        elif not obj_list.getActiveObjective():
            LOGGER.warning("No active objective in listOfObjectives")
        else:
            obj_id = obj_list.getActiveObjective()
            obj = model_fbc.getObjective(obj_id)  # type: libsbml.Objective
            obj_direction = LONG_SHORT_DIRECTION[obj.getType()]

            for (
                flux_obj
            ) in (
                obj.getListOfFluxObjectives()
            ):  # noqa: E501 type: libsbml.FluxObjective
                rid = flux_obj.getReaction()
                if f_replace and F_REACTION in f_replace:
                    rid = f_replace[F_REACTION](rid)
                try:
                    objective_reaction = geckopy_model.reactions.get_by_id(rid)
                except KeyError:
                    raise CobraSBMLError("Objective reaction '%s' " "not found" % rid)
                try:
                    coefficients[objective_reaction] = number(flux_obj.getCoefficient())
                except ValueError as e:
                    LOGGER.warning(str(e))
    else:
        # some legacy models encode objective coefficients in kinetic laws
        for reaction in model.getListOfReactions():  # type: libsbml.Reaction
            if reaction.isSetKineticLaw():
                klaw = reaction.getKineticLaw()  # type: libsbml.KineticLaw
                p_oc = klaw.getParameter(
                    "OBJECTIVE_COEFFICIENT"
                )  # type: libsbml.LocalParameter
                if p_oc:
                    rid = _check_required(reaction, reaction.getIdAttribute(), "id")
                    if f_replace and F_REACTION in f_replace:
                        rid = f_replace[F_REACTION](rid)
                    try:
                        objective_reaction = geckopy_model.reactions.get_by_id(rid)
                    except KeyError:
                        raise CobraSBMLError(
                            "Objective reaction '%s' " "not found", rid
                        )
                    try:
                        coefficients[objective_reaction] = number(p_oc.getValue())
                    except ValueError as e:
                        LOGGER.warning(str(e))

                    LOGGER.warning(
                        "Encoding OBJECTIVE_COEFFICIENT in "
                        "KineticLaw is discouraged, "
                        "use fbc:fluxObjective "
                        "instead: %s",
                        reaction,
                    )

    if len(coefficients) == 0:
        LOGGER.error(
            "No objective coefficients in model. Unclear what should " "be optimized"
        )
    set_objective(geckopy_model, coefficients)
    geckopy_model.solver.objective.direction = obj_direction

    # parse groups
    model_groups = model.getPlugin("groups")  # type: libsbml.GroupsModelPlugin
    groups = []
    if model_groups:
        # calculate hashmaps to lookup objects in O(1)
        sid_map = {}
        metaid_map = {}
        for obj_list in [
            model.getListOfCompartments(),
            model.getListOfSpecies(),
            model.getListOfReactions(),
            model_groups.getListOfGroups(),
        ]:

            for sbase in obj_list:  # type: libsbml.SBase
                if sbase.isSetId():
                    sid_map[sbase.getIdAttribute()] = sbase
                if sbase.isSetMetaId():
                    metaid_map[sbase.getMetaId()] = sbase

        # create groups
        for group in model_groups.getListOfGroups():  # type: libsbml.Group
            gid = _check_required(group, group.getIdAttribute(), "id")
            if f_replace and F_GROUP in f_replace:
                gid = f_replace[F_GROUP](gid)
            cobra_group = Group(gid)
            cobra_group.name = group.getName()
            if group.isSetKind():
                cobra_group.kind = group.getKindAsString()
            cobra_group.annotation = _parse_annotations(group)
            cobra_group.notes = _parse_notes_dict(group)

            cobra_members = []
            for member in group.getListOfMembers():  # type: libsbml.Member
                if member.isSetIdRef():
                    obj = sid_map[member.getIdRef()]
                elif member.isSetMetaIdRef():
                    obj = metaid_map[member.getMetaIdRef()]

                typecode = obj.getTypeCode()
                obj_id = _check_required(obj, obj.getIdAttribute(), "id")

                # id replacements
                cobra_member = None
                if typecode == libsbml.SBML_SPECIES:
                    if f_replace and F_SPECIE in f_replace:
                        obj_id = f_replace[F_SPECIE](obj_id)
                    try:
                        cobra_member = geckopy_model.metabolites.get_by_id(obj_id)
                    except KeyError:
                        cobra_member = geckopy_model.proteins.get_by_id(obj_id)
                elif typecode == libsbml.SBML_REACTION:
                    if f_replace and F_REACTION in f_replace:
                        obj_id = f_replace[F_REACTION](obj_id)
                    cobra_member = geckopy_model.reactions.get_by_id(obj_id)
                elif typecode == libsbml.SBML_FBC_GENEPRODUCT:
                    if f_replace and F_GENE in f_replace:
                        obj_id = f_replace[F_GENE](obj_id)
                    cobra_member = geckopy_model.genes.get_by_id(obj_id)
                else:
                    LOGGER.warning(
                        "Member %s could not be added to group %s."
                        "unsupported type code: "
                        "%s" % (member, group, typecode)
                    )

                if cobra_member:
                    cobra_members.append(cobra_member)

            cobra_group.add_members(cobra_members)
            groups.append(cobra_group)
    else:
        # parse deprecated subsystems on reactions
        groups_dict = {}
        for cobra_reaction in geckopy_model.reactions:
            if "SUBSYSTEM" in cobra_reaction.notes:
                g_name = cobra_reaction.notes["SUBSYSTEM"]
                if g_name in groups_dict:
                    groups_dict[g_name].append(cobra_reaction)
                else:
                    groups_dict[g_name] = [cobra_reaction]

        for gid, cobra_members in groups_dict.items():
            if f_replace and F_GROUP in f_replace:
                gid = f_replace[F_GROUP](gid)
            cobra_group = Group(gid, name=gid, kind="collection")
            cobra_group.add_members(cobra_members)
            groups.append(cobra_group)

    geckopy_model.add_groups(groups)

    # now add everything under group Proteins to model.proteins if it was not
    # already added based on naming conventions
    if geckopy_model.groups.query("Protein"):
        g_proteins = geckopy_model.groups.Protein.members.copy()
        g_proteins = {
            prot: {reac: reac.metabolites[prot] for reac in prot.reactions}
            for prot in g_proteins
            if prot not in geckopy_model.proteins
        }

        if g_proteins:
            geckopy_model.remove_metabolites(g_proteins.keys())
            geckopy_model.add_proteins([Protein(prot) for prot in g_proteins])
            for prot, reactions in g_proteins.items():
                for reac, stoich in reactions.items():
                    # reverse the (negative) stoichiometry coefficient to kcat
                    # we expect that Proteins that are identified only by their
                    # group are respecting the specification and do form part
                    # of reactions
                    # TODO: provide options to tune this behvior
                    reac.add_protein(prot.id, -1 / (stoich * 3600))

    # general hint for missing flux bounds
    if missing_bounds:
        LOGGER.warning(
            "Missing flux bounds on reactions set to default bounds."
            "As best practise and to avoid confusion flux bounds "
            "should be set explicitly on all reactions."
        )

    return geckopy_model


def write_sbml_ec_model(
    ec_model: Model,
    filename: str,
    f_replace=F_REPLACE,
    units=True,
    group_untyped_proteins: bool = True,
):
    """Write cobra model to filename.

    Enzyme constraint changes: proteins are written as metabolites with
    initialAmount.

    The created model is SBML level 3 version 1 (L1V3) with
    fbc package v2 (fbc-v2).
    If the given filename ends with the suffix ".gz" (for example,
    "myfile.xml.gz"), libSBML assumes the caller wants the file to be
    written compressed in gzip format. Similarly, if the given filename
    ends with ".zip" or ".bz2", libSBML assumes the caller wants the
    file to be compressed in zip or bzip2 format (respectively). Files
    whose names lack these suffixes will be written uncompressed. Special
    considerations for the zip format: If the given filename ends with
    ".zip", the file placed in the zip archive will have the suffix
    ".xml" or ".sbml".  For example, the file in the zip archive will
    be named "test.xml" if the given filename is "test.xml.zip" or
    "test.zip". Similarly, the filename in the archive will be
    "test.sbml" if the given filename is "test.sbml.zip".

    Warnings
    --------
    If :meth:`geckopy.Model.constrain_pool()` was called on the model, the
    protein pool (exchange and metabolite) will be ignored and NOT written
    to the SBML file. A `UserWarning` appears in this case.

    Parameters
    ----------
    cobra_model : geckopy.Model
        Model instance which is written to SBML
    filename : string
        path to which the model is written
    f_replace: dict of replacement functions for id replacement
    group_untyped_proteins: bool
        if True (default) the proteins whose id is not compliant with Uniprot
        and which are not part of the :code:`Protein` group will be added to it.
        This is an inplace operation!
    """
    cobra_model = ec_model
    protein_pool_metabolite = None
    protein_pool_exchange = None
    if hasattr(ec_model, "protein_pool_exchange"):
        warnings.warn("Protein pool will not be serialized to the SBML document.")
        protein_pool_metabolite = ec_model.common_protein_pool
        protein_pool_exchange = ec_model.protein_pool_exchange
    if f_replace is None:
        f_replace = {}
    if group_untyped_proteins:
        prot_to_add = [
            prot.id for prot in cobra_model.proteins if not PROT_PATTERN.match(prot.id)
        ]
        group_proteins(ec_model, prot_to_add)

    sbml_ns = libsbml.SBMLNamespaces(3, 1)  # SBML L3V1
    sbml_ns.addPackageNamespace("fbc", 2)  # fbc-v2

    doc: libsbml.SBMLDocument = libsbml.SBMLDocument(sbml_ns)
    doc.setPackageRequired("fbc", False)
    doc.setSBOTerm(SBO_FBA_FRAMEWORK)

    model: libsbml.Model = doc.createModel()
    model_fbc: libsbml.FbcModelPlugin = model.getPlugin("fbc")
    model_fbc.setStrict(True)

    if cobra_model.id is not None:
        model.setId(cobra_model.id)
        model.setMetaId("meta_" + cobra_model.id)
    else:
        model.setMetaId("meta_model")
    if cobra_model.name is not None:
        model.setName(cobra_model.name)

    # for parsing annotation corresponding to the model
    _sbase_annotations(model, cobra_model.annotation)
    # for parsing notes corresponding to the model
    _sbase_notes_dict(model, cobra_model.notes)

    # Meta information (ModelHistory) related to SBMLDocument
    if hasattr(cobra_model, "_sbml"):
        meta = cobra_model._sbml
        if "annotation" in meta:
            _sbase_annotations(doc, meta["annotation"])
        if "notes" in meta:
            _sbase_notes_dict(doc, meta["notes"])

        history: libsbml.ModelHistory = libsbml.ModelHistory()
        if "created" in meta and meta["created"]:
            history.setCreatedDate(meta["created"])
        else:
            time = datetime.datetime.now()
            timestr = time.strftime("%Y-%m-%dT%H:%M:%S")
            date = libsbml.Date(timestr)
            _check(history.setCreatedDate(date), "set creation date")
            _check(history.setModifiedDate(date), "set modified date")

        if "creators" in meta:
            for cobra_creator in meta["creators"]:  # noqa: E501
                creator: libsbml.ModelCreator = libsbml.ModelCreator()
                if cobra_creator.get("familyName", None):
                    creator.setFamilyName(cobra_creator["familyName"])
                if cobra_creator.get("givenName", None):
                    creator.setGivenName(cobra_creator["givenName"])
                if cobra_creator.get("organisation", None):
                    creator.setOrganisation(cobra_creator["organisation"])
                if cobra_creator.get("email", None):
                    creator.setEmail(cobra_creator["email"])

                _check(history.addCreator(creator), "adding creator to ModelHistory.")

        # TODO: Will be implemented as part of
        #  https://github.com/opencobra/cobrapy/issues/810
        # _check(model.setModelHistory(history), 'set model history')

    # Units
    if units:
        flux_udef = model.createUnitDefinition()  # type:libsbml.UnitDefinition
        flux_udef.setId(UNITS_FLUX[0])
        for u in UNITS_FLUX[1]:
            unit = flux_udef.createUnit()  # type: libsbml.Unit
            unit.setKind(u.kind)
            unit.setExponent(u.exponent)
            unit.setScale(u.scale)
            unit.setMultiplier(u.multiplier)

    # minimum and maximum value from model
    if len(cobra_model.reactions) > 0:
        min_value = min(cobra_model.reactions.list_attr("lower_bound"))
        max_value = max(cobra_model.reactions.list_attr("upper_bound"))
    else:
        min_value = config.lower_bound
        max_value = config.upper_bound

    _create_parameter(
        model, pid=LOWER_BOUND_ID, value=min_value, sbo=SBO_DEFAULT_FLUX_BOUND
    )
    _create_parameter(
        model, pid=UPPER_BOUND_ID, value=max_value, sbo=SBO_DEFAULT_FLUX_BOUND
    )
    _create_parameter(model, pid=ZERO_BOUND_ID, value=0, sbo=SBO_DEFAULT_FLUX_BOUND)
    _create_parameter(
        model, pid=BOUND_MINUS_INF, value=-float("Inf"), sbo=SBO_DEFAULT_FLUX_BOUND
    )
    _create_parameter(
        model, pid=BOUND_PLUS_INF, value=float("Inf"), sbo=SBO_DEFAULT_FLUX_BOUND
    )

    # Compartments
    for cid, name in cobra_model.compartments.items():
        compartment: libsbml.Compartment = model.createCompartment()
        compartment.setId(cid)
        compartment.setName(name)
        compartment.setConstant(True)

    # Species
    for metabolite in cobra_model.metabolites:
        if metabolite == protein_pool_metabolite:
            continue
        specie: libsbml.Species = model.createSpecies()
        specie.setId(
            f_replace[F_SPECIE_REV](metabolite.id)
            if f_replace and F_SPECIE_REV in f_replace
            else metabolite.id
        )
        specie.setConstant(False)
        specie.setBoundaryCondition(False)
        specie.setHasOnlySubstanceUnits(False)
        specie.setName(metabolite.name)
        specie.setCompartment(metabolite.compartment)
        s_fbc: libsbml.FbcSpeciesPlugin = specie.getPlugin("fbc")
        if metabolite.charge is not None:
            s_fbc.setCharge(int(metabolite.charge))
        if metabolite.formula is not None:
            s_fbc.setChemicalFormula(metabolite.formula)

        _sbase_annotations(specie, metabolite.annotation)
        _sbase_notes_dict(specie, metabolite.notes)

    for metabolite in ec_model.proteins:
        specie: libsbml.Species = model.createSpecies()
        specie.setId(
            f_replace[F_SPECIE_REV](metabolite.id)
            if f_replace and F_SPECIE_REV in f_replace
            else metabolite.id
        )
        specie.setConstant(False)
        specie.setBoundaryCondition(False)
        specie.setHasOnlySubstanceUnits(False)
        specie.setName(metabolite.name)
        specie.setCompartment(metabolite.compartment)
        if isinstance(metabolite.concentration, numbers.Number) and not isnan(
            metabolite.concentration
        ):
            specie.setInitialAmount(metabolite.concentration)
        s_fbc: libsbml.FbcSpeciesPlugin = specie.getPlugin("fbc")
        if metabolite.charge is not None:
            s_fbc.setCharge(int(metabolite.charge))
        if metabolite.formula is not None:
            s_fbc.setChemicalFormula(metabolite.formula)

        _sbase_annotations(specie, metabolite.annotation)

    # Genes
    for cobra_gene in cobra_model.genes:
        gp: libsbml.GeneProduct = model_fbc.createGeneProduct()
        gid = cobra_gene.id
        if f_replace and F_GENE_REV in f_replace:
            gid = f_replace[F_GENE_REV](gid)
        gp.setId(gid)
        gname = cobra_gene.name
        if gname is None or len(gname) == 0:
            gname = gid
        gp.setName(gname)
        gp.setLabel(gid)

        _sbase_annotations(gp, cobra_gene.annotation)
        _sbase_notes_dict(gp, cobra_gene.notes)

    # Objective
    objective: libsbml.Objective = model_fbc.createObjective()
    objective.setId("obj")
    objective.setType(SHORT_LONG_DIRECTION[cobra_model.objective.direction])
    model_fbc.setActiveObjectiveId("obj")

    # Reactions
    reaction_coefficients = linear_reaction_coefficients(cobra_model)
    for cobra_reaction in cobra_model.reactions:
        if cobra_reaction == protein_pool_exchange:
            continue
        rid = cobra_reaction.id
        if f_replace and F_REACTION_REV in f_replace:
            rid = f_replace[F_REACTION_REV](rid)
        reaction: libsbml.Reaction = model.createReaction()
        reaction.setId(rid)
        reaction.setName(cobra_reaction.name)
        reaction.setFast(False)
        reaction.setReversible((cobra_reaction.lower_bound < 0))
        _sbase_annotations(reaction, cobra_reaction.annotation)
        _sbase_notes_dict(reaction, cobra_reaction.notes)

        # stoichiometry
        for metabolite, stoichiometry in cobra_reaction._metabolites.items():
            sid = metabolite.id
            if f_replace and F_SPECIE_REV in f_replace:
                sid = f_replace[F_SPECIE_REV](sid)
            if stoichiometry < 0:
                sref = (
                    reaction.createReactant()
                )  # noqa: E501 type: libsbml.SpeciesReference
                sref.setSpecies(sid)
                sref.setStoichiometry(-stoichiometry)
                sref.setConstant(True)
            else:
                sref = (
                    reaction.createProduct()
                )  # noqa: E501 type: libsbml.SpeciesReference
                sref.setSpecies(sid)
                sref.setStoichiometry(stoichiometry)
                sref.setConstant(True)

        # bounds
        r_fbc: libsbml.FbcReactionPlugin = reaction.getPlugin("fbc")
        r_fbc.setLowerFluxBound(
            _create_bound(
                model,
                cobra_reaction,
                "lower_bound",
                f_replace=f_replace,
                units=units,
                flux_udef=flux_udef,
            )
        )
        r_fbc.setUpperFluxBound(
            _create_bound(
                model,
                cobra_reaction,
                "upper_bound",
                f_replace=f_replace,
                units=units,
                flux_udef=flux_udef,
            )
        )

        # GPR
        gpr = cobra_reaction.gene_reaction_rule
        if gpr is not None and len(gpr) > 0:

            # replace ids in string
            if f_replace and F_GENE_REV in f_replace:
                gpr = gpr.replace("(", "( ")
                gpr = gpr.replace(")", " )")
                tokens = gpr.split()

                for k in range(len(tokens)):
                    if tokens[k] not in ["and", "or", "(", ")"]:
                        tokens[k] = f_replace[F_GENE_REV](tokens[k])
                gpr_new = " ".join(tokens)
            else:
                gpr_new = gpr

            gpa: libsbml.GeneProductAssociation = (
                r_fbc.createGeneProductAssociation()
            )  # noqa: E501
            # uses ids to identify GeneProducts (True),
            # does not create GeneProducts (False)
            _check(gpa.setAssociation(gpr_new, True, False), "set gpr: " + gpr_new)

        # objective coefficients
        if reaction_coefficients.get(cobra_reaction, 0) != 0:
            flux_obj: libsbml.FluxObjective = (
                objective.createFluxObjective()
            )  # noqa: E501
            flux_obj.setReaction(rid)
            flux_obj.setCoefficient(cobra_reaction.objective_coefficient)

    # write groups
    if len(cobra_model.groups) > 0:
        doc.enablePackage(
            "http://www.sbml.org/sbml/level3/version1/groups/version1", "groups", True
        )
        doc.setPackageRequired("groups", False)
        model_group: libsbml.GroupsModelPlugin = model.getPlugin("groups")  # noqa: E501
        for cobra_group in cobra_model.groups:
            group: libsbml.Group = model_group.createGroup()
            if f_replace and F_GROUP_REV in f_replace:
                gid = f_replace[F_GROUP_REV](cobra_group.id)
            else:
                gid = cobra_group.id
            group.setId(gid)
            group.setName(cobra_group.name)
            group.setKind(cobra_group.kind)

            _sbase_notes_dict(group, cobra_group.notes)
            _sbase_annotations(group, cobra_group.annotation)

            for cobra_member in cobra_group.members:
                member: libsbml.Member = group.createMember()
                mid = cobra_member.id
                m_type = str(type(cobra_member))

                # id replacements
                if "Reaction" in m_type:
                    if f_replace and F_REACTION_REV in f_replace:
                        mid = f_replace[F_REACTION_REV](mid)
                if "Metabolite" in m_type or "Protein" in m_type:
                    if f_replace and F_SPECIE_REV in f_replace:
                        mid = f_replace[F_SPECIE_REV](mid)
                if "Gene" in m_type:
                    if f_replace and F_GENE_REV in f_replace:
                        mid = f_replace[F_GENE_REV](mid)

                member.setIdRef(mid)
                if cobra_member.name and len(cobra_member.name) > 0:
                    member.setName(cobra_member.name)

    if isinstance(filename, str):
        # write to path
        libsbml.writeSBMLToFile(doc, filename)

    elif hasattr(filename, "write"):
        # write to file handle
        sbml_str = libsbml.writeSBMLToString(doc)
        filename.write(sbml_str)

# Copyright 2021 <NAME>

# 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.

"""Read and write SBML directly with `geckopy.Model`s.

Proteins are Species, members of Group Protein, with:
    - initialAmount: Concentration.
    - boundary_condition: True

The MW are initialAssingments.
"""

import datetime
import logging
import numbers
import re
import warnings
from collections import defaultdict
from math import isnan
from pathlib import Path

import libsbml
from cobra import Configuration
from cobra.core.gene import parse_gpr
from cobra.io.sbml import (
    BOUND_MINUS_INF,
    BOUND_PLUS_INF,
    F_GENE,
    F_GENE_REV,
    F_GROUP,
    F_GROUP_REV,
    F_REACTION,
    F_REACTION_REV,
    F_REPLACE,
    F_SPECIE,
    F_SPECIE_REV,
    LONG_SHORT_DIRECTION,
    LOWER_BOUND_ID,
    SBO_DEFAULT_FLUX_BOUND,
    SBO_EXCHANGE_REACTION,
    SBO_FBA_FRAMEWORK,
    SHORT_LONG_DIRECTION,
    UNITS_FLUX,
    UPPER_BOUND_ID,
    ZERO_BOUND_ID,
    CobraSBMLError,
    Gene,
    Group,
    Metabolite,
    _check,
    _check_required,
    _create_bound,
    _create_parameter,
    _get_doc_from_filename,
    _parse_annotations,
    _parse_notes_dict,
    _sbase_annotations,
    _sbase_notes_dict,
    linear_reaction_coefficients,
)
from cobra.util.solver import set_objective

from geckopy.model import Model
from geckopy.protein import Protein
from geckopy.reaction import Reaction

from .standard import group_proteins


__all__ = ["read_sbml_ec_model", "write_sbml_ec_model"]

LOGGER = logging.getLogger(__name__)
config = Configuration()

PROT_PATTERN = re.compile(
    r"prot_[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}"
)
PROT_EX_PATTERN = re.compile(
    r"prot_[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}_exchange"
)


def read_sbml_ec_model(
    filename: str,
    number: float = float,
    # re.Pattern does not exist in py36 so this type hint cannot be added now
    f_replace=F_REPLACE,
    set_missing_bounds: bool = False,
    hardcoded_rev_reactions: bool = True,
    **kwargs,
) -> Model:
    """Create `geckopy.Model` from SBMLDocument.

    Parameters
    ----------
    filename: str
    number: data type of stoichiometry: {float, int}
        In which data type should the stoichiometry be parsed.
    f_replace : dict of replacement functions for id replacement
    set_missing_bounds : flag to set missing bounds
    hardcoded_rev_reactions: bool
        if reversible reaction to account for proteins being consumed on both
        directions are written explicitly for in the SBML

    Returns
    -------
    cobra.core.Model

    """
    try:
        fsanitized = str(filename) if isinstance(filename, Path) else filename
        doc = _get_doc_from_filename(fsanitized)
    except IOError as e:
        raise e
    except Exception as original_error:
        raise CobraSBMLError(
            "Something went wrong reading the SBML model. Most likely the SBML"
            " model is not valid. Please check that your model is valid using "
            "the `cobra.io.sbml.validate_sbml_model` function or via the "
            "online validator at http://sbml.org/validator .\n"
            "\t`(model, errors) = validate_sbml_model(filename)`"
            "\nIf the model is valid and cannot be read please open an issue "
            f"at https://github.com/opencobra/cobrapy/issues: {original_error}"
        )

    if f_replace is None:
        f_replace = {}

    # SBML model
    model: libsbml.Model = doc.getModel()
    if model is None:
        raise CobraSBMLError("No SBML model detected in file.")
    model_fbc: libsbml.FbcModelPlugin = model.getPlugin("fbc")

    if not model_fbc:
        LOGGER.warning("Model does not contain SBML fbc package information.")
    else:
        if not model_fbc.isSetStrict():
            LOGGER.warning('Loading SBML model without fbc:strict="true"')

        # fbc-v1 (legacy)
        doc_fbc = doc.getPlugin("fbc")  # type: libsbml.FbcSBMLDocumentPlugin
        fbc_version = doc_fbc.getPackageVersion()

        if fbc_version == 1:
            LOGGER.warning(
                "Loading SBML with fbc-v1 (models should be encoded" " using fbc-v2)"
            )
            conversion_properties = libsbml.ConversionProperties()
            conversion_properties.addOption(
                "convert fbc v1 to fbc v2", True, "Convert FBC-v1 model to FBC-v2"
            )
            result = doc.convert(conversion_properties)
            if result != libsbml.LIBSBML_OPERATION_SUCCESS:
                raise Exception("Conversion of SBML fbc v1 to fbc v2 failed")

    # Model
    model_id = model.getIdAttribute()
    if not libsbml.SyntaxChecker.isValidSBMLSId(model_id):
        LOGGER.error("'%s' is not a valid SBML 'SId'." % model_id)
    geckopy_model = Model(model_id, hardcoded_rev_reactions=hardcoded_rev_reactions)
    geckopy_model.name = model.getName()

    # meta information
    meta = {
        "model.id": model_id,
        "level": model.getLevel(),
        "version": model.getVersion(),
        "packages": [],
    }
    # History
    creators = []
    created = None
    if model.isSetModelHistory():
        history = model.getModelHistory()  # type: libsbml.ModelHistory

        if history.isSetCreatedDate():
            created = history.getCreatedDate()

        for c in history.getListCreators():  # type: libsbml.ModelCreator
            creators.append(
                {
                    "familyName": c.getFamilyName() if c.isSetFamilyName() else None,
                    "givenName": c.getGivenName() if c.isSetGivenName() else None,
                    "organisation": c.getOrganisation()
                    if c.isSetOrganisation()
                    else None,
                    "email": c.getEmail() if c.isSetEmail() else None,
                }
            )

    meta["creators"] = creators
    meta["created"] = created
    meta["notes"] = _parse_notes_dict(doc)
    meta["annotation"] = _parse_annotations(doc)

    info = "<{}> SBML L{}V{}".format(model_id, model.getLevel(), model.getVersion())
    packages = {}
    for k in range(doc.getNumPlugins()):
        plugin = doc.getPlugin(k)  # type:libsbml.SBasePlugin
        key, value = plugin.getPackageName(), plugin.getPackageVersion()
        packages[key] = value
        info += ", {}-v{}".format(key, value)
        if key not in ["fbc", "groups", "l3v2extendedmath"]:
            LOGGER.warning(
                "SBML package '%s' not supported by cobrapy, "
                "information is not parsed",
                key,
            )
    meta["info"] = info
    meta["packages"] = packages
    geckopy_model._sbml = meta

    # notes and annotations
    geckopy_model.notes = _parse_notes_dict(model)
    geckopy_model.annotation = _parse_annotations(model)

    # Compartments
    # FIXME: update with new compartments
    compartments = {}
    for (
        compartment
    ) in model.getListOfCompartments():  # noqa: E501 type: libsbml.Compartment
        cid = _check_required(compartment, compartment.getIdAttribute(), "id")
        compartments[cid] = compartment.getName()
    geckopy_model.compartments = compartments

    # Species
    metabolites = []
    # proteins that rely on the naming convention "prot_UNIPROT_ID" will be
    # catched here. Those who are annotated by groups membership will be parsed
    # after the groups are processed.
    proteins = []
    boundary_metabolites = []
    if model.getNumSpecies() == 0:
        LOGGER.warning("No metabolites in model")

    for specie in model.getListOfSpecies():  # type: libsbml.Species
        sid = _check_required(specie, specie.getIdAttribute(), "id")
        if f_replace and F_SPECIE in f_replace:
            sid = f_replace[F_SPECIE](sid)

        met = Metabolite(sid)
        met.name = specie.getName()
        met.notes = _parse_notes_dict(specie)
        met.annotation = _parse_annotations(specie)
        met.compartment = specie.getCompartment()
        initial_amount = specie.getInitialAmount()

        specie_fbc = specie.getPlugin("fbc")  # type: libsbml.FbcSpeciesPlugin
        if specie_fbc:
            met.charge = specie_fbc.getCharge()
            met.formula = specie_fbc.getChemicalFormula()
        else:
            if specie.isSetCharge():
                LOGGER.warning(
                    "Use of the species charge attribute is "
                    "discouraged, use fbc:charge "
                    "instead: %s",
                    specie,
                )
                met.charge = specie.getCharge()
            else:
                if "CHARGE" in met.notes:
                    LOGGER.warning(
                        "Use of CHARGE in the notes element is "
                        "discouraged, use fbc:charge "
                        "instead: %s",
                        specie,
                    )
                    try:
                        met.charge = int(met.notes["CHARGE"])
                    except ValueError:
                        # handle nan, na, NA, ...
                        pass

            if "FORMULA" in met.notes:
                LOGGER.warning(
                    "Use of FORMULA in the notes element is "
                    "discouraged, use fbc:chemicalFormula "
                    "instead: %s",
                    specie,
                )
                met.formula = met.notes["FORMULA"]

        # Detect boundary metabolites
        if specie.getBoundaryCondition() is True:
            boundary_metabolites.append(met)

        if not PROT_PATTERN.match(met.id):
            metabolites.append(met)
        else:
            proteins.append(Protein(met, concentration=initial_amount))

    geckopy_model.add_metabolites(metabolites)
    geckopy_model.add_proteins(proteins)

    # Add exchange reactions for boundary metabolites
    ex_reactions = []
    for met in boundary_metabolites:
        ex_rid = "EX_{}".format(met.id)
        ex_reaction = Reaction(ex_rid)
        ex_reaction.name = ex_rid
        ex_reaction.annotation = {"sbo": SBO_EXCHANGE_REACTION}
        ex_reaction.lower_bound = config.lower_bound
        ex_reaction.upper_bound = config.upper_bound
        LOGGER.warning(
            "Adding exchange reaction %s with default bounds "
            "for boundary metabolite: %s." % (ex_reaction.id, met.id)
        )
        # species is reactant
        ex_reaction.add_metabolites({met: -1})
        ex_reactions.append(ex_reaction)
    geckopy_model.add_reactions(ex_reactions)

    # Genes
    if model_fbc:
        for (
            gp
        ) in model_fbc.getListOfGeneProducts():  # noqa: E501 type: libsbml.GeneProduct
            gid = _check_required(gp, gp.getIdAttribute(), "id")
            if f_replace and F_GENE in f_replace:
                gid = f_replace[F_GENE](gid)
            cobra_gene = Gene(gid)
            cobra_gene.name = gp.getName()
            if cobra_gene.name is None:
                cobra_gene.name = gid
            cobra_gene.annotation = _parse_annotations(gp)
            cobra_gene.notes = _parse_notes_dict(gp)

            geckopy_model.genes.append(cobra_gene)
    else:
        for (
            cobra_reaction
        ) in model.getListOfReactions():  # noqa: E501 type: libsbml.Reaction
            # fallback to notes information
            notes = _parse_notes_dict(cobra_reaction)
            if "GENE ASSOCIATION" in notes:
                gpr = notes["GENE ASSOCIATION"]
            elif "GENE_ASSOCIATION" in notes:
                gpr = notes["GENE_ASSOCIATION"]
            else:
                gpr = ""

            if len(gpr) > 0:
                gpr = gpr.replace("(", ";")
                gpr = gpr.replace(")", ";")
                gpr = gpr.replace("or", ";")
                gpr = gpr.replace("and", ";")
                # Interaction of the above replacements can lead to multiple
                # ;, which results in empty gids
                gids = [t.strip() for t in gpr.split(";")]
                gids = set(gids).difference({""})

                # create missing genes
                for gid in gids:
                    if f_replace and F_GENE in f_replace:
                        gid = f_replace[F_GENE](gid)

                    if gid not in geckopy_model.genes:
                        cobra_gene = Gene(gid)
                        cobra_gene.name = gid
                        geckopy_model.genes.append(cobra_gene)

    # GPR rules
    def process_association(ass):
        """Recursively convert gpr association to a gpr string.

        Defined as inline functions to not pass the replacement dict around.
        """
        if ass.isFbcOr():
            return " ".join(
                [
                    "(",
                    " or ".join(
                        process_association(c) for c in ass.getListOfAssociations()
                    ),
                    ")",
                ]
            )
        elif ass.isFbcAnd():
            return " ".join(
                [
                    "(",
                    " and ".join(
                        process_association(c) for c in ass.getListOfAssociations()
                    ),
                    ")",
                ]
            )
        elif ass.isGeneProductRef():
            gid = ass.getGeneProduct()
            if f_replace and F_GENE in f_replace:
                return f_replace[F_GENE](gid)
            else:
                return gid

    # Reactions
    missing_bounds = False
    reactions = []
    if model.getNumReactions() == 0:
        LOGGER.warning("No reactions in model")

    for reaction in model.getListOfReactions():  # type: libsbml.Reaction
        rid = _check_required(reaction, reaction.getIdAttribute(), "id")
        # proteins are parsed based on Species, prot exchanges are ignored
        if PROT_EX_PATTERN.search(rid):
            continue
        if f_replace and F_REACTION in f_replace:
            rid = f_replace[F_REACTION](rid)
        cobra_reaction = Reaction(rid)
        cobra_reaction.name = reaction.getName()
        cobra_reaction.annotation = _parse_annotations(reaction)
        cobra_reaction.notes = _parse_notes_dict(reaction)

        # set bounds
        p_ub, p_lb = None, None
        r_fbc = reaction.getPlugin("fbc")  # type: libsbml.FbcReactionPlugin
        if r_fbc:
            # bounds in fbc
            lb_id = r_fbc.getLowerFluxBound()
            if lb_id:
                p_lb = model.getParameter(lb_id)  # type: libsbml.Parameter
                if p_lb and p_lb.getConstant() and (p_lb.getValue() is not None):
                    cobra_reaction.lower_bound = p_lb.getValue()
                else:
                    raise CobraSBMLError(
                        "No constant bound '%s' for " "reaction: %s" % (p_lb, reaction)
                    )

            ub_id = r_fbc.getUpperFluxBound()
            if ub_id:
                p_ub = model.getParameter(ub_id)  # type: libsbml.Parameter
                if p_ub and p_ub.getConstant() and (p_ub.getValue() is not None):
                    cobra_reaction.upper_bound = p_ub.getValue()
                else:
                    raise CobraSBMLError(
                        "No constant bound '%s' for " "reaction: %s" % (p_ub, reaction)
                    )

        elif reaction.isSetKineticLaw():
            # some legacy models encode bounds in kinetic laws
            klaw = reaction.getKineticLaw()  # type: libsbml.KineticLaw
            p_lb = klaw.getParameter(
                "LOWER_BOUND"
            )  # noqa: E501 type: libsbml.LocalParameter
            if p_lb:
                cobra_reaction.lower_bound = p_lb.getValue()
            p_ub = klaw.getParameter(
                "UPPER_BOUND"
            )  # noqa: E501 type: libsbml.LocalParameter
            if p_ub:
                cobra_reaction.upper_bound = p_ub.getValue()

            if p_ub is not None or p_lb is not None:
                LOGGER.warning(
                    "Encoding LOWER_BOUND and UPPER_BOUND in "
                    "KineticLaw is discouraged, "
                    "use fbc:fluxBounds instead: %s",
                    reaction,
                )

        if p_lb is None:
            missing_bounds = True
            lower_bound = config.lower_bound
            cobra_reaction.lower_bound = lower_bound
            LOGGER.warning(
                "Missing lower flux bound set to '%s' for " " reaction: '%s'",
                lower_bound,
                reaction,
            )

        if p_ub is None:
            missing_bounds = True
            upper_bound = config.upper_bound
            cobra_reaction.upper_bound = upper_bound
            LOGGER.warning(
                "Missing upper flux bound set to '%s' for " " reaction: '%s'",
                upper_bound,
                reaction,
            )

        # add reaction
        reactions.append(cobra_reaction)

        # parse equation
        stoichiometry = defaultdict(lambda: 0)
        for (
            sref
        ) in reaction.getListOfReactants():  # noqa: E501 type: libsbml.SpeciesReference
            sid = _check_required(sref, sref.getSpecies(), "species")

            if f_replace and F_SPECIE in f_replace:
                sid = f_replace[F_SPECIE](sid)
            stoichiometry[sid] -= number(
                _check_required(sref, sref.getStoichiometry(), "stoichiometry")
            )

        for (
            sref
        ) in reaction.getListOfProducts():  # noqa: E501 type: libsbml.SpeciesReference
            sid = _check_required(sref, sref.getSpecies(), "species")

            if f_replace and F_SPECIE in f_replace:
                sid = f_replace[F_SPECIE](sid)
            stoichiometry[sid] += number(
                _check_required(sref, sref.getStoichiometry(), "stoichiometry")
            )

        # convert to metabolite objects
        object_stoichiometry = {}
        for met_id in stoichiometry:
            target_set = (
                geckopy_model.proteins
                if met_id in geckopy_model.proteins
                else geckopy_model.metabolites
            )
            metabolite = target_set.get_by_id(met_id)
            object_stoichiometry[metabolite] = stoichiometry[met_id]
        cobra_reaction.add_metabolites(object_stoichiometry)

        # GPR
        if r_fbc:
            gpr = ""
            gpa = (
                r_fbc.getGeneProductAssociation()
            )  # noqa: E501 type: libsbml.GeneProductAssociation
            if gpa is not None:
                association = (
                    gpa.getAssociation()
                )  # noqa: E501 type: libsbml.FbcAssociation
                gpr = process_association(association)
        else:
            # fallback to notes information
            notes = cobra_reaction.notes
            if "GENE ASSOCIATION" in notes:
                gpr = notes["GENE ASSOCIATION"]
            elif "GENE_ASSOCIATION" in notes:
                gpr = notes["GENE_ASSOCIATION"]
            else:
                gpr = ""

            if len(gpr) > 0:
                LOGGER.warning(
                    "Use of GENE ASSOCIATION or GENE_ASSOCIATION "
                    "in the notes element is discouraged, use "
                    "fbc:gpr instead: %s",
                    reaction,
                )
                if f_replace and F_GENE in f_replace:
                    gpr = " ".join(f_replace[F_GENE](t) for t in gpr.split(" "))

        # remove outside parenthesis, if any
        if gpr.startswith("(") and gpr.endswith(")"):
            try:
                parse_gpr(gpr[1:-1].strip())
                gpr = gpr[1:-1].strip()
            except (SyntaxError, TypeError) as e:
                LOGGER.warning(
                    f"Removing parenthesis from gpr {gpr} leads to "
                    f"an error, so keeping parenthesis, error: {e}",
                )

        cobra_reaction.gene_reaction_rule = gpr

    geckopy_model.add_reactions(reactions)

    # Objective
    obj_direction = "max"
    coefficients = {}
    if model_fbc:
        obj_list = (
            model_fbc.getListOfObjectives()
        )  # noqa: E501 type: libsbml.ListOfObjectives
        if obj_list is None:
            LOGGER.warning("listOfObjectives element not found")
        elif obj_list.size() == 0:
            LOGGER.warning("No objective in listOfObjectives")
        elif not obj_list.getActiveObjective():
            LOGGER.warning("No active objective in listOfObjectives")
        else:
            obj_id = obj_list.getActiveObjective()
            obj = model_fbc.getObjective(obj_id)  # type: libsbml.Objective
            obj_direction = LONG_SHORT_DIRECTION[obj.getType()]

            for (
                flux_obj
            ) in (
                obj.getListOfFluxObjectives()
            ):  # noqa: E501 type: libsbml.FluxObjective
                rid = flux_obj.getReaction()
                if f_replace and F_REACTION in f_replace:
                    rid = f_replace[F_REACTION](rid)
                try:
                    objective_reaction = geckopy_model.reactions.get_by_id(rid)
                except KeyError:
                    raise CobraSBMLError("Objective reaction '%s' " "not found" % rid)
                try:
                    coefficients[objective_reaction] = number(flux_obj.getCoefficient())
                except ValueError as e:
                    LOGGER.warning(str(e))
    else:
        # some legacy models encode objective coefficients in kinetic laws
        for reaction in model.getListOfReactions():  # type: libsbml.Reaction
            if reaction.isSetKineticLaw():
                klaw = reaction.getKineticLaw()  # type: libsbml.KineticLaw
                p_oc = klaw.getParameter(
                    "OBJECTIVE_COEFFICIENT"
                )  # type: libsbml.LocalParameter
                if p_oc:
                    rid = _check_required(reaction, reaction.getIdAttribute(), "id")
                    if f_replace and F_REACTION in f_replace:
                        rid = f_replace[F_REACTION](rid)
                    try:
                        objective_reaction = geckopy_model.reactions.get_by_id(rid)
                    except KeyError:
                        raise CobraSBMLError(
                            "Objective reaction '%s' " "not found", rid
                        )
                    try:
                        coefficients[objective_reaction] = number(p_oc.getValue())
                    except ValueError as e:
                        LOGGER.warning(str(e))

                    LOGGER.warning(
                        "Encoding OBJECTIVE_COEFFICIENT in "
                        "KineticLaw is discouraged, "
                        "use fbc:fluxObjective "
                        "instead: %s",
                        reaction,
                    )

    if len(coefficients) == 0:
        LOGGER.error(
            "No objective coefficients in model. Unclear what should " "be optimized"
        )
    set_objective(geckopy_model, coefficients)
    geckopy_model.solver.objective.direction = obj_direction

    # parse groups
    model_groups = model.getPlugin("groups")  # type: libsbml.GroupsModelPlugin
    groups = []
    if model_groups:
        # calculate hashmaps to lookup objects in O(1)
        sid_map = {}
        metaid_map = {}
        for obj_list in [
            model.getListOfCompartments(),
            model.getListOfSpecies(),
            model.getListOfReactions(),
            model_groups.getListOfGroups(),
        ]:

            for sbase in obj_list:  # type: libsbml.SBase
                if sbase.isSetId():
                    sid_map[sbase.getIdAttribute()] = sbase
                if sbase.isSetMetaId():
                    metaid_map[sbase.getMetaId()] = sbase

        # create groups
        for group in model_groups.getListOfGroups():  # type: libsbml.Group
            gid = _check_required(group, group.getIdAttribute(), "id")
            if f_replace and F_GROUP in f_replace:
                gid = f_replace[F_GROUP](gid)
            cobra_group = Group(gid)
            cobra_group.name = group.getName()
            if group.isSetKind():
                cobra_group.kind = group.getKindAsString()
            cobra_group.annotation = _parse_annotations(group)
            cobra_group.notes = _parse_notes_dict(group)

            cobra_members = []
            for member in group.getListOfMembers():  # type: libsbml.Member
                if member.isSetIdRef():
                    obj = sid_map[member.getIdRef()]
                elif member.isSetMetaIdRef():
                    obj = metaid_map[member.getMetaIdRef()]

                typecode = obj.getTypeCode()
                obj_id = _check_required(obj, obj.getIdAttribute(), "id")

                # id replacements
                cobra_member = None
                if typecode == libsbml.SBML_SPECIES:
                    if f_replace and F_SPECIE in f_replace:
                        obj_id = f_replace[F_SPECIE](obj_id)
                    try:
                        cobra_member = geckopy_model.metabolites.get_by_id(obj_id)
                    except KeyError:
                        cobra_member = geckopy_model.proteins.get_by_id(obj_id)
                elif typecode == libsbml.SBML_REACTION:
                    if f_replace and F_REACTION in f_replace:
                        obj_id = f_replace[F_REACTION](obj_id)
                    cobra_member = geckopy_model.reactions.get_by_id(obj_id)
                elif typecode == libsbml.SBML_FBC_GENEPRODUCT:
                    if f_replace and F_GENE in f_replace:
                        obj_id = f_replace[F_GENE](obj_id)
                    cobra_member = geckopy_model.genes.get_by_id(obj_id)
                else:
                    LOGGER.warning(
                        "Member %s could not be added to group %s."
                        "unsupported type code: "
                        "%s" % (member, group, typecode)
                    )

                if cobra_member:
                    cobra_members.append(cobra_member)

            cobra_group.add_members(cobra_members)
            groups.append(cobra_group)
    else:
        # parse deprecated subsystems on reactions
        groups_dict = {}
        for cobra_reaction in geckopy_model.reactions:
            if "SUBSYSTEM" in cobra_reaction.notes:
                g_name = cobra_reaction.notes["SUBSYSTEM"]
                if g_name in groups_dict:
                    groups_dict[g_name].append(cobra_reaction)
                else:
                    groups_dict[g_name] = [cobra_reaction]

        for gid, cobra_members in groups_dict.items():
            if f_replace and F_GROUP in f_replace:
                gid = f_replace[F_GROUP](gid)
            cobra_group = Group(gid, name=gid, kind="collection")
            cobra_group.add_members(cobra_members)
            groups.append(cobra_group)

    geckopy_model.add_groups(groups)

    # now add everything under group Proteins to model.proteins if it was not
    # already added based on naming conventions
    if geckopy_model.groups.query("Protein"):
        g_proteins = geckopy_model.groups.Protein.members.copy()
        g_proteins = {
            prot: {reac: reac.metabolites[prot] for reac in prot.reactions}
            for prot in g_proteins
            if prot not in geckopy_model.proteins
        }

        if g_proteins:
            geckopy_model.remove_metabolites(g_proteins.keys())
            geckopy_model.add_proteins([Protein(prot) for prot in g_proteins])
            for prot, reactions in g_proteins.items():
                for reac, stoich in reactions.items():
                    # reverse the (negative) stoichiometry coefficient to kcat
                    # we expect that Proteins that are identified only by their
                    # group are respecting the specification and do form part
                    # of reactions
                    # TODO: provide options to tune this behvior
                    reac.add_protein(prot.id, -1 / (stoich * 3600))

    # general hint for missing flux bounds
    if missing_bounds:
        LOGGER.warning(
            "Missing flux bounds on reactions set to default bounds."
            "As best practise and to avoid confusion flux bounds "
            "should be set explicitly on all reactions."
        )

    return geckopy_model


def write_sbml_ec_model(
    ec_model: Model,
    filename: str,
    f_replace=F_REPLACE,
    units=True,
    group_untyped_proteins: bool = True,
):
    """Write cobra model to filename.

    Enzyme constraint changes: proteins are written as metabolites with
    initialAmount.

    The created model is SBML level 3 version 1 (L1V3) with
    fbc package v2 (fbc-v2).
    If the given filename ends with the suffix ".gz" (for example,
    "myfile.xml.gz"), libSBML assumes the caller wants the file to be
    written compressed in gzip format. Similarly, if the given filename
    ends with ".zip" or ".bz2", libSBML assumes the caller wants the
    file to be compressed in zip or bzip2 format (respectively). Files
    whose names lack these suffixes will be written uncompressed. Special
    considerations for the zip format: If the given filename ends with
    ".zip", the file placed in the zip archive will have the suffix
    ".xml" or ".sbml".  For example, the file in the zip archive will
    be named "test.xml" if the given filename is "test.xml.zip" or
    "test.zip". Similarly, the filename in the archive will be
    "test.sbml" if the given filename is "test.sbml.zip".

    Warnings
    --------
    If :meth:`geckopy.Model.constrain_pool()` was called on the model, the
    protein pool (exchange and metabolite) will be ignored and NOT written
    to the SBML file. A `UserWarning` appears in this case.

    Parameters
    ----------
    cobra_model : geckopy.Model
        Model instance which is written to SBML
    filename : string
        path to which the model is written
    f_replace: dict of replacement functions for id replacement
    group_untyped_proteins: bool
        if True (default) the proteins whose id is not compliant with Uniprot
        and which are not part of the :code:`Protein` group will be added to it.
        This is an inplace operation!
    """
    cobra_model = ec_model
    protein_pool_metabolite = None
    protein_pool_exchange = None
    if hasattr(ec_model, "protein_pool_exchange"):
        warnings.warn("Protein pool will not be serialized to the SBML document.")
        protein_pool_metabolite = ec_model.common_protein_pool
        protein_pool_exchange = ec_model.protein_pool_exchange
    if f_replace is None:
        f_replace = {}
    if group_untyped_proteins:
        prot_to_add = [
            prot.id for prot in cobra_model.proteins if not PROT_PATTERN.match(prot.id)
        ]
        group_proteins(ec_model, prot_to_add)

    sbml_ns = libsbml.SBMLNamespaces(3, 1)  # SBML L3V1
    sbml_ns.addPackageNamespace("fbc", 2)  # fbc-v2

    doc: libsbml.SBMLDocument = libsbml.SBMLDocument(sbml_ns)
    doc.setPackageRequired("fbc", False)
    doc.setSBOTerm(SBO_FBA_FRAMEWORK)

    model: libsbml.Model = doc.createModel()
    model_fbc: libsbml.FbcModelPlugin = model.getPlugin("fbc")
    model_fbc.setStrict(True)

    if cobra_model.id is not None:
        model.setId(cobra_model.id)
        model.setMetaId("meta_" + cobra_model.id)
    else:
        model.setMetaId("meta_model")
    if cobra_model.name is not None:
        model.setName(cobra_model.name)

    # for parsing annotation corresponding to the model
    _sbase_annotations(model, cobra_model.annotation)
    # for parsing notes corresponding to the model
    _sbase_notes_dict(model, cobra_model.notes)

    # Meta information (ModelHistory) related to SBMLDocument
    if hasattr(cobra_model, "_sbml"):
        meta = cobra_model._sbml
        if "annotation" in meta:
            _sbase_annotations(doc, meta["annotation"])
        if "notes" in meta:
            _sbase_notes_dict(doc, meta["notes"])

        history: libsbml.ModelHistory = libsbml.ModelHistory()
        if "created" in meta and meta["created"]:
            history.setCreatedDate(meta["created"])
        else:
            time = datetime.datetime.now()
            timestr = time.strftime("%Y-%m-%dT%H:%M:%S")
            date = libsbml.Date(timestr)
            _check(history.setCreatedDate(date), "set creation date")
            _check(history.setModifiedDate(date), "set modified date")

        if "creators" in meta:
            for cobra_creator in meta["creators"]:  # noqa: E501
                creator: libsbml.ModelCreator = libsbml.ModelCreator()
                if cobra_creator.get("familyName", None):
                    creator.setFamilyName(cobra_creator["familyName"])
                if cobra_creator.get("givenName", None):
                    creator.setGivenName(cobra_creator["givenName"])
                if cobra_creator.get("organisation", None):
                    creator.setOrganisation(cobra_creator["organisation"])
                if cobra_creator.get("email", None):
                    creator.setEmail(cobra_creator["email"])

                _check(history.addCreator(creator), "adding creator to ModelHistory.")

        # TODO: Will be implemented as part of
        #  https://github.com/opencobra/cobrapy/issues/810
        # _check(model.setModelHistory(history), 'set model history')

    # Units
    if units:
        flux_udef = model.createUnitDefinition()  # type:libsbml.UnitDefinition
        flux_udef.setId(UNITS_FLUX[0])
        for u in UNITS_FLUX[1]:
            unit = flux_udef.createUnit()  # type: libsbml.Unit
            unit.setKind(u.kind)
            unit.setExponent(u.exponent)
            unit.setScale(u.scale)
            unit.setMultiplier(u.multiplier)

    # minimum and maximum value from model
    if len(cobra_model.reactions) > 0:
        min_value = min(cobra_model.reactions.list_attr("lower_bound"))
        max_value = max(cobra_model.reactions.list_attr("upper_bound"))
    else:
        min_value = config.lower_bound
        max_value = config.upper_bound

    _create_parameter(
        model, pid=LOWER_BOUND_ID, value=min_value, sbo=SBO_DEFAULT_FLUX_BOUND
    )
    _create_parameter(
        model, pid=UPPER_BOUND_ID, value=max_value, sbo=SBO_DEFAULT_FLUX_BOUND
    )
    _create_parameter(model, pid=ZERO_BOUND_ID, value=0, sbo=SBO_DEFAULT_FLUX_BOUND)
    _create_parameter(
        model, pid=BOUND_MINUS_INF, value=-float("Inf"), sbo=SBO_DEFAULT_FLUX_BOUND
    )
    _create_parameter(
        model, pid=BOUND_PLUS_INF, value=float("Inf"), sbo=SBO_DEFAULT_FLUX_BOUND
    )

    # Compartments
    for cid, name in cobra_model.compartments.items():
        compartment: libsbml.Compartment = model.createCompartment()
        compartment.setId(cid)
        compartment.setName(name)
        compartment.setConstant(True)

    # Species
    for metabolite in cobra_model.metabolites:
        if metabolite == protein_pool_metabolite:
            continue
        specie: libsbml.Species = model.createSpecies()
        specie.setId(
            f_replace[F_SPECIE_REV](metabolite.id)
            if f_replace and F_SPECIE_REV in f_replace
            else metabolite.id
        )
        specie.setConstant(False)
        specie.setBoundaryCondition(False)
        specie.setHasOnlySubstanceUnits(False)
        specie.setName(metabolite.name)
        specie.setCompartment(metabolite.compartment)
        s_fbc: libsbml.FbcSpeciesPlugin = specie.getPlugin("fbc")
        if metabolite.charge is not None:
            s_fbc.setCharge(int(metabolite.charge))
        if metabolite.formula is not None:
            s_fbc.setChemicalFormula(metabolite.formula)

        _sbase_annotations(specie, metabolite.annotation)
        _sbase_notes_dict(specie, metabolite.notes)

    for metabolite in ec_model.proteins:
        specie: libsbml.Species = model.createSpecies()
        specie.setId(
            f_replace[F_SPECIE_REV](metabolite.id)
            if f_replace and F_SPECIE_REV in f_replace
            else metabolite.id
        )
        specie.setConstant(False)
        specie.setBoundaryCondition(False)
        specie.setHasOnlySubstanceUnits(False)
        specie.setName(metabolite.name)
        specie.setCompartment(metabolite.compartment)
        if isinstance(metabolite.concentration, numbers.Number) and not isnan(
            metabolite.concentration
        ):
            specie.setInitialAmount(metabolite.concentration)
        s_fbc: libsbml.FbcSpeciesPlugin = specie.getPlugin("fbc")
        if metabolite.charge is not None:
            s_fbc.setCharge(int(metabolite.charge))
        if metabolite.formula is not None:
            s_fbc.setChemicalFormula(metabolite.formula)

        _sbase_annotations(specie, metabolite.annotation)

    # Genes
    for cobra_gene in cobra_model.genes:
        gp: libsbml.GeneProduct = model_fbc.createGeneProduct()
        gid = cobra_gene.id
        if f_replace and F_GENE_REV in f_replace:
            gid = f_replace[F_GENE_REV](gid)
        gp.setId(gid)
        gname = cobra_gene.name
        if gname is None or len(gname) == 0:
            gname = gid
        gp.setName(gname)
        gp.setLabel(gid)

        _sbase_annotations(gp, cobra_gene.annotation)
        _sbase_notes_dict(gp, cobra_gene.notes)

    # Objective
    objective: libsbml.Objective = model_fbc.createObjective()
    objective.setId("obj")
    objective.setType(SHORT_LONG_DIRECTION[cobra_model.objective.direction])
    model_fbc.setActiveObjectiveId("obj")

    # Reactions
    reaction_coefficients = linear_reaction_coefficients(cobra_model)
    for cobra_reaction in cobra_model.reactions:
        if cobra_reaction == protein_pool_exchange:
            continue
        rid = cobra_reaction.id
        if f_replace and F_REACTION_REV in f_replace:
            rid = f_replace[F_REACTION_REV](rid)
        reaction: libsbml.Reaction = model.createReaction()
        reaction.setId(rid)
        reaction.setName(cobra_reaction.name)
        reaction.setFast(False)
        reaction.setReversible((cobra_reaction.lower_bound < 0))
        _sbase_annotations(reaction, cobra_reaction.annotation)
        _sbase_notes_dict(reaction, cobra_reaction.notes)

        # stoichiometry
        for metabolite, stoichiometry in cobra_reaction._metabolites.items():
            sid = metabolite.id
            if f_replace and F_SPECIE_REV in f_replace:
                sid = f_replace[F_SPECIE_REV](sid)
            if stoichiometry < 0:
                sref = (
                    reaction.createReactant()
                )  # noqa: E501 type: libsbml.SpeciesReference
                sref.setSpecies(sid)
                sref.setStoichiometry(-stoichiometry)
                sref.setConstant(True)
            else:
                sref = (
                    reaction.createProduct()
                )  # noqa: E501 type: libsbml.SpeciesReference
                sref.setSpecies(sid)
                sref.setStoichiometry(stoichiometry)
                sref.setConstant(True)

        # bounds
        r_fbc: libsbml.FbcReactionPlugin = reaction.getPlugin("fbc")
        r_fbc.setLowerFluxBound(
            _create_bound(
                model,
                cobra_reaction,
                "lower_bound",
                f_replace=f_replace,
                units=units,
                flux_udef=flux_udef,
            )
        )
        r_fbc.setUpperFluxBound(
            _create_bound(
                model,
                cobra_reaction,
                "upper_bound",
                f_replace=f_replace,
                units=units,
                flux_udef=flux_udef,
            )
        )

        # GPR
        gpr = cobra_reaction.gene_reaction_rule
        if gpr is not None and len(gpr) > 0:

            # replace ids in string
            if f_replace and F_GENE_REV in f_replace:
                gpr = gpr.replace("(", "( ")
                gpr = gpr.replace(")", " )")
                tokens = gpr.split()

                for k in range(len(tokens)):
                    if tokens[k] not in ["and", "or", "(", ")"]:
                        tokens[k] = f_replace[F_GENE_REV](tokens[k])
                gpr_new = " ".join(tokens)
            else:
                gpr_new = gpr

            gpa: libsbml.GeneProductAssociation = (
                r_fbc.createGeneProductAssociation()
            )  # noqa: E501
            # uses ids to identify GeneProducts (True),
            # does not create GeneProducts (False)
            _check(gpa.setAssociation(gpr_new, True, False), "set gpr: " + gpr_new)

        # objective coefficients
        if reaction_coefficients.get(cobra_reaction, 0) != 0:
            flux_obj: libsbml.FluxObjective = (
                objective.createFluxObjective()
            )  # noqa: E501
            flux_obj.setReaction(rid)
            flux_obj.setCoefficient(cobra_reaction.objective_coefficient)

    # write groups
    if len(cobra_model.groups) > 0:
        doc.enablePackage(
            "http://www.sbml.org/sbml/level3/version1/groups/version1", "groups", True
        )
        doc.setPackageRequired("groups", False)
        model_group: libsbml.GroupsModelPlugin = model.getPlugin("groups")  # noqa: E501
        for cobra_group in cobra_model.groups:
            group: libsbml.Group = model_group.createGroup()
            if f_replace and F_GROUP_REV in f_replace:
                gid = f_replace[F_GROUP_REV](cobra_group.id)
            else:
                gid = cobra_group.id
            group.setId(gid)
            group.setName(cobra_group.name)
            group.setKind(cobra_group.kind)

            _sbase_notes_dict(group, cobra_group.notes)
            _sbase_annotations(group, cobra_group.annotation)

            for cobra_member in cobra_group.members:
                member: libsbml.Member = group.createMember()
                mid = cobra_member.id
                m_type = str(type(cobra_member))

                # id replacements
                if "Reaction" in m_type:
                    if f_replace and F_REACTION_REV in f_replace:
                        mid = f_replace[F_REACTION_REV](mid)
                if "Metabolite" in m_type or "Protein" in m_type:
                    if f_replace and F_SPECIE_REV in f_replace:
                        mid = f_replace[F_SPECIE_REV](mid)
                if "Gene" in m_type:
                    if f_replace and F_GENE_REV in f_replace:
                        mid = f_replace[F_GENE_REV](mid)

                member.setIdRef(mid)
                if cobra_member.name and len(cobra_member.name) > 0:
                    member.setName(cobra_member.name)

    if isinstance(filename, str):
        # write to path
        libsbml.writeSBMLToFile(doc, filename)

    elif hasattr(filename, "write"):
        # write to file handle
        sbml_str = libsbml.writeSBMLToString(doc)
        filename.write(sbml_str)

# Copyright 2021 <NAME> # 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. Read and write SBML directly with `geckopy.Model`s.  Proteins are Species, members of Group Protein, with:     - initialAmount: Concentration.     - boundary_condition: True  The MW are initialAssingments. # re.Pattern does not exist in py36 so this type hint cannot be added now Create `geckopy.Model` from SBMLDocument.      Parameters     ----------     filename: str     number: data type of stoichiometry: {float, int}         In which data type should the stoichiometry be parsed.     f_replace : dict of replacement functions for id replacement     set_missing_bounds : flag to set missing bounds     hardcoded_rev_reactions: bool         if reversible reaction to account for proteins being consumed on both         directions are written explicitly for in the SBML      Returns     -------     cobra.core.Model # SBML model # fbc-v1 (legacy) # type: libsbml.FbcSBMLDocumentPlugin # Model # meta information # History # type: libsbml.ModelHistory # type: libsbml.ModelCreator # type:libsbml.SBasePlugin # notes and annotations # Compartments # FIXME: update with new compartments # noqa: E501 type: libsbml.Compartment # Species # proteins that rely on the naming convention "prot_UNIPROT_ID" will be # catched here. Those who are annotated by groups membership will be parsed # after the groups are processed. # type: libsbml.Species # type: libsbml.FbcSpeciesPlugin # handle nan, na, NA, ... # Detect boundary metabolites # Add exchange reactions for boundary metabolites # species is reactant # Genes # noqa: E501 type: libsbml.GeneProduct # noqa: E501 type: libsbml.Reaction # fallback to notes information # Interaction of the above replacements can lead to multiple # ;, which results in empty gids # create missing genes # GPR rules Recursively convert gpr association to a gpr string.          Defined as inline functions to not pass the replacement dict around. # Reactions # type: libsbml.Reaction # proteins are parsed based on Species, prot exchanges are ignored # set bounds # type: libsbml.FbcReactionPlugin # bounds in fbc # type: libsbml.Parameter # type: libsbml.Parameter # some legacy models encode bounds in kinetic laws # type: libsbml.KineticLaw # noqa: E501 type: libsbml.LocalParameter # noqa: E501 type: libsbml.LocalParameter # add reaction # parse equation # noqa: E501 type: libsbml.SpeciesReference # noqa: E501 type: libsbml.SpeciesReference # convert to metabolite objects # GPR # noqa: E501 type: libsbml.GeneProductAssociation # noqa: E501 type: libsbml.FbcAssociation # fallback to notes information # remove outside parenthesis, if any # Objective # noqa: E501 type: libsbml.ListOfObjectives # type: libsbml.Objective # noqa: E501 type: libsbml.FluxObjective # some legacy models encode objective coefficients in kinetic laws # type: libsbml.Reaction # type: libsbml.KineticLaw # type: libsbml.LocalParameter # parse groups # type: libsbml.GroupsModelPlugin # calculate hashmaps to lookup objects in O(1) # type: libsbml.SBase # create groups # type: libsbml.Group # type: libsbml.Member # id replacements # parse deprecated subsystems on reactions # now add everything under group Proteins to model.proteins if it was not # already added based on naming conventions # reverse the (negative) stoichiometry coefficient to kcat # we expect that Proteins that are identified only by their # group are respecting the specification and do form part # of reactions # TODO: provide options to tune this behvior # general hint for missing flux bounds Write cobra model to filename.      Enzyme constraint changes: proteins are written as metabolites with     initialAmount.      The created model is SBML level 3 version 1 (L1V3) with     fbc package v2 (fbc-v2).     If the given filename ends with the suffix ".gz" (for example,     "myfile.xml.gz"), libSBML assumes the caller wants the file to be     written compressed in gzip format. Similarly, if the given filename     ends with ".zip" or ".bz2", libSBML assumes the caller wants the     file to be compressed in zip or bzip2 format (respectively). Files     whose names lack these suffixes will be written uncompressed. Special     considerations for the zip format: If the given filename ends with     ".zip", the file placed in the zip archive will have the suffix     ".xml" or ".sbml".  For example, the file in the zip archive will     be named "test.xml" if the given filename is "test.xml.zip" or     "test.zip". Similarly, the filename in the archive will be     "test.sbml" if the given filename is "test.sbml.zip".      Warnings     --------     If :meth:`geckopy.Model.constrain_pool()` was called on the model, the     protein pool (exchange and metabolite) will be ignored and NOT written     to the SBML file. A `UserWarning` appears in this case.      Parameters     ----------     cobra_model : geckopy.Model         Model instance which is written to SBML     filename : string         path to which the model is written     f_replace: dict of replacement functions for id replacement     group_untyped_proteins: bool         if True (default) the proteins whose id is not compliant with Uniprot         and which are not part of the :code:`Protein` group will be added to it.         This is an inplace operation! # SBML L3V1 # fbc-v2 # for parsing annotation corresponding to the model # for parsing notes corresponding to the model # Meta information (ModelHistory) related to SBMLDocument # noqa: E501 # TODO: Will be implemented as part of #  https://github.com/opencobra/cobrapy/issues/810 # _check(model.setModelHistory(history), 'set model history') # Units # type:libsbml.UnitDefinition # type: libsbml.Unit # minimum and maximum value from model # Compartments # Species # Genes # Objective # Reactions # stoichiometry # noqa: E501 type: libsbml.SpeciesReference # noqa: E501 type: libsbml.SpeciesReference # bounds # GPR # replace ids in string # noqa: E501 # uses ids to identify GeneProducts (True), # does not create GeneProducts (False) # objective coefficients # noqa: E501 # write groups # noqa: E501 # id replacements # write to path # write to file handle