Find a set of reactions to add such that no compounds are blocked.

    Returns two iterators: first an iterator of reactions not in
    core, that were added to resolve the model. Second, an
    iterator of reactions in core that had flux bounds expanded (i.e.
    irreversible reactions become reversible). Similarly to
    GapFind, this method assumes, by default, 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. This means that the resulting model will not necessarily be
    flux consistent.

    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.
        core: The set of core (already present) reactions in the model.
        blocked: The compounds to unblock.
        exclude: Set of reactions in core to be excluded from gap-filling (e.g.
            biomass reaction).
        solver: MILP solver instance.
        epsilon: Threshold amount of a compound produced for it to not be
            considered blocked.
        v_max: Maximum flux.
        weights: Dictionary of weights for reactions. Weight is the penalty
            score for adding the reaction (non-core reactions) or expanding the
            flux bounds (all reactions).
        implicit_sinks: Whether implicit sinks for all compounds are included
            when gap-filling (traditional GapFill uses implicit sinks).
        allow_bounds_expansion: Allow flux bounds to be expanded at the cost
            of a penalty which can be specified using weights (traditional
            GapFill does not allow this). This includes turning irreversible
            reactions reversible.

def gapfill(
        model, core, blocked, exclude, solver, epsilon=0.001, v_max=1000,
        weights={}, implicit_sinks=True, allow_bounds_expansion=False):
    """Find a set of reactions to add such that no compounds are blocked.

    Returns two iterators: first an iterator of reactions not in
    core, that were added to resolve the model. Second, an
    iterator of reactions in core that had flux bounds expanded (i.e.
    irreversible reactions become reversible). Similarly to
    GapFind, this method assumes, by default, 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. This means that the resulting model will not necessarily be
    flux consistent.

    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.
        core: The set of core (already present) reactions in the model.
        blocked: The compounds to unblock.
        exclude: Set of reactions in core to be excluded from gap-filling (e.g.
            biomass reaction).
        solver: MILP solver instance.
        epsilon: Threshold amount of a compound produced for it to not be
            considered blocked.
        v_max: Maximum flux.
        weights: Dictionary of weights for reactions. Weight is the penalty
            score for adding the reaction (non-core reactions) or expanding the
            flux bounds (all reactions).
        implicit_sinks: Whether implicit sinks for all compounds are included
            when gap-filling (traditional GapFill uses implicit sinks).
        allow_bounds_expansion: Allow flux bounds to be expanded at the cost
            of a penalty which can be specified using weights (traditional
            GapFill does not allow this). This includes turning irreversible
            reactions reversible.
    """
    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(model.reactions, lower=-v_max, upper=v_max)

    # Add binary indicator variables
    database_reactions = set(model.reactions).difference(core, exclude)
    ym = prob.namespace(model.reactions, types=lp.VariableType.Binary)
    yd = prob.namespace(database_reactions, types=lp.VariableType.Binary)

    objective = ym.expr(
        (rxnid, weights.get(rxnid, 1)) for rxnid in model.reactions)
    objective += yd.expr(
        (rxnid, weights.get(rxnid, 1)) for rxnid in database_reactions)
    prob.set_objective(objective)

    # Add constraints on all reactions
    for reaction_id in model.reactions:
        lower, upper = (float(x) for x in model.limits[reaction_id])

        if reaction_id in exclude or not allow_bounds_expansion:
            prob.add_linear_constraints(
                upper >= v(reaction_id), v(reaction_id) >= lower)
        else:
            # Allow flux bounds to expand up to v_max with penalty
            delta_lower = min(0, -v_max - lower)
            delta_upper = max(0, v_max - upper)
            prob.add_linear_constraints(
                v(reaction_id) >= lower + ym(reaction_id) * delta_lower,
                v(reaction_id) <= upper + ym(reaction_id) * delta_upper)

    # Add constraints on database reactions
    for reaction_id in database_reactions:
        lower, upper = model.limits[reaction_id]
        prob.add_linear_constraints(
            v(reaction_id) >= yd(reaction_id) * -v_max,
            v(reaction_id) <= yd(reaction_id) * v_max)

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

            dv = v(reaction_id) if value > 0 else -v(reaction_id)
            prob.add_linear_constraints(
                dv <= v_max * w_var,
                dv >= epsilon + (-v_max - epsilon) * (1 - w_var))

            binary_cons_lhs[compound] += w_var

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

    # Define mass balance constraints
    massbalance_lhs = {compound: 0 for compound in model.compounds}
    for (compound, reaction_id), value in iteritems(model.matrix):
        if reaction_id not in exclude:
            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:
        prob.solve(lp.ObjectiveSense.Minimize)
    except lp.SolverError as e:
        raise_from(GapFillError('Failed to solve gapfill: {}'.format(e)), e)

    def added_iter():
        for reaction_id in database_reactions:
            if yd.value(reaction_id) > 0.5:
                yield reaction_id

    def no_bounds_iter():
        for reaction_id in model.reactions:
            if ym.value(reaction_id) > 0.5:
                yield reaction_id

    return added_iter(), no_bounds_iter()