Yield each compound in the database with assigned mass

    Each compound will be assigned a mass and the number of compounds having a
    positive mass will be approximately maximized.

    This is an implementation of the solution originally proposed by
    [Gevorgyan08]_  but using the new method proposed by [Thiele14]_ to avoid
    MILP constraints. This is similar to the way Fastcore avoids MILP
    contraints.

def check_compound_consistency(database, solver, exchange=set(),
                               zeromass=set()):
    """Yield each compound in the database with assigned mass

    Each compound will be assigned a mass and the number of compounds having a
    positive mass will be approximately maximized.

    This is an implementation of the solution originally proposed by
    [Gevorgyan08]_  but using the new method proposed by [Thiele14]_ to avoid
    MILP constraints. This is similar to the way Fastcore avoids MILP
    contraints.
    """

    # Create mass 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=0)

    # Define z variables
    z = prob.namespace(mass_compounds, lower=0, upper=1)
    prob.set_objective(z.sum(mass_compounds))

    prob.add_linear_constraints(m.set(mass_compounds) >= z.set(mass_compounds))

    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:
            prob.add_linear_constraints(lhs == 0)

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

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