class Result(object):
    def __init__(self, *args, **kwargs):
        super(Result, self).__init__(*args, **kwargs)
        self._meta_information = MetaInformation()

    @property
    def meta_information(self):
        return self._meta_information

    @property
    def data_frame(self):
        raise NotImplementedError

    def plot(self, grid=None, width=None, height=None, title=None, *args, **kwargs):
        raise NotImplementedError def fba(model, objective=None, reactions=None, *args, **kwargs):
def pfba(model, objective=None, reactions=None, fraction_of_optimum=1, *args, **kwargs):
def moma(model, reference=None, cache=None, reactions=None, *args, **kwargs):
            def create_variable(model, variable_id):
            def create_constraint(model, constraint_id, var, reaction, flux_value):
            def update_constraint(model, constraint, var, reaction, flux_value):
        def create_objective(model, variables):
def lmoma(model, reference=None, cache=None, reactions=None, *args, **kwargs):
            def create_variable(model, var_id, lb):
            def update_upper_constraint(model, constraint, var, reaction, flux_value):
            def create_upper_constraint(model, constraint_id, var, reaction, flux_value):
            def update_lower_constraint(model, constraint, var, reaction, flux_value):
            def create_lower_constraint(model, constraint_id, var, reaction, flux_value):
        def create_objective(model, variables):
def room(model, reference=None, cache=None, delta=0.03, epsilon=0.001, reactions=None, *args, **kwargs):
            def create_variable(model, var_id):
            def create_upper_constraint(model, constraint_id, reaction, variable, flux_value, epsilon):
            def update_upper_constraint(model, constraint, reaction, variable, flux_value, epsilon):
            def create_lower_constraint(model, constraint_id, reaction, variable, flux_value, epsilon):
            def update_lower_constraint(model, constraint, reaction, variable, flux_value, epsilon):
    def from_solution(cls, solution, *args, **kwargs):
    def __init__(self, fluxes, objective_value, *args, **kwargs):
    def __getitem__(self, item):
    def data_frame(self):
    def fluxes(self):
    def objective_value(self):
    def plot(self, grid=None, width=None, height=None, title=None):
    def iteritems(self):
    def items(self):
    def keys(self):
    def values(self):
    def _repr_html_(self):
    def plot_scale(self, values, palette="YlGnBu"):
    def display_on_map(self, map_name=None, palette="YlGnBu"):
class FluxDistributionResult(Result): def find_blocked_reactions_nullspace(model, ns=None, tol=1e-10):
    """Identify reactions that can't carry flux based on the nullspace of the stoichiometric matrix.
    A blocked reaction will correspond to an all-zero row in N(S)."""
    if ns is None:
        ns = nullspace(create_stoichiometric_array(model))
    mask = (np.abs(ns) <= tol).all(1)
    blocked = frozenset(reac for reac, b in zip(model.reactions, mask) if b is True)
    return blocked def find_coupled_reactions_nullspace(model, ns=None, tol=1e-10):
    """
    Find groups of reactions whose fluxes are forced to be multiples of each other.

    Parameters
    ----------
    model: SolverBasedModel
        A constraint-based model.
    ns: numpy.array
        The nullspace if already computed somewhere else.
    tol: float
        Value after which a value is treated as zero.

    Returns
    -------
    list:
        A list with groups: dictionaries {reaction: relative_coefficient}

    """
    if ns is None:
        ns = nullspace(create_stoichiometric_array(model))
    mask = (np.abs(ns) <= tol).all(axis=1)  # Mask for blocked reactions
    non_blocked_ns = ns[~mask]
    non_blocked_reactions = np.array(list(model.reactions))[~mask]

    corr_mat = np.corrcoef(non_blocked_ns)
    dist_mat = 1 - np.abs(corr_mat)
    groups = []

    reaction_index = {r: i for i, r in enumerate(non_blocked_reactions)}

    for i, reaction_i in enumerate(non_blocked_reactions):
        left = non_blocked_ns[i]
        group = next((g for g in groups if reaction_i in g), None)
        if group:
            reaction_i = next(reac for reac, c in group.items() if c == 1)
            left = non_blocked_ns[reaction_index[reaction_i]]
        else:
            group = {reaction_i: 1}
            groups.append(group)

        good_corr_index = np.argwhere(dist_mat[i] < tol).reshape(-1)
        for j in good_corr_index:
            if j == i:
                continue
            reaction_j = non_blocked_reactions[j]
            right = non_blocked_ns[j]
            ratio = np.apply_along_axis(lambda x: x[0] / x[1] if abs(x[1]) > 0. else np.inf, 0, np.array([left, right]))

            # special case:
            # if ratio is 1 (a/b == 1) then a == b.
            # but if a = 0 and b = 0, then a/b = np.inf.
            # solution:
            # mask inf from ratio
            # check if non-inf elements ratio is ~1
            # check if left and right values are 0 for indices with inf ratio
            # if yes, replace with 1
            inf_mask = np.isinf(ratio)
            non_inf = ratio[~inf_mask]

            if (abs((non_inf - 1)) < tol * 100).all():
                right_is_zero = (abs(right[inf_mask]) < tol).all()
                left_is_zero = (abs(left[inf_mask]) < tol).all()

                if right_is_zero and left_is_zero:
                    ratio[inf_mask] = 1

            if abs(max(ratio) - min(ratio)) < tol * 100:
                group[reaction_j] = round(ratio.mean(), 10)

    groups = [g for g in groups if len(g) > 1]

    return groups def nullspace(matrix, atol=1e-13, rtol=0):
    """Compute an approximate basis for the nullspace of A.

    The algorithm used by this function is based on the singular value
    decomposition of `A`.

    Parameters
    ----------
    matrix : ndarray
        A should be at most 2-D.  A 1-D array with length k will be treated
        as a 2-D with shape (1, k)
    atol : float
        The absolute tolerance for a zero singular value.  Singular values
        smaller than `atol` are considered to be zero.
    rtol : float
        The relative tolerance.  Singular values less than rtol*smax are
        considered to be zero, where smax is the largest singular value.

    If both `atol` and `rtol` are positive, the combined tolerance is the
    maximum of the two; that is::
        tol = max(atol, rtol * smax)
    Singular values smaller than `tol` are considered to be zero.

    Return value
    ------------
    ns : ndarray
        If `A` is an array with shape (m, k), then `ns` will be an array
        with shape (k, n), where n is the estimated dimension of the
        nullspace of `A`.  The columns of `ns` are a basis for the
        nullspace; each element in numpy.dot(A, ns) will be approximately
        zero.
    """
    matrix = np.atleast_2d(matrix)
    u, s, vh = svd(matrix)
    tol = max(atol, rtol * s[0])
    nnz = (s >= tol).sum()
    ns = vh[nnz:].conj().T
    return ns def create_stoichiometric_array(model, array_type='dense', dtype=None):
    """Return a stoichiometric array representation of the given model.
    The the columns represent the reactions and rows represent
    metabolites. S[i,j] therefore contains the quantity of metabolite `i`
    produced (negative for consumed) by reaction `j`.
    Parameters
    ----------
    model : cobra.Model
        The cobra model to construct the matrix for.
    array_type : string
        The type of array to construct. if 'dense', return a standard
        numpy.array, 'dok', or 'lil' will construct a sparse array using
        scipy of the corresponding type and 'data_frame' will give a
        pandas `DataFrame` with metabolite and reaction identifiers as indices.
    dtype : data-type
        The desired data-type for the array. If not given, defaults to float.
    Returns
    -------
    matrix of class `dtype`
        The stoichiometric matrix for the given model.
    """

    if dtype is None:
        dtype = np.float64

    def data_frame(_, dtype):
        metabolite_ids = [met.id for met in model.metabolites]
        reaction_ids = [rxn.id for rxn in model.reactions]
        index = pandas.MultiIndex.from_tuples(
            list(product(metabolite_ids, reaction_ids)))
        return pandas.DataFrame(data=0, index=index, columns=['stoichiometry'], dtype=dtype)

    array_constructor = {
        'dense': np.zeros, 'dok': dok_matrix, 'lil': lil_matrix,
        'data_frame': data_frame
    }

    n_metabolites = len(model.metabolites)
    n_reactions = len(model.reactions)
    array = array_constructor[array_type]((n_metabolites, n_reactions), dtype=dtype)

    m_ind = model.metabolites.index
    r_ind = model.reactions.index

    for reaction in model.reactions:
        for metabolite, stoich in reaction.metabolites.items():
            if array_type == 'data_frame':
                array.set_value((metabolite.id, reaction.id), 'stoichiometry', stoich)
            else:
                array[m_ind(metabolite), r_ind(reaction)] = stoich

    return array class BestSolutionArchive(object):
class ProductionStrainArchive(BestSolutionArchive):
class Individual(OriginalIndividual):
    def __init__(self):
    def __call__(self, random, population, archive, args):
    def add(self, candidate, fitness, birthdate, maximize, max_archive_size):
    def length(self):
    def get(self, index):
    def __iter__(self):
    def __len__(self):
    def reset(self):
    def __call__(self, random, population, archive, args):
    def __init__(self, candidate, fitness, maximize=True, birthdate=time.time()):
    def __iter__(self):
    def __eq__(self, other):
    def __cmp__(self, other):
    def __lt__(self, other):
    def __ge__(self, other):
    def __le__(self, other):
    def __gt__(self, other):
    def __str__(self):
    def __repr__(self):
    def issubset(self, other):
    def symmetric_difference(self, other):
    def improves(self, other): class SetDecoder(object):
class ReactionSetDecoder(SetDecoder):
class GeneSetDecoder(SetDecoder):
    def __init__(self, representation, model, *args, **kwargs):
    def __call__(self, individual, flat=False, decompose=False):
    def __init__(self, representation, model, groups=None, *args, **kwargs):
    def __call__(self, individual, flat=False, decompose=False):
    def __init__(self, representation, model, groups=None, *args, **kwargs):
    def __call__(self, individual, flat=False, decompose=False): class Evaluator(object):
class TargetEvaluator(Evaluator):
class KnockoutEvaluator(TargetEvaluator):
class SwapEvaluator(TargetEvaluator):
class KnockinKnockoutEvaluator(KnockoutEvaluator):
    def __call__(self, population):
    def evaluate_individual(self, individual):
    def reset(self):
    def __init__(self, model, decoder, objective_function, simulation_method, simulation_kwargs):
    def __call__(self, population):
    def reset(self):
    def evaluate_individual(self, individual):
    def __init__(self, swap_pair=None, *args, **kwargs):
    def evaluate_individual(self, individual):
    def __init__(self, database, *args, **kwargs): def set_generator(random, args):
def multiple_chromosome_set_generator(random, args):
def linear_set_generator(random, args): class ProgressObserver(object):
    def __init__(self):
    def __call__(self, population, num_generations, num_evaluations, args):
    def reset(self):
    def end(self): class IPythonBokehFitnessPlotter(object):
class IPythonBokehParetoPlotter(object):
class GeneFrequencyPlotter(object):
    def __init__(self, window_size=1000):
    def _set_plot(self):
    def __call__(self, population, num_generations, num_evaluations, args):
    def _update(self):
    def reset(self):
    def end(self):
    def __init__(self, objective_function=None, x=0, y=1):
    def _set_plot(self):
    def __call__(self, population, num_generations, num_evaluations, args):
    def _update(self):
    def reset(self):
    def end(self):
    def __init__(self, solutions, url='default'):
    def frequencies(self):
    def plot(self): class GenericStatsData(object):
    class CLIStatsData(GenericStatsData):
class BokehStatsData(GenericStatsData):
    def __init__(self, solution, *args, **kwargs):
    def display(self):
        def __init__(self, *args, **kwargs):
        def display(self):
    def __init__(self, *args, **kwargs):
    def display(self): def _subtract(list_a, list_b):
def _do_set_n_point_crossover(representation, mom, dad, points, random, max_size):
def set_n_point_crossover(random, mom, dad, args):
def set_mutation(random, individual, args):
def set_indel(random, individual, args):
def multiple_chromosome_set_mutation(random, individual, args):
def multiple_chromosome_set_indel(random, individual, args):
def multiple_chromosome_n_point_crossover(random, mom, dad, args): class Individual(OriginalIndividual):
    def __init__(self, candidate, fitness, maximize=True, birthdate=time.time()):
        super(Individual, self).__init__(set(candidate), maximize)
        self.fitness = fitness
        self.birthdate = birthdate

    def __iter__(self):
        for v in self.candidate:
            yield v

    def __eq__(self, other):
        return self.candidate == other.candidate and self.fitness == other.fitness

    def __cmp__(self, other):
        if self.fitness > other.fitness:
            return -1 if self.maximize else 1
        elif self.fitness == other.fitness:
            if self.improves(other):
                return -1
            elif self == other:
                return 0
            else:
                return 1
        else:
            return 1 if self.maximize else -1

    def __lt__(self, other):
        if self.fitness > other.fitness:
            return self.maximize
        elif self.fitness == other.fitness:
            if self.improves(other):
                return True
            elif self == other:
                return False
            else:
                return False
        else:
            return not self.maximize

    def __ge__(self, other):
        return self.__eq__(other) or self.__gt__(other)

    def __le__(self, other):
        return self.__eq__(other) or self.__lt__(other)

    def __gt__(self, other):
        if self.fitness > other.fitness:
            return not self.maximize
        elif self.fitness == other.fitness:
            if self.improves(other):
                return False
            elif self == other:
                return False
            else:
                return True
        else:
            return self.maximize

    def __str__(self):
        sense = "max" if self.maximize else "min"
        return "%s - %s sense: %s" % (list(self.candidate), self.fitness, sense)

    def __repr__(self):
        return "Individual #%s: %s" % (id(self), self.__str__())

    def issubset(self, other):
        return self.candidate.issubset(other.candidate)

    def symmetric_difference(self, other):
        return self.candidate.symmetric_difference(other.candidate)

    def improves(self, other):
        assert isinstance(other, Individual)
        if self.maximize:
            return self.issubset(other) and len(self.symmetric_difference(other)) > 0 and self.fitness >= other.fitness
        else:
            return self.issubset(other) and len(self.symmetric_difference(other)) > 0 and self.fitness <= other.fitness class MultiObjectiveFunction(ObjectiveFunction):
    def __init__(self, objectives, *args, **kwargs):
        super(MultiObjectiveFunction, self).__init__(*args, **kwargs)
        failed = []
        for i, objective in enumerate(objectives):
            if not isinstance(objective, ObjectiveFunction):
                failed.append(i)

        if len(failed) > 0:
            raise ValueError("Objectives %s are not instance of ObjectiveFunction" %
                             ", ".join(str(objectives[i]) for i in failed))

        self.objectives = objectives

    def __getitem__(self, item):
        return self.objectives[item]

    def __call__(self, model, solution, targets):
        return Pareto(values=[o(model, solution, targets) for o in self.objectives])

    def worst_fitness(self, maximize=True):
        return Pareto(values=[o.worst_fitness(maximize) for o in self.objectives], maximize=maximize)

    @property
    def reactions(self):
        return set(sum([o.reactions for o in self.objectives], []))

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

    @property
    def name(self):
        return "MO: " + "|".join([of.name for of in self.objectives])

    def _repr_latex_(self):
        return "\\begin{gather}\n" + "\\\\".join(o._repr_latex_() for o in self.objectives) + "\n\\end{gather}\n" class ObjectiveFunction(object):
    """
    Blueprint for objective function.

    All objective functions must extend this class and override the __call__ method.
    Because they implement the __call__ method, they will behave like functions.

    Methods
    -------
    __call__(model, solution, decoded_representation)
        Calculates the fitness of the solution

    """

    def __init__(self, *args, **kwargs):
        super(ObjectiveFunction, self).__init__(*args, **kwargs)

    def __call__(self, model, solution, decoded_representation):
        raise NotImplementedError

    @property
    def reactions(self):
        """
        Returns the reactions ids that require flux data for this evaluation function.
        """
        return []

    def _repr_latex_(self):
        return self.name

    @property
    def name(self):
        return self.__class__.__name__

    def worst_fitness(self, maximize=True):
        """

        Parameters
        ----------
        maximize: bool
            The sense of the optimization problem.

        Returns
        -------
        float
            The worst fitness value that a given objective can have based on the sense of the optimization.
        """
        raise NotImplementedError

    def __len__(self):
        return 1 class RandomGenerator(object):
    def __init__(self, seed=None):
        self._random = RandomState(seed=seed)

    def seed(self, seed):
        self._random = RandomState(seed=seed)

    def random(self):
        return self._random.rand()

    def randint(self, a, b=None):
        if b is None:
            b = a
            a = 0
        r = self._random.randint(a, high=b, size=1)
        return r[0]

    def sample(self, population, k):
        if k == 0:
            return []
        return list(self._random.choice(population, size=k, replace=False))

    def __getattr__(self, attr):
        return getattr(self._random, attr)

    def __getstate__(self):
        return {'_random': self._random}

    def __setstate__(self, d):
        self._random = d['_random']

    def uniform(self, low=0.0, high=1.0, size=None):
        return self._random.uniform(low, high, size) def reduce_reaction_set(reaction_set, groups):
    """
    Reduces a set of reactions according to a number of groups of reactions.
    The reduction will be performed so that the resulting set will contain no more than 1 reaction from each group.
    Reactions that are not in any of the groups will remain in the set.

    Parameters
    ----------
    reaction_set: Set
    groups: Iterable of sets

    Returns
    -------
    Set
    """
    reaction_set = set(reaction_set)  # Make a shallow copy
    result = []
    for group in groups:
        intersection = group & reaction_set
        if intersection:  # If any elements of group are in reaction_set, add one of these to result
            result.append(pick_one(intersection))
            reaction_set = reaction_set - intersection
    result = set(result) | reaction_set  # Add the remaining reactions to result
    return result def in_ipnb():
    """
    Check if it is running inside an IPython Notebook (updated for new notebooks)
    """
    return pandas.io.formats.console.in_ipython_frontend() def partition(ite, n):
    """Partition an iterable into n bite size chunks."""
    try:
        length = len(ite)
    except TypeError:
        ite = list(ite)
        length = len(ite)
    division = length / float(n)
    iterator = iter(ite)
    return [list(islice(iterator, 0, round(division * (i + 1)) - round(division * i))) for i in range(n)]
import collections
import logging
import time
import types
import inspyred
import numpy
from functools import reduce
from inspyred.ec.generators import diversify as diversify_function
from pandas import DataFrame
from cameo import config
from cameo.core.result import Result
from cameo.flux_analysis.simulation import pfba, lmoma, moma, room, logger as simulation_logger
from cameo.flux_analysis.structural import (find_blocked_reactions_nullspace, find_coupled_reactions_nullspace,
                                            nullspace,
                                            create_stoichiometric_array)
from cobra.flux_analysis import find_essential_genes, find_essential_reactions
from cameo.strain_design.heuristic.evolutionary import archives
from cameo.strain_design.heuristic.evolutionary import decoders
from cameo.strain_design.heuristic.evolutionary import evaluators
from cameo.strain_design.heuristic.evolutionary import generators
from cameo.strain_design.heuristic.evolutionary import observers
from cameo.strain_design.heuristic.evolutionary import plotters
from cameo.strain_design.heuristic.evolutionary import stats
from cameo.strain_design.heuristic.evolutionary import variators
from cameo.strain_design.heuristic.evolutionary.archives import Individual
from cameo.strain_design.heuristic.evolutionary.objective_functions import MultiObjectiveFunction, ObjectiveFunction
from cameo.util import RandomGenerator as Random, reduce_reaction_set
from cameo.util import in_ipnb
from cameo.util import partition
        wt_reference : dict, cameo.flux_analysis.simulation.FluxDistributionResult
            A dict (dict-like) object with flux values from a reference state.
        simulation_method : method
           the simulation method to use for evaluating results
        evaluator : TargetEvaluator
           the class used to evaluate results
        """
        super(TargetOptimization, self).__init__(*args, **kwargs)
        self._simulation_kwargs = dict()
        self._simulation_kwargs['reference'] = wt_reference
        self._simulation_method = None
        self.simulation_method = simulation_method
        self.representation = None
        self._evaluator = None
        self._target_type = None
        self._decoder = None
        self._metadata = {}

    @property
    def metadata(self):
        return self._metadata

    @property
    def simulation_method(self):
        return self._simulation_method

    @simulation_method.setter
    def simulation_method(self, simulation_method):
        if simulation_method in [lmoma, moma, room] and self._simulation_kwargs.get("reference", None) is None:
            logger.warning("No WT reference found, generating using pfba.")
            self._simulation_kwargs['reference'] = pfba(self.model).fluxes
            logger.warning("Reference successfully computed.")
        self._simulation_method = simulation_method

    @property
    def simulation_kwargs(self):
        return self._simulation_kwargs

    @simulation_kwargs.setter
    def simulation_kwargs(self, simulation_kwargs):
        if self.simulation_method in [lmoma, moma, room] and simulation_kwargs.get("reference", None) is None:
            logger.warning("No WT reference found, generating using pfba.")
            simulation_kwargs['reference'] = pfba(self.model).fluxes
            logger.warning("Reference successfully computed.")
        self._simulation_kwargs = simulation_kwargs

    @HeuristicOptimization.heuristic_method.setter
    def heuristic_method(self, heuristic_method):
        HeuristicOptimization.heuristic_method.fset(self, heuristic_method)
        self._set_observer()
        try:
            configuration = PRE_CONFIGURED[heuristic_method]
            self._setup(*configuration)
        except KeyError:
            logger.warning("Please verify the variator is compatible with set representation")

    def _setup(self, variator, selector, replacer, archiver):
        logger.debug("Setting up algorithm: %s" % self.heuristic_method)
        self.heuristic_method.variator = variator
        self.heuristic_method.selector = selector
        self.heuristic_method.replacer = replacer
        self.heuristic_method.archiver = archiver
        self.heuristic_method.terminator = self.termination

    @HeuristicOptimization.objective_function.setter
    def objective_function(self, objective_function):
        HeuristicOptimization.objective_function.fset(self, objective_function)
        self._set_observer()

    def _set_observer(self):
        self.observers = []

        if in_ipnb() and self.plot:
            if config.use_bokeh:
                if len(self.objective_function) > 1:
                    self.observers.append(plotters.IPythonBokehParetoPlotter(self.objective_function))
                else:
                    self.observers.append(plotters.IPythonBokehFitnessPlotter())

        else:
            if config.use_bokeh:
                pass
            else:
                pass
        if self.progress:
            self.observers.append(observers.ProgressObserver())

    def run(self, max_size=10, variable_size=True, diversify=False, view=config.default_view, **kwargs):
        """
        Runs the evolutionary algorithm.

        Parameters
        ----------
        max_size : int
            Maximum size of a solution, e.g., the maximum number of reactions or genes to knock-out or swap.
        variable_size : boolean
            If true, the solution size can change meaning that the combination of knockouts can have different sizes up
            to max_size. Otherwise it only produces knockout solutions with a fixed number of knockouts.
        diversify : bool
            It true, the generator will not be allowed to generate repeated candidates in the initial population.
        view : cameo.parallel.SequentialView, cameo.parallel.MultiprocessingView
            A view for single or multiprocessing.

        Returns
        -------
        TargetOptimizationResult
            The result of the optimization.
        """

        if kwargs.get('seed', None) is None:
            kwargs['seed'] = int(time.time())

        self.heuristic_method.observer = self.observers

        log_level = simulation_logger.level
        simulation_logger.setLevel(logging.CRITICAL)

        if diversify:
            generator = diversify_function(generators.set_generator)
        else:

            generator = generators.set_generator