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