Check inconsistent reactions by minimizing mass residuals

    Return a reaction iterable, and compound iterable. The reaction iterable
    yields reaction ids and mass residuals. The compound iterable yields
    compound ids and mass assignments.

    Each compound is assigned a mass of at least one, and the masses are
    balanced using the stoichiometric matrix. In addition, each reaction has a
    residual mass that is included in the mass balance equations. The L1-norm
    of the residuals is minimized. Reactions in the checked set are assumed to
    have been manually checked and therefore have the residual fixed at zero.

def check_reaction_consistency(database, solver, exchange=set(),
                               checked=set(), zeromass=set(), weights={}):
    """Check inconsistent reactions by minimizing mass residuals

    Return a reaction iterable, and compound iterable. The reaction iterable
    yields reaction ids and mass residuals. The compound iterable yields
    compound ids and mass assignments.

    Each compound is assigned a mass of at least one, and the masses are
    balanced using the stoichiometric matrix. In addition, each reaction has a
    residual mass that is included in the mass balance equations. The L1-norm
    of the residuals is minimized. Reactions in the checked set are assumed to
    have been manually checked and therefore have the residual fixed at zero.
    """

    # Create Flux balance problem
    prob = solver.create_problem()
    compound_set = _non_localized_compounds(database)
    mass_compounds = compound_set.difference(zeromass)

    # Define mass variables
    m = prob.namespace(mass_compounds, lower=1)

    # Define residual mass variables and objective constriants
    z = prob.namespace(database.reactions, lower=0)
    r = prob.namespace(database.reactions)

    objective = z.expr((reaction_id, weights.get(reaction_id, 1))
                       for reaction_id in database.reactions)
    prob.set_objective(objective)

    rs = r.set(database.reactions)
    zs = z.set(database.reactions)
    prob.add_linear_constraints(zs >= rs, rs >= -zs)

    massbalance_lhs = {reaction_id: 0 for reaction_id in database.reactions}
    for (compound, reaction_id), value in iteritems(database.matrix):
        if compound not in zeromass:
            mass_var = m(compound.in_compartment(None))
            massbalance_lhs[reaction_id] += mass_var * value
    for reaction_id, lhs in iteritems(massbalance_lhs):
        if reaction_id not in exchange:
            if reaction_id not in checked:
                prob.add_linear_constraints(lhs + r(reaction_id) == 0)
            else:
                prob.add_linear_constraints(lhs == 0)

    # Solve
    try:
        prob.solve(lp.ObjectiveSense.Minimize)
    except lp.SolverError as e:
        raise_from(
            MassConsistencyError('Failed to solve mass consistency: {}'.format(
                e)), e)

    def iterate_reactions():
        for reaction_id in database.reactions:
            residual = r.value(reaction_id)
            yield reaction_id, residual

    def iterate_compounds():
        for compound in mass_compounds:
            yield compound, m.value(compound)

    return iterate_reactions(), iterate_compounds()