docs: def check_compound_consistency(database, solver, exchange=set(),
                                 zeromass=set()):
      """Yield each compound in the database with assigned mass
  
      Each compound will be assigned a mass and the number of compounds having a
      positive mass will be approximately maximized.
  
      This is an implementation of the solution originally proposed by
      [Gevorgyan08]_  but using the new method proposed by [Thiele14]_ to avoid
      MILP constraints. This is similar to the way Fastcore avoids MILP
      contraints.
      """
  
      # Create mass balance problem
      prob = solver.create_problem()
      compound_set = _non_localized_compounds(database)
      mass_compounds = compound_set.difference(zeromass)
  
      # Define mass variables
      m = prob.namespace(mass_compounds, lower=0)
  
      # Define z variables
      z = prob.namespace(mass_compounds, lower=0, upper=1)
      prob.set_objective(z.sum(mass_compounds))
  
      prob.add_linear_constraints(m.set(mass_compounds) >= z.set(mass_compounds))
  
      massbalance_lhs = {reaction_id: 0 for reaction_id in database.reactions}
      for (compound, reaction_id), value in iteritems(database.matrix):
          if compound not in zeromass:
              mass_var = m(compound.in_compartment(None))
              massbalance_lhs[reaction_id] += mass_var * value
      for reaction_id, lhs in iteritems(massbalance_lhs):
          if reaction_id not in exchange:
              prob.add_linear_constraints(lhs == 0)
  
      # Solve
      try:
          prob.solve(lp.ObjectiveSense.Maximize)
      except lp.SolverError as e:
          raise_from(
              MassConsistencyError('Failed to solve mass consistency: {}'.format(
                  e)), e)
  
      for compound in mass_compounds:
          yield compound, m.value(compound)