docs: 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()