Apply thermodynamic constraints to the model.

        Adding these constraints restricts the solution space to only
        contain solutions that have no internal loops [Schilling00]_. This is
        solved as a MILP problem as described in [Muller13]_. The time to solve
        a problem with thermodynamic constraints is usually much longer than a
        normal FBA problem.

        The ``em`` parameter is the upper bound on the delta mu reaction
        variables. This parameter has to be balanced based on the model size
        since setting the value too low can result in the correct solutions
        being infeasible and setting the value too high can result in
        numerical instability which again makes the correct solutions
        infeasible. The default value should work in all cases as long as the
        model is not unusually large.

def add_thermodynamic(self, em=1000):
        """Apply thermodynamic constraints to the model.

        Adding these constraints restricts the solution space to only
        contain solutions that have no internal loops [Schilling00]_. This is
        solved as a MILP problem as described in [Muller13]_. The time to solve
        a problem with thermodynamic constraints is usually much longer than a
        normal FBA problem.

        The ``em`` parameter is the upper bound on the delta mu reaction
        variables. This parameter has to be balanced based on the model size
        since setting the value too low can result in the correct solutions
        being infeasible and setting the value too high can result in
        numerical instability which again makes the correct solutions
        infeasible. The default value should work in all cases as long as the
        model is not unusually large.
        """

        internal = set(r for r in self._model.reactions
                       if not self._model.is_exchange(r))

        # Reaction fluxes
        v = self._v

        # Indicator variable
        alpha = self._prob.namespace(internal, types=lp.VariableType.Binary)

        # Delta mu is the stoichiometrically weighted sum of the compound mus.
        dmu = self._prob.namespace(internal)

        for reaction_id in self._model.reactions:
            if not self._model.is_exchange(reaction_id):
                flux = v(reaction_id)
                alpha_r = alpha(reaction_id)
                dmu_r = dmu(reaction_id)

                lower, upper = self._model.limits[reaction_id]

                # Constrain the reaction to a direction determined by alpha
                # and contrain the delta mu to a value in [-em; -1] if
                # alpha is one, otherwise in [1; em].
                self._prob.add_linear_constraints(
                    flux >= lower * (1 - alpha_r),
                    flux <= upper * alpha_r,
                    dmu_r >= -em * alpha_r + (1 - alpha_r),
                    dmu_r <= em * (1 - alpha_r) - alpha_r)

        # Define mu variables
        mu = self._prob.namespace(self._model.compounds)

        tdbalance_lhs = {reaction_id: 0
                         for reaction_id in self._model.reactions}
        for spec, value in iteritems(self._model.matrix):
            compound, reaction_id = spec
            if not self._model.is_exchange(reaction_id):
                tdbalance_lhs[reaction_id] += mu(compound) * value
        for reaction_id, lhs in iteritems(tdbalance_lhs):
            if not self._model.is_exchange(reaction_id):
                self._prob.add_linear_constraints(lhs == dmu(reaction_id))