Identify compounds in the model that cannot be produced.

    Yields all compounds that cannot be produced. This method
    assumes implicit sinks for all compounds in the model so
    the only factor that influences whether a compound can be
    produced is the presence of the compounds needed to produce it.

    Epsilon indicates the threshold amount of reaction flux for the products
    to be considered non-blocked. V_max indicates the maximum flux.

    This method is implemented as a MILP-program. Therefore it may
    not be efficient for larger models.

    Args:
        model: :class:`MetabolicModel` containing core reactions and reactions
            that can be added for gap-filling.
        solver: MILP solver instance.
        epsilon: Threshold amount of a compound produced for it to not be
            considered blocked.
        v_max: Maximum flux.
        implicit_sinks: Whether implicit sinks for all compounds are included
            when gap-filling (traditional GapFill uses implicit sinks).

def gapfind(model, solver, epsilon=0.001, v_max=1000, implicit_sinks=True):
    """Identify compounds in the model that cannot be produced.

    Yields all compounds that cannot be produced. This method
    assumes implicit sinks for all compounds in the model so
    the only factor that influences whether a compound can be
    produced is the presence of the compounds needed to produce it.

    Epsilon indicates the threshold amount of reaction flux for the products
    to be considered non-blocked. V_max indicates the maximum flux.

    This method is implemented as a MILP-program. Therefore it may
    not be efficient for larger models.

    Args:
        model: :class:`MetabolicModel` containing core reactions and reactions
            that can be added for gap-filling.
        solver: MILP solver instance.
        epsilon: Threshold amount of a compound produced for it to not be
            considered blocked.
        v_max: Maximum flux.
        implicit_sinks: Whether implicit sinks for all compounds are included
            when gap-filling (traditional GapFill uses implicit sinks).
    """
    prob = solver.create_problem()

    # Set integrality tolerance such that w constraints are correct
    min_tol = prob.integrality_tolerance.min
    int_tol = _find_integer_tolerance(epsilon, v_max, min_tol)
    if int_tol < prob.integrality_tolerance.value:
        prob.integrality_tolerance.value = int_tol

    # Define flux variables
    v = prob.namespace()
    for reaction_id in model.reactions:
        lower, upper = model.limits[reaction_id]
        v.define([reaction_id], lower=lower, upper=upper)

    # Define constraints on production of metabolites in reaction
    w = prob.namespace(types=lp.VariableType.Binary)
    binary_cons_lhs = {compound: 0 for compound in model.compounds}
    for spec, value in iteritems(model.matrix):
        compound, reaction_id = spec
        if value != 0:
            w.define([spec])
            w_var = w(spec)

            lower, upper = (float(x) for x in model.limits[reaction_id])
            if value > 0:
                dv = v(reaction_id)
            else:
                dv = -v(reaction_id)
                lower, upper = -upper, -lower

            prob.add_linear_constraints(
                dv <= upper * w_var,
                dv >= epsilon + (lower - epsilon) * (1 - w_var))

            binary_cons_lhs[compound] += w_var

    xp = prob.namespace(model.compounds, types=lp.VariableType.Binary)
    objective = xp.sum(model.compounds)
    prob.set_objective(objective)

    for compound, lhs in iteritems(binary_cons_lhs):
        prob.add_linear_constraints(lhs >= xp(compound))

    # Define mass balance constraints
    massbalance_lhs = {compound: 0 for compound in model.compounds}
    for spec, value in iteritems(model.matrix):
        compound, reaction_id = spec
        massbalance_lhs[compound] += v(reaction_id) * value
    for compound, lhs in iteritems(massbalance_lhs):
        if implicit_sinks:
            # The constraint is merely >0 meaning that we have implicit sinks
            # for all compounds.
            prob.add_linear_constraints(lhs >= 0)
        else:
            prob.add_linear_constraints(lhs == 0)

    # Solve
    try:
        result = prob.solve(lp.ObjectiveSense.Maximize)
    except lp.SolverError as e:
        raise_from(GapFillError('Failed to solve gapfill: {}'.format(e), e))

    for compound in model.compounds:
        if result.get_value(xp(compound)) < 0.5:
            yield compound