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 quantity_simplify(expr):
    if expr.is_Atom:
        return expr
    if not expr.is_Mul:
        return expr.func(*map(quantity_simplify, expr.args))

    if expr.has(Prefix):
        coeff, args = expr.as_coeff_mul(Prefix)
        args = list(args)
        for arg in args:
            if isinstance(arg, Pow):
                coeff = coeff * (arg.base.scale_factor ** arg.exp)
            else:
                coeff = coeff * arg.scale_factor
        expr = coeff

    coeff, args = expr.as_coeff_mul(Quantity)
    args_pow = [arg.as_base_exp() for arg in args]
    quantity_pow, other_pow = sift(args_pow, lambda x: isinstance(x[0], Quantity), binary=True)
    quantity_pow_by_dim = sift(quantity_pow, lambda x: x[0].dimension)
    # Just pick the first quantity:
    ref_quantities = [i[0][0] for i in quantity_pow_by_dim.values()]
    new_quantities = [
        Mul.fromiter(
            (quantity*i.scale_factor/quantity.scale_factor)**p for i, p in v)
            if len(v) > 1 else v[0][0]**v[0][1]
        for quantity, (k, v) in zip(ref_quantities, quantity_pow_by_dim.items())]
    return coeff*Mul.fromiter(other_pow)*Mul.fromiter(new_quantities)

 def t(a, b, arg, n):
     from sympy import Mul
     m1 = meijerg(a, b, arg)
     m2 = Mul(*_inflate_g(m1, n))
     # NOTE: (the random number)**9 must still be on the principal sheet.
     # Thus make b&d small to create random numbers of small imaginary part.
     return verify_numerically(m1.subs(subs), m2.subs(subs), x, b=0.1, d=-0.1)

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 test_as_ordered_factors():
    f, g = symbols("f,g", cls=Function)

    assert x.as_ordered_factors() == [x]
    assert (2 * x * x ** n * sin(x) * cos(x)).as_ordered_factors() == [Integer(2), x, x ** n, sin(x), cos(x)]

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

    assert expr.as_ordered_factors() == args

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 remove_infeasible_cycles(model, fluxes, fix=()):
    """Remove thermodynamically infeasible cycles from a flux distribution.

    Arguments
    ---------
    model : cobra.Model
        The model that generated the flux distribution.
    fluxes : dict
        The flux distribution containing infeasible loops.

    Returns
    -------
    dict
        A cycle free flux distribution.

    References
    ----------
    .. [1]	A. A. Desouki, F. Jarre, G. Gelius-Dietrich, and M. J. Lercher, “CycleFreeFlux: efficient removal of
            thermodynamically infeasible loops from flux distributions.”
    """
    with model:
        # make sure the original object is restored
        exchange_reactions = model.boundary
        exchange_ids = [exchange.id for exchange in exchange_reactions]
        internal_reactions = [reaction for reaction in model.reactions if reaction.id not in exchange_ids]
        for exchange in exchange_reactions:
            exchange_flux = fluxes[exchange.id]
            exchange.bounds = (exchange_flux, exchange_flux)
        cycle_free_objective_list = []
        for internal_reaction in internal_reactions:
            internal_flux = fluxes[internal_reaction.id]
            if internal_flux >= 0:
                cycle_free_objective_list.append(Mul._from_args((FloatOne, internal_reaction.forward_variable)))
                internal_reaction.bounds = (0, internal_flux)
            else:  # internal_flux < 0:
                cycle_free_objective_list.append(Mul._from_args((FloatOne, internal_reaction.reverse_variable)))
                internal_reaction.bounds = (internal_flux, 0)
        cycle_free_objective = model.solver.interface.Objective(
            Add._from_args(cycle_free_objective_list), direction="min", sloppy=True
        )
        model.objective = cycle_free_objective

        for reaction_id in fix:
            reaction_to_fix = model.reactions.get_by_id(reaction_id)
            reaction_to_fix.bounds = (fluxes[reaction_id], fluxes[reaction_id])
        try:
            solution = model.optimize(raise_error=True)
        except OptimizationError as e:
            logger.warning("Couldn't remove cycles from reference flux distribution.")
            raise e
        result = solution.fluxes
        return result

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

    assert x.as_ordered_factors() == [x]
    assert (2 * x * x ** n * sin(x) * cos(x)).as_ordered_factors() == [Integer(2), x, x ** n, sin(x), cos(x)]

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

    assert expr.as_ordered_factors() == args

    A, B = symbols("A,B", commutative=False)

    assert (A * B).as_ordered_factors() == [A, B]
    assert (B * A).as_ordered_factors() == [B, A]

def linear_expand(expr):
    """
    If a sympy 'Expr' is of the form:

    expr = expr_0 + expr_1*a_1 + ... + expr_n*a_n

    where all the a_j are noncommuting symbols in basis then

    (expr_0, ..., expr_n) and (1, a_1, ..., a_n) are returned.  Note that
    expr_j*a_j does not have to be of that form, but rather can be any
    Mul with a_j as a factor (it doen not have to be a postmultiplier).
    expr_0 is the scalar part of the expression.
    """
    expr = expand(expr)
    if expr.is_commutative:  # commutative expr only contains expr_0
        return (expr, ), (S.One, )

    if isinstance(expr, Mul):  # expr only contains one term
        (coefs, bases) = expr.args_cnc()
        coefs = Mul(*coefs)
        bases = bases[0]
    elif isinstance(expr, Symbol):  # term is Symbol
        coefs = S.One
        bases = expr
    elif isinstance(expr, Add):  # expr has multiple terms
        coefs = []
        bases = []
        for arg in expr.args:
            term = arg.args_cnc()
            coef = Mul(*term[0])
            base = term[1][0]
            if base in bases:  # increment coefficient of base
                ibase = list(bases).index(base)  # Python 2.5
                coefs[ibase] += coef
            else:  # add base to list
                coefs.append(coef)
                bases.append(base)
    else:
        raise NotImplementedError("linear_expand for type %s" % type(expr))


    if not isinstance(coefs, list):  # convert single coef to list
        coefs = [coefs]
    if not isinstance(bases, list):  # convert single base to list
        bases = [bases]
    coefs = tuple(coefs)
    bases = tuple(bases)
    return coefs, bases

        def multiply(expr, mrow):
            from sympy.simplify import fraction
            numer, denom = fraction(expr)
            if denom is not S.One:
                frac = self.dom.createElement('mfrac')
                if self._settings["fold_short_frac"] and len(str(expr)) < 7:
                    frac.setAttribute('bevelled', 'true')
                xnum = self._print(numer)
                xden = self._print(denom)
                frac.appendChild(xnum)
                frac.appendChild(xden)
                mrow.appendChild(frac)
                return mrow

            coeff, terms = expr.as_coeff_mul()
            if coeff is S.One and len(terms) == 1:
                mrow.appendChild(self._print(terms[0]))
                return mrow
            if self.order != 'old':
                terms = Mul._from_args(terms).as_ordered_factors()

            if coeff != 1:
                x = self._print(coeff)
                y = self.dom.createElement('mo')
                y.appendChild(self.dom.createTextNode(self.mathml_tag(expr)))
                mrow.appendChild(x)
                mrow.appendChild(y)
            for term in terms:
                x = self._print(term)
                mrow.appendChild(x)
                if not term == terms[-1]:
                    y = self.dom.createElement('mo')
                    y.appendChild(self.dom.createTextNode(self.mathml_tag(expr)))
                    mrow.appendChild(y)
            return mrow

        def multiply(expr, mrow):
            from sympy.simplify import fraction
            numer, denom = fraction(expr)

            if denom is not S.One:
                frac = self.dom.createElement('mfrac')
                xnum = self._print(numer)
                xden = self._print(denom)
                frac.appendChild(xnum)
                frac.appendChild(xden)
                return frac

            coeff, terms = expr.as_coeff_mul()
            if coeff is S.One and len(terms) == 1:
                return self._print(terms[0])

            if self.order != 'old':
                terms = Mul._from_args(terms).as_ordered_factors()

            if(coeff != 1):
                x = self._print(coeff)
                y = self.dom.createElement('mo')
                y.appendChild(self.dom.createTextNode(self.mathml_tag(expr)))
                mrow.appendChild(x)
                mrow.appendChild(y)
            for term in terms:
                x = self._print(term)
                mrow.appendChild(x)
                if not term == terms[-1]:
                    y = self.dom.createElement('mo')
                    y.appendChild(self.dom.createTextNode(self.mathml_tag(expr)))
                    mrow.appendChild(y)
            return mrow

    def trig_replace(self, M, angle, name):
        """Replaces trigonometric expressions cos(x)
        and sin(x) by CX and SX

        Parameters
        ==========
        M: var or Matrix
            Object of substitution
        angle: var
            symbol that stands for the angle value
        name: int or string
            brief name X for the angle

        Notes
        =====
        The cos(x) and sin(x) will be replaced by CX and SX,
        where X is the name and x is the angle
        """
        if not isinstance(angle, Expr) or angle.is_number:
            return M
        cos_sym, sin_sym = tools.cos_sin_syms(name)
        sym_list = [(cos_sym, cos(angle)), (sin_sym, sin(angle))]
        subs_dict = {}
        for sym, sym_old in sym_list:
            if -1 in Mul.make_args(sym_old):
                sym_old = -sym_old
            subs_dict[sym_old] = sym
            self.add_to_dict(sym, sym_old)
        for i1 in xrange(M.shape[0]):
            for i2 in xrange(M.shape[1]):
                M[i1, i2] = M[i1, i2].subs(subs_dict)
        return M

def test_issue_3268():
    z = -5*sqrt(2)/(2*sqrt(2*sqrt(29) + 29)) + sqrt(-sqrt(29)/29 + S(1)/2)
    assert Mul(*[powsimp(a) for a in Mul.make_args(z.normal())]) == 0
    assert powsimp(z.normal()) == 0
    assert simplify(z) == 0
    assert powsimp(sqrt(2 + sqrt(3))*sqrt(2 - sqrt(3)) + 1) == 2
    assert powsimp(z) != 0

    def _print_Mul(self, expr):

        if _coeff_isneg(expr):
            x = self.dom.createElement('apply')
            x.appendChild(self.dom.createElement('minus'))
            x.appendChild(self._print_Mul(-expr))
            return x

        from sympy.simplify import fraction
        numer, denom = fraction(expr)

        if denom is not S.One:
            x = self.dom.createElement('apply')
            x.appendChild(self.dom.createElement('divide'))
            x.appendChild(self._print(numer))
            x.appendChild(self._print(denom))
            return x

        coeff, terms = expr.as_coeff_mul()
        if coeff is S.One and len(terms) == 1:
            # XXX since the negative coefficient has been handled, I don't
            # think a coeff of 1 can remain
            return self._print(terms[0])

        if self.order != 'old':
            terms = Mul._from_args(terms).as_ordered_factors()

        x = self.dom.createElement('apply')
        x.appendChild(self.dom.createElement('times'))
        if(coeff != 1):
            x.appendChild(self._print(coeff))
        for term in terms:
            x.appendChild(self._print(term))
        return x

 def _print_Mul(self, expr):
     if len(expr.args) >= 2:
         return "mul({}, {})".format(
             self._print(expr.args[0]),
             self._print(Mul.fromiter(expr.args[1:])),
         )
     else:
         return self._print(expr.args[0])

def test_zoo():
    b = Symbol('b', bounded=True)
    nz = Symbol('nz', nonzero=True)
    p = Symbol('p', positive=True)
    n = Symbol('n', negative=True)
    im = Symbol('i', imaginary=True)
    c = Symbol('c', complex=True)
    pb = Symbol('pb', positive=True, bounded=True)
    nb = Symbol('nb', negative=True, bounded=True)
    imb = Symbol('ib', imaginary=True, bounded=True)
    for i in [I, S.Infinity, S.NegativeInfinity, S.Zero, S.One, S.Pi, S.Half, S(3), log(3),
              b, nz, p, n, im, pb, nb, imb, c]:
        if i.is_bounded and (i.is_real or i.is_imaginary):
            assert i + zoo is zoo
            assert i - zoo is zoo
            assert zoo + i is zoo
            assert zoo - i is zoo
        elif i.is_bounded is not False:
            assert (i + zoo).is_Add
            assert (i - zoo).is_Add
            assert (zoo + i).is_Add
            assert (zoo - i).is_Add
        else:
            assert (i + zoo) is S.NaN
            assert (i - zoo) is S.NaN
            assert (zoo + i) is S.NaN
            assert (zoo - i) is S.NaN

        if i.is_nonzero and (i.is_real or i.is_imaginary):
            assert i*zoo is zoo
            assert zoo*i is zoo
        elif i.is_zero:
            assert i*zoo is S.NaN
            assert zoo*i is S.NaN
        else:
            assert (i*zoo).is_Mul
            assert (zoo*i).is_Mul

        if (1/i).is_nonzero and (i.is_real or i.is_imaginary):
            assert zoo/i is zoo
        elif (1/i).is_zero:
            assert zoo/i is S.NaN
        elif i.is_zero:
            assert zoo/i is zoo
        else:
            assert (zoo/i).is_Mul

    assert (I*oo).is_Mul  # allow directed infinity
    assert zoo + zoo is S.NaN
    assert zoo * zoo is zoo
    assert zoo - zoo is S.NaN
    assert zoo/zoo is S.NaN
    assert zoo**zoo is S.NaN
    assert zoo**0 is S.One
    assert zoo**2 is zoo
    assert 1/zoo is S.Zero

    assert Mul.flatten([S(-1), oo, S(0)]) == ([S.NaN], [], None)

 def flatten(cls, args):
     # TODO: disallow nested TensorProducts.
     c_part = []
     nc_parts = []
     for arg in args:
         cp, ncp = arg.args_cnc()
         c_part.extend(list(cp))
         nc_parts.append(Mul._from_args(ncp))
     return c_part, nc_parts

    def eval(cls, a, b):
        if not (a and b): return S.Zero
        if a == b: return Integer(2)*a**2
        if a.is_commutative or b.is_commutative:
            return Integer(2)*a*b

        # [xA,yB]  ->  xy*[A,B]
        # from sympy.physics.qmul import QMul
        ca, nca = a.args_cnc()
        cb, ncb = b.args_cnc()
        c_part = ca + cb
        if c_part:
            return Mul(Mul(*c_part), cls(Mul._from_args(nca), Mul._from_args(ncb)))

        # Canonical ordering of arguments
        #The Commutator [A,B] is on canonical form if A < B.
        if a.compare(b) == 1:
            return cls(b,a)

    def eval(cls, a, b):
        """The Commutator [A,B] is on canonical form if A < B.
        """
        if not (a and b): return S.Zero
        if a == b: return S.Zero
        if a.is_commutative or b.is_commutative:
            return S.Zero

        # [xA,yB]  ->  xy*[A,B]
        # from sympy.physics.qmul import QMul
        ca, nca = a.args_cnc()
        cb, ncb = b.args_cnc()
        c_part = list(ca) + list(cb)
        if c_part:
            return Mul(Mul(*c_part), cls(Mul._from_args(nca), Mul._from_args(ncb)))

        # Canonical ordering of arguments
        if a.compare(b) == 1:
            return S.NegativeOne*cls(b,a)