class StrainDesignMethod(object):
    def __init__(self, *args, **kwargs):
        super(StrainDesignMethod, self).__init__(*args, **kwargs)

    def __call__(self, *args, **kwargs):
        self.run(*args, **kwargs)

    def run(self, *args, **kwargs):
        raise NotImplementedError class StrainDesignMethodResult(Result):
    __method_name__ = None

    _aggreate_functions_ = {
        "method": lambda series: tuple(sum(series.values.tolist(), []))
    }

    def __init__(self, designs, *args, **kwargs):
        super(StrainDesignMethodResult, self).__init__(*args, **kwargs)
        self._designs = designs

    def __len__(self):
        """
        Return the number of designs found.
        """
        return len(self._designs)

    def __iter__(self):
        """
        Returns an iterator that yields StrainDesign objects.
        """
        return iter(self._designs)

    @property
    def data_frame(self):
        df = DataFrame(columns=['targets'])
        for i, design in enumerate(self._designs):
            df.loc[i] = [design.targets]
        return df

    def display_on_map(self, index, map_name):
        raise NotImplementedError

    def __add__(self, other):
        df = DataFrame([design.targets for design in self._designs], columns=["targets"])
        df['method'] = self.__method_name__

        i = len(df)
        for j, design in enumerate(other):
            df.loc[i + j] = [design, [self.__method_name__]]

        df = df.groupby(["design"]).aggregate(self._aggreate_functions_)

        return StrainDesignMethodEnsemble([StrainDesign(targets) for targets in df.targets], df.method.tolist())

    def _repr_html_(self):
        return self.data_frame._repr_html_()

    def plot(self, plotter, grid=None, width=None, height=None, title=None, *args, **kwargs):
        if title is None:
            title = "Target frequency plot for %s result" % self.__method_name__

        counts = DataFrame(columns=['count'])
        for design in self._designs:
            for target in design:
                if target.id not in counts.index:
                    counts.loc[target.id, 'count'] = 0
                counts.loc[target.id, 'count'] += 1
            counts['frequency'] = counts['count'].apply(lambda c: c / len(self._designs))

        plotter.frequency(counts, title=title, width=width, height=height, **kwargs) class StrainDesign(object):
    """
    A StrainDesign is a collection of targets in a COBRA model. The targets, identified by a StrainDesignMethod,
    map elements in the model that need to be modified to achieve the objective of the design method.
    """
    def __init__(self, targets):
        self.targets = DictList(targets)

    def __str__(self):
        return ", ".join(str(t) for t in self.targets)

    def __repr__(self):
        return "<StrainDesign [" + ";".join(repr(t) for t in self.targets) + "]>"

    def __iter__(self):
        return iter(self.targets)

    def __len__(self):
        return len(self.targets)

    def __contains__(self, item):
        if isinstance(item, Target):
            if item in self.targets:
                return item == self.targets.get_by_id(item.id)
            else:
                return False
        elif isinstance(item, str):
            return item in self.targets
        else:
            return False

    def __eq__(self, other):
        if isinstance(other, StrainDesign):
            if len(self) != len(other):
                return False
            else:
                return all(t in other for t in self.targets) and all(t in self for t in other.targets)
        else:
            return False

    def apply(self, model):
        for target in self.targets:
            target.apply(model)

    def __add__(self, other):
        if not isinstance(other, StrainDesign):
            raise AssertionError("Only instances of StrainDesign can be added together")

        targets = {}
        for target in self.targets:
            if target.id not in targets:
                targets[target.id] = set()
            targets[target.id].add(target)

        for target in other.targets:
            if target.id not in targets:
                targets[target.id] = set()
            targets[target.id].add(target)

        targets = [next(iter(t)) if len(t) == 1 else EnsembleTarget(id, t) for id, t in targets.items()]

        return StrainDesign(targets)

    def __iadd__(self, other):
        if not isinstance(other, StrainDesign):
            raise AssertionError("Only instances of StrainDesign can be added together")

        targets = {}
        for target in self.targets:
            if target.id not in targets:
                targets[target.id] = set()
            targets[target.id].add(target)

        for target in other.targets:
            if target.id not in targets:
                targets[target.id] = set()
            targets[target.id].add(target)

        targets = [next(iter(t)) if len(t) == 1 else EnsembleTarget(id, t) for id, t in targets.items()]

        self.targets = DictList(targets)

        return self

    def _repr_html_(self):
        return " ".join(t._repr_html_() for t in self.targets)

    def to_gnomic(self):
        return Genotype([target.to_gnomic() for target in self.targets]) class ReactionKnockoutTarget(ReactionModulationTarget):
    """
    Reaction Knockout Target. Knockout a reaction present in a COBRA model.
    """

    def __init__(self, id):
        super(ReactionKnockoutTarget, self).__init__(id, 0, None)

    def apply(self, model):
        target = self.get_model_target(model)
        target.knock_out()

    def __gt__(self, other):
        if self.id == other.id:
            if isinstance(other, ReactionModulationTarget):
                if other._value == 0:
                    return True
                else:
                    raise IncompatibleTargets(self, other)
            elif isinstance(other, ReactionCofactorSwapTarget):
                raise IncompatibleTargets(self, other)
            elif isinstance(other, EnsembleTarget):
                return not other > self
            else:
                raise IncompatibleTargets(self, other)
        else:
            return self.id > other.id

    def __eq__(self, other):
        if isinstance(other, ReactionKnockoutTarget):
            return self.id == other.id
        elif isinstance(other, ReactionModulationTarget):
            return self.id == other.id and other._value == 0
        else:
            return False

    def __repr__(self):
        return "<ReactionKnockout %s>" % self.id

    def __hash__(self):
        return hash(str(self)) class GeneKnockoutTarget(GeneModulationTarget):
    """
    Gene Knockout Target. Knockout a gene present in a COBRA model.
    """

    def __init__(self, id, *args, **kwargs):
        super(GeneKnockoutTarget, self).__init__(id, 0, None, *args, **kwargs)

    def apply(self, model):
        target = self.get_model_target(model)
        target.knock_out()

    def __gt__(self, other):
        if self.id == other.id:
            if isinstance(other, GeneModulationTarget) and not isinstance(other, GeneKnockoutTarget):
                return self.fold_change > other.fold_change
            elif isinstance(other, GeneKnockoutTarget):
                return False
            elif isinstance(other, EnsembleTarget):
                return not other > self
            else:
                raise IncompatibleTargets(self, other)
        else:
            return self.id > other.id

    def __eq__(self, other):
        if isinstance(other, GeneKnockoutTarget):
            return self.id == other.id
        elif isinstance(other, GeneModulationTarget):
            return self.id == other.id and other._value == 0
        else:
            return False

    def __repr__(self):
        return "<GeneKnockout %s>" % self.id

    def __hash__(self):
        return hash(str(self))

    def to_gnomic(self):
        accession = Target.to_gnomic(self)
        feature = Feature(name=self.id, accession=accession)
        return Change(before=feature) class ReactionCofactorSwapTarget(Target):
    """
    Swap cofactors of a given reaction.
    """

    def __init__(self, id, swap_pairs, *args, **kwargs):
        super(ReactionCofactorSwapTarget, self).__init__(id, *args, **kwargs)
        self.swap_pairs = swap_pairs

    def apply(self, model):
        reaction = model.reactions.get_by_id(self.id)
        swap_cofactors(reaction, model, self.swap_pairs)

    @property
    def swap_str(self):
        return "+".join(str(s) for s in self.swap_pairs[0]) + "<->" + "+".join(str(s) for s in self.swap_pairs[1])

    def to_gnomic(self):
        accession = Target.to_gnomic(self)
        pairs = [(m.id for m in pair) for pair in self.swap_pairs]
        wt_feature = Feature(name=self.id, accession=accession, type='reaction',
                             variant=["cofactors=%s" % ",".join(pairs[0])])

        feature = Feature(name=self.id, accession=accession, type='reaction',
                          variant=["cofactors=%s" % ",".join(pairs[1])])
        return Change(before=wt_feature, after=feature)

    def __repr__(self):
        return "<ReactionCofactorSwap %s swap=%s>" % (self.id, self.swap_str)

    def __str__(self):
        return "%s[%s]" % (self.id, self.swap_str)

    def __gt__(self, other):
        if self.id == other.id:
            if isinstance(other, ReactionKnockoutTarget):
                raise IncompatibleTargets(self, other)
            if isinstance(other, ReactionModulationTarget):
                return False
            if isinstance(other, ReactionKnockinTarget):
                return True
            elif isinstance(other, EnsembleTarget):
                return not other > self
            else:
                raise IncompatibleTargets(self, other)
        else:
            return self.id > other.id

    def __eq__(self, other):
        if isinstance(other, ReactionCofactorSwapTarget):
            return self.id == other.id and self.swap_pairs == other.swap_pairs
        else:
            return False

    def _repr_html_(self):
        return self.id + "|" + self.swap_str.replace("<->", "&rlarr;")

    def __hash__(self):
        return hash(str(self)) def swap_cofactors(reaction, model, swap_pairs, inplace=True):
    """
    Swaps the cofactors of a reaction. For speed, it can be done inplace which just changes the coefficients.
    If not done inplace, it will create a new Reaction, add it to the model, and knockout the original reaction.

    Parameters
    ----------
    reaction: cobra.Reaction
        The reaction to swap.
    model: cameo.cobra.Model
        A constraint-based model.
    swap_pairs: tuple
        A tuple of (cofactors, equivalent_cofactors)
    inplace: bool
        If replace is done inplace, it changes the coefficients in the matrix. Otherwise, it creates a new reaction
        with the other cofactors and adds it to the model.

    Returns
    -------
        Reaction
            A reaction with swapped cofactors (the same if inplace).
    """
    if all(reaction.metabolites.get(met, False) for met in swap_pairs[0]):
        new_coefficients = {met: -reaction.metabolites[met] for met in swap_pairs[0]}
        new_coefficients.update({new_met: reaction.metabolites[met] for met, new_met in zip(*swap_pairs)})
    elif all(reaction.metabolites.get(met, False) for met in swap_pairs[1]):
        new_coefficients = {met: -reaction.metabolites[met] for met in swap_pairs[1]}
        new_coefficients.update({new_met: reaction.metabolites[met] for new_met, met in zip(*swap_pairs)})
    else:
        raise ValueError("%s: Invalid swap pairs %s (%s)" % (reaction.id, str(swap_pairs), reaction.reaction))

    def _inplace(rxn, stoichiometry):
        rxn.add_metabolites(stoichiometry, combine=True)

    def _replace(rxn, stoichiometry):
        new_reaction = Reaction(id="%s_swap" % rxn.id, name=rxn.name,
                                lower_bound=rxn.lower_bound, upper_bound=rxn.upper_bound)
        new_reaction.stoichiometry = rxn
        new_reaction.add_metabolites(stoichiometry)
        return new_reaction

    if inplace:
        _inplace(reaction, new_coefficients)
        return reaction
    else:
        new_reaction = _replace(reaction, new_coefficients)
        model.add_reactions([new_reaction])
        reaction.knock_out()
        return new_reaction def phenotypic_phase_plane(model, variables, objective=None, source=None, points=20, view=None):
    """Phenotypic phase plane analysis [1].

    Implements a phenotypic phase plan analysis with interpretation same as
    that presented in [1] but calculated by optimizing the model for all
    steps of the indicated variables (instead of using shadow prices).

    Parameters
    ----------
    model: cobra.Model
    variables: str or reaction or iterable
        A reaction ID, reaction, or list of reactions to be varied.
    objective: str or reaction or optlang.Objective or Metabolite, optional
        An objective, a reaction's flux, or a metabolite's production to be minimized/maximized
        (defaults to the current model objective).
    source: Reaction or reaction identifier
        The reaction to use as the source when calculating mass and carbon yield. Set to the medium reaction with the
        highest input carbon flux if left None.
    points: int or iterable
        Number of points to be interspersed between the variable bounds.
        A list of same same dimensions as `variables` can be used to specify
        variable specific numbers of points.
    view: SequentialView or MultiprocessingView or ipython.cluster.DirectView
        A parallelization view.

    Returns
    -------
    PhenotypicPhasePlaneResult
        The phenotypic phase plane with flux, mol carbon input / mol carbon
        output, gram input / gram output

    References
    ----------
    [1] Edwards, J. S., Ramakrishna, R. and Palsson, B. O. (2002). Characterizing the metabolic phenotype: a phenotype
        phase plane analysis. Biotechnology and Bioengineering, 77(1), 27–36. doi:10.1002/bit.10047
    """

    if isinstance(variables, str):
        variables = [variables]
    elif isinstance(variables, Reaction):
        variables = [variables]
    variable_ids = [var if isinstance(var, str) else var.id for var in variables]

    if view is None:
        view = config.default_view
    with model:
        if objective is not None:
            if isinstance(objective, Metabolite):
                try:
                    objective = model.reactions.get_by_id("DM_%s" % objective.id)
                except KeyError:
                    objective = model.add_boundary(objective, type='demand')
            # try:
            #     objective = model.reaction_for(objective, time_machine=tm)
            # except KeyError:
            #     pass

            model.objective = objective

        if source is not None:
            source_reaction = get_reaction_for(model, source)
        else:
            source_reaction = get_c_source_reaction(model)

        variable_reactions = model.reactions.get_by_any(variables)
        variables_min_max = flux_variability_analysis(model, reactions=variable_reactions, view=SequentialView())
        grid = [numpy.linspace(lower_bound, upper_bound, points, endpoint=True) for
                reaction_id, lower_bound, upper_bound in
                variables_min_max.data_frame.loc[variable_ids].itertuples()]
        grid_generator = itertools.product(*grid)
        chunks_of_points = partition(list(grid_generator), len(view))
        evaluator = _PhenotypicPhasePlaneChunkEvaluator(model, variable_reactions, source_reaction)
        chunk_results = view.map(evaluator, chunks_of_points)
        envelope = reduce(list.__add__, chunk_results)

    nice_variable_ids = [_nice_id(reaction) for reaction in variable_reactions]
    variable_reactions_ids = [reaction.id for reaction in variable_reactions]
    phase_plane = pandas.DataFrame(
        envelope, columns=(variable_reactions_ids + [
            'objective_lower_bound',
            'objective_upper_bound',
            'c_yield_lower_bound',
            'c_yield_upper_bound',
            'mass_yield_lower_bound',
            'mass_yield_upper_bound'
        ])
    )

    if objective is None:
        objective = model.objective
    nice_objective_id = _nice_id(objective)

    objective = objective.id if isinstance(objective, Reaction) else str(objective)
    return PhenotypicPhasePlaneResult(phase_plane, variable_reactions_ids, objective,
                                      nice_variable_ids=nice_variable_ids,
                                      source_reaction=_nice_id(source_reaction),
                                      nice_objective_id=nice_objective_id) def fba(model, objective=None, reactions=None, *args, **kwargs):
    """Flux Balance Analysis.

    Parameters
    ----------
    model: cobra.Model
    objective: a valid objective - see SolverBaseModel.objective (optional)

    Returns
    -------
    FluxDistributionResult
        Contains the result of the linear solver.

    """
    with model:
        if objective is not None:
            model.objective = objective
        solution = model.optimize(raise_error=True)
        if reactions is not None:
            result = FluxDistributionResult({r: solution[r] for r in reactions}, solution.objective_value)
        else:
            result = FluxDistributionResult.from_solution(solution)
        return result class ProductionStrainArchive(BestSolutionArchive):
    def __call__(self, random, population, archive, args):
        self.archive = archive
        max_archive_size = args.get('max_archive_size', 100)
        [self.add(i.candidate, i.fitness, i.birthdate, True, max_archive_size) for i in population
         if round(i.fitness, ndecimals) > 0]
        return self.archive class biomass_product_coupled_min_yield(biomass_product_coupled_yield):
    """
    Biomass-Product Coupled Minimum Yield: (v[biomass] * min(v[product])) / v[substrate] [1]

    Parameters
    ----------
    biomass: str or Reaction
        biomass reaction identifier
    product: str or Reaction
        product reaction identifier
    substrate: str or Reaction
        substrate reaction identifier

    Returns
    -------
    float
        fitness value

    """

    def __call__(self, model, solution, targets):
        biomass_flux = round(solution.fluxes[self.biomass], config.ndecimals)
        fva_res = flux_variability_analysis(model, reactions=[self.product], fraction_of_optimum=1)
        min_product_flux = round(fva_res["lower_bound"][self.product], config.ndecimals)
        if self.carbon_yield:
            product = model.reactions.get_by_id(self.product)
            if product.boundary:
                product_flux = min_product_flux * n_carbon(product)
            else:
                product_flux = min_product_flux * n_carbon(product) / 2
            substrate_flux = 0
            for substrate_id in self.substrates:
                substrate = model.reactions.get_by_id(substrate_id)
                if substrate.boundary:
                    substrate_flux += abs(solution.fluxes[substrate_id]) * n_carbon(substrate)
                else:
                    substrate_flux += abs(solution.fluxes[substrate_id]) * n_carbon(substrate) / 2
            substrate_flux = round(substrate_flux, config.ndecimals)
        else:
            product_flux = min_product_flux
            substrate_flux = sum(round(abs(solution.fluxes[s]), config.ndecimals) for s in self.substrates)

        try:
            return (biomass_flux * product_flux) / substrate_flux
        except ZeroDivisionError:
            return 0.0

    def _repr_latex_(self):
        if self.carbon_yield:
            substrates = " + ".join("C(%s)" % s for s in self.substrates)
            return "$$bpcy = \\frac{(%s * Cmin(%s))}{%s}$$" % (self.biomass.replace("_", "\\_"),
                                                               self.product.replace("_", "\\_"),
                                                               substrates)
        else:
            return "$$bpcy = \\frac{(%s * min(%s))}{%s}$$" % (self.biomass.replace("_", "\\_"),
                                                              self.product.replace("_", "\\_"),
                                                              " + ".join(
                                                                  s.replace("_", "\\_") for s in self.substrates))

    @property
    def name(self):
        if len(self.substrates) == 1:
            substrate = self.substrates[0]
        else:
            substrate = "(" + " + ".join(self.substrates) + ")"
        return "bpcy = (%s * min(%s)) / %s" % (self.biomass, self.product, substrate) class biomass_product_coupled_yield(YieldFunction):
    """
    Biomass-Product Coupled Yield: (v[biomass] * v[product]) / v[substrate] [1]

    Parameters
    ----------
    biomass: str or Reaction
        biomass reaction identifier
    product: str or Reaction
        product reaction identifier
    substrate: str or Reaction
        substrate reaction identifier

    Returns
    -------
    float
        fitness value

    References
    ----------
    [1] Patil, K. R., Rocha, I., Förster, J., & Nielsen, J. (2005). "Evolutionary programming as a
    platform for in silico metabolic engineering". BMC Bioinformatics, 6, 308.
    doi:10.1186/1471-2105-6-308
    """

    def __init__(self, biomass, product, substrate, carbon_yield=False, *args, **kwargs):
        super(biomass_product_coupled_yield, self).__init__(product, substrate,
                                                            carbon_yield=carbon_yield, *args, **kwargs)
        if isinstance(biomass, Reaction):
            biomass = biomass.id
        self.biomass = biomass

    def __call__(self, model, solution, targets):
        biomass_flux = round(solution.fluxes[self.biomass], config.ndecimals)
        if self.carbon_yield:
            product = model.reaction.get_by_id(self.product)
            if product.boundary:
                product_flux = round(solution.fluxes[self.product], config.ndecimals) * n_carbon(product)
            else:
                product_flux = round(solution.fluxes[self.product], config.ndecimals) * n_carbon(product) / 2
            substrate_flux = 0
            for substrate_id in self.substrates:
                substrate = model.reactions.get_by_id(substrate_id)
                if substrate.boundary:
                    substrate_flux += abs(solution.fluxes[substrate_id]) * n_carbon(substrate)
                else:
                    substrate_flux += abs(solution.fluxes[substrate_id]) * n_carbon(substrate) / 2
            substrate_flux = round(substrate_flux, config.ndecimals)
        else:
            product_flux = round(solution.fluxes[self.product], config.ndecimals)
            substrate_flux = round(sum(abs(solution.fluxes[s]) for s in self.substrates), config.ndecimals)
        if substrate_flux > config.non_zero_flux_threshold:
            return (biomass_flux * product_flux) / substrate_flux
        else:
            return 0.

    def _repr_latex_(self):
        if self.carbon_yield:
            substrates = " + ".join("C(%s)" % s for s in self.substrates)
            return "$$bpcy = \\frac{(%s * C(%s))}{%s}$$" % (self.biomass.replace("_", "\\_"),
                                                            self.product.replace("_", "\\_"),
                                                            substrates)

        else:
            return "$$bpcy = \\frac{(%s * %s)}{%s}$$" % (self.biomass.replace("_", "\\_"),
                                                         self.product.replace("_", "\\_"),
                                                         " + ".join(s.replace("_", "\\_") for s in self.substrates))

    @property
    def name(self):
        if len(self.substrates) == 1:
            substrate = self.substrates[0]
        else:
            substrate = "(" + " + ".join(self.substrates) + ")"
        return "bpcy = (%s * %s) / %s" % (self.biomass, self.product, substrate)

    @property
    def reactions(self):
        return [self.biomass, self.product] + self.substrates class GeneKnockoutOptimization(KnockoutOptimization):
    """
    Knockout optimization using genes.

    Attributes
    ----------
    model : cobra.Model
        A constraint-based model.
    heuristic_method : inspyred.ec.EvolutionaryComputation
        An evolutionary algorithm.
    objective_function : objective function or list(objective function)
        The objectives for the algorithm to maximize.
    seed : int
        A seed for random. It is auto-generated if None is given.
    termination : inspyred.ec.terminators
        A termination criteria for the algorithm. The default is inspyred.ec.terminators.evaluation_termination.
    simulation_method: flux_analysis.simulation
        The method used to simulate the model.
    wt_reference: dict
        A reference initial state for the optimization. It is required for flux_analysis.simulation.lmoma and
        flux_analysis.simulation.room. If not given, it will be computed using flux_analysis.simulation.pfba
    genes: list
        A list of valid genes to knockout. If None, then all genes in the model will be knockout candidates except the
         ones defined in essential_genes
    essential_genes: list
        A list of genes that cannot be knocked out. If None, then all essential genes will be removed from the valid
        genes set.
    use_nullspace_simplification: Boolean (default True)
        Use a basis for the nullspace dead end reactions. Gene present only in dead end reactions will be ignored.

    Methods
    -------
    run(view=config.default_view, maximize=True, **kwargs)


    See Also
    --------
    *inspyred.ec
    *cameo.config.default_view

    Examples
    --------
    >>> from cameo import models
    >>> model = models.bigg.iJO1366
    >>> from cameo.strain_design.heuristic.evolutionary.objective_functions import biomass_product_coupled_yield
    >>> bpcy = biomass_product_coupled_yield(model.reactions.Ec_biomass_iJO1366_core_53p95,
    >>>                                      model.reactions.EX_succ_e),
    >>>                                      model.reactions.EX_glc__D_e)
    >>> knockout_optimization = GeneKnockoutOptimization(model=model, objective_function=bpcy)
    >>> knockout_optimization.run(max_evaluations=50000)

    """

    def __init__(self, genes=None, essential_genes=None, use_nullspace_simplification=True, *args, **kwargs):
        super(GeneKnockoutOptimization, self).__init__(*args, **kwargs)
        if genes is None:
            self.genes = set([g.id for g in self.model.genes])
        else:
            self.genes = genes
        if essential_genes is None:
            self.essential_genes = {g.id for g in find_essential_genes(self.model, processes=1)}
        else:
            self.essential_genes = set([g.id for g in find_essential_genes(self.model, processes=1)] + essential_genes)

        # TODO: use genes from groups
        if use_nullspace_simplification:
            ns = nullspace(create_stoichiometric_array(self.model))
            dead_end_reactions = find_blocked_reactions_nullspace(self.model, ns=ns)
            dead_end_genes = {g.id for g in self.model.genes if all(r in dead_end_reactions for r in g.reactions)}
            exclude_genes = self.essential_genes.union(dead_end_genes)
            genes = [g for g in self.model.genes if g.id not in exclude_genes]
            self.representation = [g.id for g in genes]
        else:
            self.representation = list(self.genes.difference(self.essential_genes))

        self._target_type = GENE_KNOCKOUT_TYPE
        self._decoder = decoders.GeneSetDecoder(self.representation, self.model)
        self._evaluator = evaluators.KnockoutEvaluator(model=self.model,
                                                       decoder=self._decoder,
                                                       objective_function=self.objective_function,
                                                       simulation_method=self._simulation_method,
                                                       simulation_kwargs=self._simulation_kwargs) class ReactionKnockoutOptimization(KnockoutOptimization):
    """
    Knockout optimization using reactions.

    Attributes
    ----------
    model : cobra.Model
        A constraint-based model.
    heuristic_method : inspyred.ec.EvolutionaryComputation
        An evolutionary algorithm.
    objective_function : objective function or list(objective function)
        The objectives for the algorithm to maximize.
    seed : int
        A seed for random. It is auto-generated if None is given.
    termination : inspyred.ec.terminators
        A termination criteria for the algorithm. The default is inspyred.ec.terminators.evaluation_termination.
    simulation_method: flux_analysis.simulation
        The method used to simulate the model.
    wt_reference: dict
        A reference initial state for the optimization. It is required for flux_analysis.simulation.lmoma and
        flux_analysis.simulation.room. If not given, it will be computed using flux_analysis.simulation.pfba
    reactions: list
        A list of valid reactions to knockout. If None, then all reactions in the model will be knockout candidates
        except the ones defined in essential_reactions
    essential_reactions: list
        A list of reactions that cannot be knocked out. If None, then all essential reactions will be removed from
        the valid reactions set.
    use_nullspace_simplification: Boolean (default True)
        Use a basis for the nullspace to find groups of reactions whose fluxes are multiples of each other and dead
        end reactions. From each of these groups only 1 reaction will be included as a possible knockout.

    Methods
    -------
    run(view=config.default_view, maximize=True, **kwargs)


    See Also
    --------
    *inspyred.ec
    *cameo.config.default_view

    Examples
    --------
    >>> from cameo import models
    >>> model = models.bigg.iJO1366
    >>> from cameo.strain_design.heuristic.evolutionary.objective_functions import biomass_product_coupled_yield
    >>> bpcy = biomass_product_coupled_yield(model.reactions.Ec_biomass_iJO1366_core_53p95,
    >>>                                      model.reactions.EX_succ_e),
    >>>                                      model.reactions.EX_glc__D_e)
    >>> knockout_optimization = ReactionKnockoutOptimization(model=model, objective_function=bpcy,
    >>>                                                      essential_reactions=["ATPM"])
    >>> knockout_optimization.run(max_evaluations=50000)
    """

    def __init__(self, reactions=None, essential_reactions=None, use_nullspace_simplification=True, *args, **kwargs):
        super(ReactionKnockoutOptimization, self).__init__(*args, **kwargs)
        if reactions is None:
            self.reactions = set([r.id for r in self.model.reactions])
        else:
            self.reactions = reactions
        logger.debug("Computing essential reactions...")
        if essential_reactions is None:
            self.essential_reactions = set(r.id for r in find_essential_reactions(self.model, processes=1))
        else:
            self.essential_reactions = set([r.id for r in find_essential_reactions(self.model, processes=1)])
            self.essential_reactions.update(essential_reactions)

        if use_nullspace_simplification:
            ns = nullspace(create_stoichiometric_array(self.model))
            dead_ends = set(find_blocked_reactions_nullspace(self.model, ns=ns))
            exchanges = set(self.model.boundary)
            reactions = [
                r for r in self.model.reactions
                if (r not in exchanges) and (
                    r not in dead_ends) and (
                        r.id not in self.essential_reactions)
            ]

            groups = find_coupled_reactions_nullspace(self.model, ns=ns)
            groups_keys = [set(group) for group in groups if any(r.id in reactions for r in group)]
            reduced_set = reduce_reaction_set(reactions, groups_keys)
            to_keep = [r.id for r in reduced_set]

        else:
            groups = None
            to_keep = set(r.id for r in self.model.reactions)
            to_keep.difference_update(r.id for r in self.model.boundary)
            to_keep.difference_update(self.essential_reactions)
            to_keep = list(to_keep)

        self.representation = to_keep
        self._target_type = REACTION_KNOCKOUT_TYPE
        self._decoder = decoders.ReactionSetDecoder(self.representation, self.model, groups=groups)
        self._evaluator = evaluators.KnockoutEvaluator(model=self.model,
                                                       decoder=self._decoder,
                                                       objective_function=self.objective_function,
                                                       simulation_method=self._simulation_method,
                                                       simulation_kwargs=self._simulation_kwargs) class CofactorSwapOptimization(TargetOptimization):
    """
    Optimize co-factor swapping

    As suggested in [1]_, flux through a given reaction can sometimes be optimized by swapping complementary
    co-factor. This class implements a search for reactions when swapped improve the given objective. Briefly,
    the approach is to

    - find reactions that have all the targeted co-factor pairs e.g. (nad_c -> nadp_c, nadh_c -> nadph_c)

    - add reactions that have the co-factors swapped and then by a search algorithm switching one off in favor of the
      other

    The implementation here differs from that in [1]_ in that we use a general purpose search algorithm rather than
    formulating the search as a mixed integer linear programming problem.

    References
    ----------
    .. [1] King, Zachary A., and Adam M. Feist. "Optimizing Cofactor Specificity of Oxidoreductase Enzymes for the
      Generation of Microbial Production Strains - OptSwap." Industrial Biotechnology 9, no. 4 (August 1,
      2013): 236-46. - doi:10.1089/ind.2013.0005.

    Parameters
    ----------
    model : cobra.Model
       the model to operator on
    cofactor_id_swaps : tuple
       a tuple of length 2 that defines two lists of metabolite identifiers that should be interchanged during the
       swap optimization see e.g. `NADH_NADPH` which is also the default.
    candidate_reactions : list
       reactions to consider for co-factor swap - if not given then search for all reactions that include the given
       cofactors
    skip_reactions : list
       reactions to not consider for co-factor swap, defaults to the objective function if not provided
    args, kwargs : keyword arguments
       passed on to super-classes, see in particular `objective_function`, `heuristic_method`, `termination`,
       `simulation_method`, `wt_reference`, of `HeuristicOptimization` and `max_size` of `HeuristicOptimization.run`

    Examples
    --------
    >>> from cameo import models
    >>> from cameo.strain_design.heuristic.evolutionary.objective_functions import product_yield
    >>> model = models.bigg.iJO1366
    >>> model.objective = model.reactions.EX_thr__L_e
    >>> model.reactions.BIOMASS_Ec_iJO1366_core_53p95M.lower_bound = 0.1
    >>> py = product_yield(model.reactions.EX_thr__L_e, model.reactions.EX_glc__D_e)
    >>> swap_optimization = CofactorSwapOptimization(model=model, objective_function=py)
    >>> swap_optimization.run(max_evaluations=2000, max_size=2)
    """

    def __init__(self, cofactor_id_swaps=NADH_NADPH, candidate_reactions=None, skip_reactions=None, *args, **kwargs):
        super(CofactorSwapOptimization, self).__init__(*args, **kwargs)
        self._target_type = SWAP_TYPE
        swap_pairs = ([self.model.metabolites.get_by_id(m) for m in cofactor_id_swaps[0]],
                      [self.model.metabolites.get_by_id(m) for m in cofactor_id_swaps[1]])
        self.metadata['swap_pairs'] = swap_pairs
        self.representation = candidate_reactions or self.find_swappable_reactions(self.model, swap_pairs)
        if skip_reactions:
            self.representation -= skip_reactions
        self._decoder = decoders.ReactionSetDecoder(self.representation, self.model)

        self._evaluator = evaluators.SwapEvaluator(model=self.model,
                                                   decoder=self._decoder,
                                                   objective_function=self.objective_function,
                                                   simulation_method=self._simulation_method,
                                                   simulation_kwargs=self._simulation_kwargs,
                                                   swap_pair=swap_pairs)

    @staticmethod
    def find_swappable_reactions(model, swaps):
        """
        Get all reactions that can undergo co-factor swapping

        Find reactions that have one set of the cofactors targeted for swapping and are mass balanced and updates the
        `candidate_reactions` attribute

        Parameters
        ----------
        model: cobra.Model
            A model with reactions to search on.
        swaps: tuple
            Pair of cofactors to swap.

        """

        def swap_search(mets):
            has_pairs = all(mets.get(m, False) for m in swaps[0]) or all(mets.get(m, False) for m in swaps[1])

            contains_all = all(mets.get(m, False) for m in swaps[0]) and all(mets.get(m, False) for m in swaps[1])

            return has_pairs and not contains_all

        candidates = model.reactions.query(search_function=swap_search, attribute='metabolites')
        return [r.id for r in candidates] NADH_NADPH = (['nad_c', 'nadh_c'], ['nadp_c', 'nadph_c']) def process_reaction_knockout_solution(model, solution, simulation_method, simulation_kwargs,
                                       biomass, target, substrate, objective_function):
    """

    Parameters
    ----------

    model: cobra.Model
        A constraint-based model
    solution: tuple - (reactions, knockouts)
        The output of a decoder
    simulation_method: function
        See see cameo.flux_analysis.simulation
    simulation_kwargs: dict
        Keyword arguments to run the simulation method
    biomass: Reaction
        Cellular biomass reaction
    target: Reaction
        The strain design target
    substrate: Reaction
        The main carbon source uptake rate
    objective_function: cameo.strain_design.heuristic.evolutionary.objective_functions.ObjectiveFunction
        The objective function used for evaluation.

    Returns
    -------
    list
        A list with: reactions, size, fva_min, fva_max, target flux, biomass flux, yield, fitness

    """

    with model:
        reactions = [model.reactions.get_by_id(rid) for rid in solution]
        for reaction in reactions:
            reaction.knock_out()

        flux_dist = simulation_method(model, reactions=objective_function.reactions,
                                      objective=biomass, **simulation_kwargs)
        model.objective = biomass
        fva = flux_variability_analysis(model, fraction_of_optimum=0.99, reactions=[target])
        target_yield = flux_dist[target] / abs(flux_dist[substrate])
        return [solution, len(solution), fva.lower_bound(target),
                fva.upper_bound(target), flux_dist[target], flux_dist[biomass],
                target_yield, objective_function(model, flux_dist, reactions)] def process_gene_knockout_solution(model, solution, simulation_method, simulation_kwargs,
                                   biomass, target, substrate, objective_function):
    """

    Parameters
    ----------
    model: cobra.Model
        A constraint-based model
    solution: tuple
        The genes
    simulation_method: function
        See see cameo.flux_analysis.simulation
    simulation_kwargs: dict
        Keyword arguments to run the simulation method
    biomass: Reaction
        Cellular biomass reaction
    target: Reaction
        The strain design target
    substrate: Reaction
        The main carbon source uptake rate
    objective_function: cameo.strain_design.heuristic.evolutionary.objective_functions.ObjectiveFunction
        The objective function used for evaluation.

    Returns
    -------
    list
        A list with: reactions, genes, size, fva_min, fva_max, target flux, biomass flux, yield, fitness

    """

    with model:
        genes = [model.genes.get_by_id(gid) for gid in solution]
        reactions = find_gene_knockout_reactions(model, solution)
        for reaction in reactions:
            reaction.knock_out()

        reaction_ids = [r.id for r in reactions]
        flux_dist = simulation_method(model, reactions=objective_function.reactions,
                                      objective=biomass, **simulation_kwargs)
        model.objective = biomass

        fva = flux_variability_analysis(model, fraction_of_optimum=0.99, reactions=[target])
        target_yield = flux_dist[target] / abs(flux_dist[substrate])

        return [tuple(reaction_ids), solution, len(solution), fva.lower_bound(target), fva.upper_bound(target),
                flux_dist[target], flux_dist[biomass], target_yield, objective_function(model, flux_dist, genes)] def process_reaction_swap_solution(model, solution, simulation_method, simulation_kwargs, biomass,
                                   target, substrate, objective_function, swap_pairs):
    """

    Parameters
    ----------
    model: cobra.Model
        A constraint-based model
    solution: tuple - (reactions, knockouts)
        The output of a decoder
    simulation_method: function
        See see cameo.flux_analysis.simulation
    simulation_kwargs: dict
        Keyword arguments to run the simulation method
    biomass: Reaction
        Cellular biomass reaction
    target: Reaction
        The strain design target
    substrate: Reaction
        The main carbon source uptake rate
    objective_function: cameo.strain_design.heuristic.evolutionary.objective_functions.ObjectiveFunction
        The objective function used for evaluation.
    swap_pairs:
        The metabolites to swap

    Returns
    -------
    list
        A list with: reactions, size, fva_min, fva_max, target flux, biomass flux, yield, fitness,
        [fitness, [fitness]]

    """

    with model:
        reactions = [model.reactions.get_by_id(rid) for rid in solution]
        for reaction in reactions:
            swap_cofactors(reaction, model, swap_pairs)

        flux_dist = simulation_method(model, reactions=objective_function.reactions,
                                      objective=biomass, **simulation_kwargs)
        model.objective = biomass
        fva = flux_variability_analysis(model, fraction_of_optimum=0.99, reactions=[target])
        target_yield = flux_dist[target] / abs(flux_dist[substrate])
        return [solution, fva.lower_bound(target),
                fva.upper_bound(target), flux_dist[target], flux_dist[biomass],
                target_yield, objective_function(model, flux_dist, reactions)] class TimeMachine(object):
    """Travel back and forth in time."""

    def __init__(self):
        super(TimeMachine, self).__init__()
        self.history = OrderedDict()

    def __call__(self, do=None, undo=None, bookmark=None):
        output = do()
        current_time = time()
        if bookmark is None:
            entry_id = uuid1()
        else:
            entry_id = bookmark
        # make sure that entry is added to the end of history
        self.history.pop(entry_id, None)
        self.history[entry_id] = {'unix_epoch': current_time, 'undo': undo, 'redo': do}
        return entry_id, output

    def __str__(self):
        info = '\n'
        for item in self.history.items():
            info += self._history_item_to_str(item)
        return info

    def __enter__(self):
        return self

    def __exit__(self, type, value, traceback):
        self.reset()

    @staticmethod
    def _history_item_to_str(item):
        info = ''
        uuid, entry = item
        info += datetime.fromtimestamp(entry['unix_epoch']).strftime('%Y-%m-%d %H:%M:%S') + '\n'
        undo_entry = entry['undo']
        try:
            # partial  (if .keywords is None print {} instead)
            elements = undo_entry.func, undo_entry.args, undo_entry.keywords or {}
            info += 'undo: ' + ' '.join([str(elem) for elem in elements]) + '\n'
        except AttributeError:  # normal python function
            info += 'undo: ' + undo_entry.__name__ + '\n'

        redo_entry = entry['redo']
        try:
            elements = redo_entry.func, redo_entry.args, redo_entry.keywords or {}  # partial
            info += 'redo: ' + ' '.join([str(elem) for elem in elements]) + '\n'
        except AttributeError:
            info += 'redo: ' + redo_entry.__name__ + '\n'
        return info

    def undo(self, bookmark=None):
        if bookmark is None:
            try:
                (uuid, entry) = self.history.popitem()
                entry['undo']()
            except KeyError:  # history is empty
                pass
        elif bookmark in list(self.history.keys()):
            uuid = False
            while uuid is not bookmark:
                (uuid, entry) = self.history.popitem()
                entry['undo']()
        else:
            raise Exception(
                'Provided bookmark %s cannot be found in the time machine.')

    def redo(self):
        raise NotImplementedError

    def reset(self):
        if self.history:  # history is not empty
            self.undo(bookmark=list(self.history.keys())[0]) def get_reaction_for(model, value, add=True):
    """Get or create a reaction for a metabolite or a reaction.

    If value is a Metabolite or a Metabolite id, return any already existing demand or exchange reaction.
    If *add* is true, add a demand reaction if it does not already exist.

    Parameters
    ----------
    model : cobra.Model
        The model to for which to get / create a reaction
    value: str, Reaction or Metabolite
        A reaction identifier, a Reaction or a Metabolite for which an exchange reaction is to be created.
    add: bool
        Adds a demand reaction for a metabolite if a metabolite is found for *value*

    Returns
    -------
    Reaction

    Raises
    ------
    KeyError
        If *value* does not match any Reaction or Metabolite

    """
    try:
        reactions = model.reactions.get_by_any(value)
    except (ValueError, KeyError, TypeError):
        metabolite = model.metabolites.get_by_any(value)[0]
        reactions = model.reactions.query("^(EX|DM)_{}$".format(metabolite.id))
        if len(reactions) == 0:
            if add:
                reactions = [model.add_boundary(metabolite, type='demand')]
            else:
                raise KeyError('Invalid target %s' % value)
    return reactions[0]
import logging
import inspyred
import numpy
from IProgress.progressbar import ProgressBar
from IProgress.widgets import Bar, Percentage
from pandas import DataFrame
from cobra import Model
from cameo.core.strain_design import StrainDesignMethod, StrainDesignMethodResult, StrainDesign
from cameo.core.target import ReactionKnockoutTarget, GeneKnockoutTarget, ReactionCofactorSwapTarget
from cameo.core.manipulation import swap_cofactors
from cobra.exceptions import OptimizationError
from cameo.flux_analysis.analysis import phenotypic_phase_plane
from cameo.flux_analysis.simulation import fba
from cameo.strain_design.heuristic.evolutionary.archives import ProductionStrainArchive
from cameo.strain_design.heuristic.evolutionary.objective_functions import biomass_product_coupled_min_yield, \
    biomass_product_coupled_yield
from cameo.strain_design.heuristic.evolutionary.optimization import GeneKnockoutOptimization, \
    ReactionKnockoutOptimization, CofactorSwapOptimization, NADH_NADPH
from cameo.strain_design.heuristic.evolutionary.processing import process_reaction_knockout_solution, \
    process_gene_knockout_solution, process_reaction_swap_solution
from cameo.util import TimeMachine
from cameo.core.utils import get_reaction_for
# 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.


from __future__ import absolute_import, print_function




__all__ = ["OptGene"]


logger = logging.getLogger(__name__)


class OptGene(StrainDesignMethod):
    def __init__(self, model, evolutionary_algorithm=inspyred.ec.GA, manipulation_type="genes", essential_genes=None,
                 essential_reactions=None, plot=True, exclude_non_gene_reactions=True, *args, **kwargs):
        if not isinstance(model, Model):
            raise TypeError("Argument 'model' should be of type 'cobra.Model'.")

        super(OptGene, self).__init__(*args, **kwargs)

        if exclude_non_gene_reactions:
            essential_reactions = essential_reactions or []
            essential_reactions += [r for r in model.reactions if not r.genes]

        self._model = model
        self._algorithm = evolutionary_algorithm
        self._optimization_algorithm = None
        self._manipulation_type = None
        self._essential_genes = essential_genes
        self._essential_reactions = essential_reactions
        self._plot = plot
        self._manipulation_type = manipulation_type

    @property
    def manipulation_type(self):
        return self._manipulation_type

    @property
    def plot(self):
        return self._plot

    @plot.setter
    def plot(self, plot):
        self._plot = plot
        if self._optimization_algorithm is not None:
            self._optimization_algorithm.plot = plot

    @manipulation_type.setter
    def manipulation_type(self, manipulation_type):
        self._manipulation_type = manipulation_type

    def run(self, target=None, biomass=None, substrate=None, max_knockouts=5, variable_size=True,
            simulation_method=fba, growth_coupled=False, max_evaluations=20000, population_size=200,
            max_results=50, use_nullspace_simplification=True, seed=None, **kwargs):
        """
        Parameters
        ----------
        target : str, Metabolite or Reaction
            The design target
        biomass : str, Metabolite or Reaction
            The biomass definition in the model
        substrate : str, Metabolite or Reaction
            The main carbon source
        max_knockouts : int
            Max number of knockouts allowed
        variable_size : bool
            If true, all candidates have the same size. Otherwise the candidate size can be from 1 to max_knockouts.
        simulation_method : function
            Any method from cameo.flux_analysis.simulation or equivalent
        growth_coupled : bool
            If true will use the minimum flux rate to compute the fitness
        max_evaluations : int
            Number of evaluations before stop
        population_size : int
            Number of individuals in each generation
        max_results : int
            Max number of different designs to return if found.
        kwargs : dict
            Arguments for the simulation method.
        seed : int
            A seed for random.
        use_nullspace_simplification : Boolean (default True)
            Use a basis for the nullspace to find groups of reactions whose fluxes are multiples of each other and dead
            end reactions. From each of these groups only 1 reaction will be included as a possible knockout.


        Returns
        -------
        OptGeneResult
        """

        target = get_reaction_for(self._model, target)
        biomass = get_reaction_for(self._model, biomass)
        substrate = get_reaction_for(self._model, substrate)

        if growth_coupled:
            objective_function = biomass_product_coupled_min_yield(biomass, target, substrate)
        else:
            objective_function = biomass_product_coupled_yield(biomass, target, substrate)
        if self.manipulation_type == "genes":
            optimization_algorithm = GeneKnockoutOptimization(
                model=self._model,
                heuristic_method=self._algorithm,
                essential_genes=self._essential_genes,
                plot=self.plot,
                objective_function=objective_function,
                use_nullspace_simplification=use_nullspace_simplification)
        elif self.manipulation_type == "reactions":

            optimization_algorithm = ReactionKnockoutOptimization(