def _add_switches(self, reactions):
        logger.info("Adding switches.")
        y_vars = list()
        switches = list()
        self._exchanges = list()
        for reaction in reactions:
            if reaction.id.startswith('DM_'):
                # demand reactions don't need integer switches
                self._exchanges.append(reaction)
                continue

            y = self.model.solver.interface.Variable('y_' + reaction.id, lb=0, ub=1, type='binary')
            y_vars.append(y)
            # The following is a complicated but efficient way to write the following constraints

            # switch_lb = self.model.solver.interface.Constraint(y * reaction.lower_bound - reaction.flux_expression,
            #                                                    name='switch_lb_' + reaction.id, ub=0)
            # switch_ub = self.model.solver.interface.Constraint(y * reaction.upper_bound - reaction.flux_expression,
            #                                                    name='switch_ub_' + reaction.id, lb=0)
            forward_var_term = Mul._from_args((RealNumber(-1), reaction.forward_variable))
            reverse_var_term = Mul._from_args((RealNumber(-1), reaction.reverse_variable))
            switch_lb_y_term = Mul._from_args((RealNumber(reaction.lower_bound), y))
            switch_ub_y_term = Mul._from_args((RealNumber(reaction.upper_bound), y))
            switch_lb = self.model.solver.interface.Constraint(
                Add._from_args((switch_lb_y_term, forward_var_term, reverse_var_term)), name='switch_lb_' + reaction.id,
                ub=0, sloppy=True)
            switch_ub = self.model.solver.interface.Constraint(
                Add._from_args((switch_ub_y_term, forward_var_term, reverse_var_term)), name='switch_ub_' + reaction.id,
                lb=0, sloppy=True)
            switches.extend([switch_lb, switch_ub])
        self.model.solver.add(y_vars)
        self.model.solver.add(switches, sloppy=True)
        logger.info("Setting minimization of switch variables as objective.")
        self.model.objective = self.model.solver.interface.Objective(Add(*y_vars), direction='min')
        self._y_vars_ids = [var.name for var in y_vars]

def test_lookup_table():
    from random import uniform, randrange
    from sympy import Add
    from sympy.integrals.meijerint import z as z_dummy

    table = {}
    _create_lookup_table(table)
    for _, l in sorted(table.items()):
        for formula, terms, cond, hint in sorted(l, key=default_sort_key):
            subs = {}
            for a in list(formula.free_symbols) + [z_dummy]:
                if hasattr(a, "properties") and a.properties:
                    # these Wilds match positive integers
                    subs[a] = randrange(1, 10)
                else:
                    subs[a] = uniform(1.5, 2.0)
            if not isinstance(terms, list):
                terms = terms(subs)

            # First test that hyperexpand can do this.
            expanded = [hyperexpand(g) for (_, g) in terms]
            assert all(x.is_Piecewise or not x.has(meijerg) for x in expanded)

            # Now test that the meijer g-function is indeed as advertised.
            expanded = Add(*[f * x for (f, x) in terms])
            a, b = formula.n(subs=subs), expanded.n(subs=subs)
            r = min(abs(a), abs(b))
            if r < 1:
                assert abs(a - b).n() <= 1e-10
            else:
                assert (abs(a - b) / r).n() <= 1e-10

def test_gcd_terms():
    f = 2*(x + 1)*(x + 4)/(5*x**2 + 5) + (2*x + 2)*(x + 5)/(x**2 + 1)/5 + (2*x + 2)*(x + 6)/(5*x**2 + 5)

    assert _gcd_terms(f) == ((S(6)/5)*((1 + x)/(1 + x**2)), 5 + x, 1)
    assert _gcd_terms(Add.make_args(f)) == ((S(6)/5)*((1 + x)/(1 + x**2)), 5 + x, 1)

    assert gcd_terms(f) == (S(6)/5)*((1 + x)*(5 + x)/(1 + x**2))
    assert gcd_terms(Add.make_args(f)) == (S(6)/5)*((1 + x)*(5 + x)/(1 + x**2))

    assert gcd_terms((2*x + 2)**3 + (2*x + 2)**2) == 4*(x + 1)**2*(2*x + 3)

    assert gcd_terms(0) == 0
    assert gcd_terms(1) == 1
    assert gcd_terms(x) == x
    assert gcd_terms(2 + 2*x) == Mul(2, 1 + x, evaluate=False)
    arg = x*(2*x + 4*y)
    garg = 2*x*(x + 2*y)
    assert gcd_terms(arg) == garg
    assert gcd_terms(sin(arg)) == sin(garg)

    # issue 3040-like
    alpha, alpha1, alpha2, alpha3 = symbols('alpha:4')
    a = alpha**2 - alpha*x**2 + alpha + x**3 - x*(alpha + 1)
    rep = (alpha, (1 + sqrt(5))/2 + alpha1*x + alpha2*x**2 + alpha3*x**3)
    s = (a/(x - alpha)).subs(*rep).series(x, 0, 1)
    assert simplify(collect(s, x)) == -sqrt(5)/2 - S(3)/2 + O(x)

    # issue 2818
    assert _gcd_terms([S.Zero, S.Zero]) == (0, 0, 1)
    assert _gcd_terms([2*x + 4]) == (2, x + 2, 1)

def test_as_ordered_terms():
    f, g = symbols("f,g", cls=Function)

    assert x.as_ordered_terms() == [x]
    assert (sin(x) ** 2 * cos(x) + sin(x) * cos(x) ** 2 + 1).as_ordered_terms() == [
        sin(x) ** 2 * cos(x),
        sin(x) * cos(x) ** 2,
        1,
    ]

    args = [f(1), f(2), f(3), f(1, 2, 3), g(1), g(2), g(3), g(1, 2, 3)]
    expr = Add(*args)

    assert expr.as_ordered_terms() == args

    assert (1 + 4 * sqrt(3) * pi * x).as_ordered_terms() == [4 * pi * x * sqrt(3), 1]

    assert (2 + 3 * I).as_ordered_terms() == [2, 3 * I]
    assert (-2 + 3 * I).as_ordered_terms() == [-2, 3 * I]
    assert (2 - 3 * I).as_ordered_terms() == [2, -3 * I]
    assert (-2 - 3 * I).as_ordered_terms() == [-2, -3 * I]

    assert (4 + 3 * I).as_ordered_terms() == [4, 3 * I]
    assert (-4 + 3 * I).as_ordered_terms() == [-4, 3 * I]
    assert (4 - 3 * I).as_ordered_terms() == [4, -3 * I]
    assert (-4 - 3 * I).as_ordered_terms() == [-4, -3 * I]

    f = x ** 2 * y ** 2 + x * y ** 4 + y + 2

    assert f.as_ordered_terms(order="lex") == [x ** 2 * y ** 2, x * y ** 4, y, 2]
    assert f.as_ordered_terms(order="grlex") == [x * y ** 4, x ** 2 * y ** 2, y, 2]
    assert f.as_ordered_terms(order="rev-lex") == [2, y, x * y ** 4, x ** 2 * y ** 2]
    assert f.as_ordered_terms(order="rev-grlex") == [2, y, x ** 2 * y ** 2, x * y ** 4]

def cancel_terms(sym, x_term, coef):
    if coef.is_Add:
        for arg_c in coef.args:
            sym = cancel_terms(sym, x_term, arg_c)
    else:
        terms = Add.make_args(sym)
        return Add.fromiter(t for t in terms if t != x_term*coef)

def test_combine_inverse():
    x, y = symbols("x y")
    assert Mul._combine_inverse(x*I*y, x*I) == y
    assert Mul._combine_inverse(x*I*y, y*I) == x
    assert Mul._combine_inverse(oo*I*y, y*I) == oo
    assert Mul._combine_inverse(oo*I*y, oo*I) == y
    assert Add._combine_inverse(oo, oo) == S(0)
    assert Add._combine_inverse(oo*I, oo*I) == S(0)

    def atomic_ordering_energy(self, dbe):
        """
        Return the atomic ordering contribution in symbolic form.
        Description follows Servant and Ansara, Calphad, 2001.
        """
        phase = dbe.phases[self.phase_name]
        ordered_phase_name = phase.model_hints.get("ordered_phase", None)
        disordered_phase_name = phase.model_hints.get("disordered_phase", None)
        if phase.name != ordered_phase_name:
            return S.Zero
        disordered_model = self.__class__(dbe, self.components, disordered_phase_name)
        constituents = [
            sorted(set(c).intersection(self.components)) for c in dbe.phases[ordered_phase_name].constituents
        ]

        # Fix variable names
        variable_rename_dict = {}
        for atom in disordered_model.energy.atoms(v.SiteFraction):
            # Replace disordered phase site fractions with mole fractions of
            # ordered phase site fractions.
            # Special case: Pure vacancy sublattices
            all_species_in_sublattice = dbe.phases[disordered_phase_name].constituents[atom.sublattice_index]
            if atom.species == "VA" and len(all_species_in_sublattice) == 1:
                # Assume: Pure vacancy sublattices are always last
                vacancy_subl_index = len(dbe.phases[ordered_phase_name].constituents) - 1
                variable_rename_dict[atom] = v.SiteFraction(ordered_phase_name, vacancy_subl_index, atom.species)
            else:
                # All other cases: replace site fraction with mole fraction
                variable_rename_dict[atom] = self.mole_fraction(
                    atom.species, ordered_phase_name, constituents, dbe.phases[ordered_phase_name].sublattices
                )
        # Save all of the ordered energy contributions
        # This step is why this routine must be called _last_ in build_phase
        ordered_energy = Add(*list(self.models.values()))
        self.models.clear()
        # Copy the disordered energy contributions into the correct bins
        for name, value in disordered_model.models.items():
            self.models[name] = value.xreplace(variable_rename_dict)
        # All magnetic parameters will be defined in the disordered model
        self.TC = self.curie_temperature = disordered_model.TC
        self.TC = self.curie_temperature = self.TC.xreplace(variable_rename_dict)

        molefraction_dict = {}

        # Construct a dictionary that replaces every site fraction with its
        # corresponding mole fraction in the disordered state
        for sitefrac in ordered_energy.atoms(v.SiteFraction):
            all_species_in_sublattice = dbe.phases[ordered_phase_name].constituents[sitefrac.sublattice_index]
            if sitefrac.species == "VA" and len(all_species_in_sublattice) == 1:
                # pure-vacancy sublattices should not be replaced
                # this handles cases like AL,NI,VA:AL,NI,VA:VA and
                # ensures the VA's don't get mixed up
                continue
            molefraction_dict[sitefrac] = self.mole_fraction(
                sitefrac.species, ordered_phase_name, constituents, dbe.phases[ordered_phase_name].sublattices
            )

        return ordered_energy - ordered_energy.subs(molefraction_dict, simultaneous=True)

def test_gcd_terms():
    f = 2*(x + 1)*(x + 4)/(5*x**2 + 5) + (2*x + 2)*(x + 5)/(x**2 + 1)/5 + \
        (2*x + 2)*(x + 6)/(5*x**2 + 5)

    assert _gcd_terms(f) == ((S(6)/5)*((1 + x)/(1 + x**2)), 5 + x, 1)
    assert _gcd_terms(Add.make_args(f)) == \
        ((S(6)/5)*((1 + x)/(1 + x**2)), 5 + x, 1)

    newf = (S(6)/5)*((1 + x)*(5 + x)/(1 + x**2))
    assert gcd_terms(f) == newf
    args = Add.make_args(f)
    # non-Basic sequences of terms treated as terms of Add
    assert gcd_terms(list(args)) == newf
    assert gcd_terms(tuple(args)) == newf
    assert gcd_terms(set(args)) == newf
    # but a Basic sequence is treated as a container
    assert gcd_terms(Tuple(*args)) != newf
    assert gcd_terms(Basic(Tuple(1, 3*y + 3*x*y), Tuple(1, 3))) == \
        Basic((1, 3*y*(x + 1)), (1, 3))
    # but we shouldn't change keys of a dictionary or some may be lost
    assert gcd_terms(Dict((x*(1 + y), 2), (x + x*y, y + x*y))) == \
        Dict({x*(y + 1): 2, x + x*y: y*(1 + x)})

    assert gcd_terms((2*x + 2)**3 + (2*x + 2)**2) == 4*(x + 1)**2*(2*x + 3)

    assert gcd_terms(0) == 0
    assert gcd_terms(1) == 1
    assert gcd_terms(x) == x
    assert gcd_terms(2 + 2*x) == Mul(2, 1 + x, evaluate=False)
    arg = x*(2*x + 4*y)
    garg = 2*x*(x + 2*y)
    assert gcd_terms(arg) == garg
    assert gcd_terms(sin(arg)) == sin(garg)

    # issue 6139-like
    alpha, alpha1, alpha2, alpha3 = symbols('alpha:4')
    a = alpha**2 - alpha*x**2 + alpha + x**3 - x*(alpha + 1)
    rep = (alpha, (1 + sqrt(5))/2 + alpha1*x + alpha2*x**2 + alpha3*x**3)
    s = (a/(x - alpha)).subs(*rep).series(x, 0, 1)
    assert simplify(collect(s, x)) == -sqrt(5)/2 - S(3)/2 + O(x)

    # issue 5917
    assert _gcd_terms([S.Zero, S.Zero]) == (0, 0, 1)
    assert _gcd_terms([2*x + 4]) == (2, x + 2, 1)

    eq = x/(x + 1/x)
    assert gcd_terms(eq, fraction=False) == eq
    eq = x/2/y + 1/x/y
    assert gcd_terms(eq, fraction=True, clear=True) == \
        (x**2 + 2)/(2*x*y)
    assert gcd_terms(eq, fraction=True, clear=False) == \
        (x**2/2 + 1)/(x*y)
    assert gcd_terms(eq, fraction=False, clear=True) == \
        (x + 2/x)/(2*y)
    assert gcd_terms(eq, fraction=False, clear=False) == \
        (x/2 + 1/x)/y

def test__combine_inverse():
    x, y = symbols("x y")
    assert Mul._combine_inverse(x*I*y, x*I) == y
    assert Mul._combine_inverse(x*I*y, y*I) == x
    assert Mul._combine_inverse(oo*I*y, y*I) == oo
    assert Mul._combine_inverse(oo*I*y, oo*I) == y
    assert Add._combine_inverse(oo, oo) == S(0)
    assert Add._combine_inverse(oo*I, oo*I) == S(0)
    assert Add._combine_inverse(x*oo, x*oo) == S(0)
    assert Add._combine_inverse(-x*oo, -x*oo) == S(0)
    assert Add._combine_inverse((x - oo)*(x + oo), -oo)

def test_make_args():
    assert Add.make_args(x) == (x,)
    assert Mul.make_args(x) == (x,)

    assert Add.make_args(x*y*z) == (x*y*z,)
    assert Mul.make_args(x*y*z) == (x*y*z).args

    assert Add.make_args(x + y + z) == (x + y + z).args
    assert Mul.make_args(x + y + z) == (x + y + z,)

    assert Add.make_args((x + y)**z) == ((x + y)**z,)
    assert Mul.make_args((x + y)**z) == ((x + y)**z,)

def test_as_ordered_terms():
    f, g = symbols('f,g', cls=Function)

    assert x.as_ordered_terms() == [x]
    assert (sin(x)**2*cos(x) + sin(x)*cos(x)**2 + 1).as_ordered_terms() == [sin(x)**2*cos(x), sin(x)*cos(x)**2, 1]

    expr = Add(*[f(1), f(2), f(3), f(1, 2, 3), g(1), g(2), g(3), g(1, 2, 3)])

    assert expr.as_ordered_terms() == \
        [f(1), f(2), f(3), f(1, 2, 3), g(1), g(2), g(3), g(1, 2, 3)]

    assert (1 + 4*sqrt(3)*pi*x).as_ordered_terms() == [4*pi*x*sqrt(3), 1]

def test_gcd_terms():
    f = 2*(x + 1)*(x + 4)/(5*x**2 + 5) + (2*x + 2)*(x + 5)/(x**2 + 1)/5 + (2*x + 2)*(x + 6)/(5*x**2 + 5)

    assert _gcd_terms(f) == ((S(6)/5)*((1 + x)/(1 + x**2)), 5 + x, 1)
    assert _gcd_terms(Add.make_args(f)) == ((S(6)/5)*((1 + x)/(1 + x**2)), 5 + x, 1)

    assert gcd_terms(f) == (S(6)/5)*((1 + x)*(5 + x)/(1 + x**2))
    assert gcd_terms(Add.make_args(f)) == (S(6)/5)*((1 + x)*(5 + x)/(1 + x**2))

    assert gcd_terms(0) == 0
    assert gcd_terms(1) == 1
    assert gcd_terms(x) == x

 def _add_reaction_dual_constraint(self, reaction, coefficient, maximization, prefix):
     """Add a dual constraint corresponding to the reaction's objective coefficient"""
     stoichiometry = {self.solver.variables["lambda_" + m.id]: c for m, c in six.iteritems(reaction.metabolites)}
     if maximization:
         constraint = self.solver.interface.Constraint(
             Add._from_args(tuple(c * v for v, c in six.iteritems(stoichiometry))),
             name="r_%s_%s" % (reaction.id, prefix),
             lb=coefficient)
     else:
         constraint = self._dual_solver.interface.Constraint(
             Add._from_args(tuple(c * v for v, c in six.iteritems(stoichiometry))),
             name="r_%s_%s" % (reaction.id, prefix),
             ub=coefficient)
     self.solver._add_constraint(constraint)

def _inverse_mellin_transform(F, s, x_, strip, as_meijerg=False):
    """ A helper for the real inverse_mellin_transform function, this one here
        assumes x to be real and positive. """
    from sympy import (expand, expand_mul, hyperexpand, meijerg, And, Or,
                       arg, pi, re, factor, Heaviside, gamma, Add)
    x = _dummy('t', 'inverse-mellin-transform', F, positive=True)
    # Actually, we won't try integration at all. Instead we use the definition
    # of the Meijer G function as a fairly general inverse mellin transform.
    F = F.rewrite(gamma)
    for g in [factor(F), expand_mul(F), expand(F)]:
        if g.is_Add:
            # do all terms separately
            ress = [_inverse_mellin_transform(G, s, x, strip, as_meijerg,
                                              noconds=False) \
                    for G in g.args]
            conds = [p[1] for p in ress]
            ress = [p[0] for p in ress]
            res = Add(*ress)
            if not as_meijerg:
                res = factor(res, gens=res.atoms(Heaviside))
            return res.subs(x, x_), And(*conds)

        try:
            a, b, C, e, fac = _rewrite_gamma(g, s, strip[0], strip[1])
        except IntegralTransformError:
            continue
        G = meijerg(a, b, C/x**e)
        if as_meijerg:
            h = G
        else:
            h = hyperexpand(G)
            if h.is_Piecewise and len(h.args) == 3:
                # XXX we break modularity here!
                h = Heaviside(x - abs(C))*h.args[0].args[0] \
                  + Heaviside(abs(C) - x)*h.args[1].args[0]
        # We must ensure that the intgral along the line we want converges,
        # and return that value.
        # See [L], 5.2
        cond = [abs(arg(G.argument)) < G.delta*pi]
        # Note: we allow ">=" here, this corresponds to convergence if we let
        # limits go to oo symetrically. ">" corresponds to absolute convergence.
        cond += [And(Or(len(G.ap) != len(G.bq), 0 >= re(G.nu) + 1),
                     abs(arg(G.argument)) == G.delta*pi)]
        cond = Or(*cond)
        if cond is False:
            raise IntegralTransformError('Inverse Mellin', F, 'does not converge')
        return (h*fac).subs(x, x_), cond

    raise IntegralTransformError('Inverse Mellin', F, '')

def sampling_E(condition, given=None, numsamples=1, evalf=True, **kwargs):
    """
    Sampling version of E

    See Also
    ========
    P
    sampling_P
    """

    samples = sample_iter(condition, given, numsamples=numsamples, **kwargs)

    result = Add(*list(samples)) / numsamples
    if evalf:   return result.evalf()
    else:       return result

def test_gcd_terms():
    f = 2*(x + 1)*(x + 4)/(5*x**2 + 5) + (2*x + 2)*(x + 5)/(x**2 + 1)/5 + (2*x + 2)*(x + 6)/(5*x**2 + 5)

    assert _gcd_terms(f) == ((S(6)/5)*((1 + x)/(1 + x**2)), 5 + x, 1)
    assert _gcd_terms(Add.make_args(f)) == ((S(6)/5)*((1 + x)/(1 + x**2)), 5 + x, 1)

    assert gcd_terms(f) == (S(6)/5)*((1 + x)*(5 + x)/(1 + x**2))
    assert gcd_terms(Add.make_args(f)) == (S(6)/5)*((1 + x)*(5 + x)/(1 + x**2))

    assert gcd_terms((2*x + 2)**3 + (2*x + 2)**2) == 4*(x + 1)**2*(2*x + 3)

    assert gcd_terms(0) == 0
    assert gcd_terms(1) == 1
    assert gcd_terms(x) == x
    assert gcd_terms(2 + 2*x) == Mul(2, 1 + x, evaluate=False)

 def calc_part(expr, nexpr):
     from sympy import Add
     nint = int(to_int(nexpr, rnd))
     n, c, p, b = nexpr
     if (c != 1 and p != 0) or p < 0:
         expr = Add(expr, -nint, evaluate=False)
         x, _, x_acc, _ = evalf(expr, 10, options)
         try:
             check_target(expr, (x, None, x_acc, None), 3)
         except PrecisionExhausted:
             if not expr.equals(0):
                 raise PrecisionExhausted
             x = fzero
         nint += int(no*(mpf_cmp(x or fzero, fzero) == no))
     nint = from_int(nint)
     return nint, fastlog(nint) + 10

    def __new__(cls, *args):

        args = map(matrixify, args)

        args = [arg for arg in args if arg!=0]

        if not all(arg.is_Matrix for arg in args):
            raise ValueError("Mix of Matrix and Scalar symbols")

        # Check that the shape of the args is consistent
        A = args[0]
        for B in args[1:]:
            if A.shape != B.shape:
                raise ShapeError("Matrices %s and %s are not aligned"%(A,B))

        expr = Add.__new__(cls, *args)
        if expr == S.Zero:
            return ZeroMatrix(*args[0].shape)
        expr = matrixify(expr)

        if expr.is_Mul:
            return MatMul(*expr.args)

        # Clear out Identities
        # Any zeros around?
        if expr.is_Add and any(M.is_ZeroMatrix for M in expr.args):
            newargs = [M for M in expr.args if not M.is_ZeroMatrix] # clear out
            if len(newargs)==0: # Did we lose everything?
                return ZeroMatrix(*args[0].shape)
            if expr.args != newargs: # Removed some 0's but not everything?
                return MatAdd(*newargs) # Repeat with simpler expr

        return expr

    def _as_ordered_terms(self, expr, order=None):
        """A compatibility function for ordering terms in Add. """
        order = order or self.order

        if order == 'old':
            return sorted(Add.make_args(expr), key=cmp_to_key(Basic._compare_pretty))
        else:
            return expr.as_ordered_terms(order=order)

 def _print_Add(self, expr):
     if len(expr.args) != 2:
         return "add({}, {})".format(
             self._print(expr.args[0]),
             self._print(Add.fromiter(expr.args[1:]))
         )
     return "add({}, {})".format(
         self._print(expr.args[0]),
         self._print(expr.args[1]),
     )