def gapfill(<EOL>model, core, blocked, exclude, solver, epsilon=<NUM_LIT>, v_max=<NUM_LIT:1000>,<EOL>weights={}, implicit_sinks=True, allow_bounds_expansion=False):
prob = solver.create_problem()<EOL>min_tol = prob.integrality_tolerance.min<EOL>int_tol = _find_integer_tolerance(epsilon, v_max, min_tol)<EOL>if int_tol < prob.integrality_tolerance.value:<EOL><INDENT>prob.integrality_tolerance.value = int_tol<EOL><DEDENT>v = prob.namespace(model.reactions, lower=-v_max, upper=v_max)<EOL>database_reactions = set(model.reactions).difference(core, exclude)<EOL>ym = prob.namespace(model.reactions, types=lp.VariableType.Binary)<EOL>yd = prob.namespace(database_reactions, types=lp.VariableType.Binary)<EOL>objective = ym.expr(<EOL>(rxnid, weights.get(rxnid, <NUM_LIT:1>)) for rxnid in model.reactions)<EOL>objective += yd.expr(<EOL>(rxnid, weights.get(rxnid, <NUM_LIT:1>)) for rxnid in database_reactions)<EOL>prob.set_objective(objective)<EOL>for reaction_id in model.reactions:<EOL><INDENT>lower, upper = (float(x) for x in model.limits[reaction_id])<EOL>if reaction_id in exclude or not allow_bounds_expansion:<EOL><INDENT>prob.add_linear_constraints(<EOL>upper >= v(reaction_id), v(reaction_id) >= lower)<EOL><DEDENT>else:<EOL><INDENT>delta_lower = min(<NUM_LIT:0>, -v_max - lower)<EOL>delta_upper = max(<NUM_LIT:0>, v_max - upper)<EOL>prob.add_linear_constraints(<EOL>v(reaction_id) >= lower + ym(reaction_id) * delta_lower,<EOL>v(reaction_id) <= upper + ym(reaction_id) * delta_upper)<EOL><DEDENT><DEDENT>for reaction_id in database_reactions:<EOL><INDENT>lower, upper = model.limits[reaction_id]<EOL>prob.add_linear_constraints(<EOL>v(reaction_id) >= yd(reaction_id) * -v_max,<EOL>v(reaction_id) <= yd(reaction_id) * v_max)<EOL><DEDENT>w = prob.namespace(types=lp.VariableType.Binary)<EOL>binary_cons_lhs = {compound: <NUM_LIT:0> for compound in blocked}<EOL>for (compound, reaction_id), value in iteritems(model.matrix):<EOL><INDENT>if reaction_id not in exclude and compound in blocked and value != <NUM_LIT:0>:<EOL><INDENT>w.define([(compound, reaction_id)])<EOL>w_var = w((compound, reaction_id))<EOL>dv = v(reaction_id) if value > <NUM_LIT:0> else -v(reaction_id)<EOL>prob.add_linear_constraints(<EOL>dv <= v_max * w_var,<EOL>dv >= epsilon + (-v_max - epsilon) * (<NUM_LIT:1> - w_var))<EOL>binary_cons_lhs[compound] += w_var<EOL><DEDENT><DEDENT>for compound, lhs in iteritems(binary_cons_lhs):<EOL><INDENT>prob.add_linear_constraints(lhs >= <NUM_LIT:1>)<EOL><DEDENT>massbalance_lhs = {compound: <NUM_LIT:0> for compound in model.compounds}<EOL>for (compound, reaction_id), value in iteritems(model.matrix):<EOL><INDENT>if reaction_id not in exclude:<EOL><INDENT>massbalance_lhs[compound] += v(reaction_id) * value<EOL><DEDENT><DEDENT>for compound, lhs in iteritems(massbalance_lhs):<EOL><INDENT>if implicit_sinks:<EOL><INDENT>prob.add_linear_constraints(lhs >= <NUM_LIT:0>)<EOL><DEDENT>else:<EOL><INDENT>prob.add_linear_constraints(lhs == <NUM_LIT:0>)<EOL><DEDENT><DEDENT>try:<EOL><INDENT>prob.solve(lp.ObjectiveSense.Minimize)<EOL><DEDENT>except lp.SolverError as e:<EOL><INDENT>raise_from(GapFillError('<STR_LIT>'.format(e)), e)<EOL><DEDENT>def added_iter():<EOL><INDENT>for reaction_id in database_reactions:<EOL><INDENT>if yd.value(reaction_id) > <NUM_LIT:0.5>:<EOL><INDENT>yield reaction_id<EOL><DEDENT><DEDENT><DEDENT>def no_bounds_iter():<EOL><INDENT>for reaction_id in model.reactions:<EOL><INDENT>if ym.value(reaction_id) > <NUM_LIT:0.5>:<EOL><INDENT>yield reaction_id<EOL><DEDENT><DEDENT><DEDENT>return added_iter(), no_bounds_iter()<EOL>
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.