docs: def check_reaction_consistency(database, solver, exchange=set(),
                                 checked=set(), zeromass=set(), weights={}):
      """Check inconsistent reactions by minimizing mass residuals
  
      Return a reaction iterable, and compound iterable. The reaction iterable
      yields reaction ids and mass residuals. The compound iterable yields
      compound ids and mass assignments.
  
      Each compound is assigned a mass of at least one, and the masses are
      balanced using the stoichiometric matrix. In addition, each reaction has a
      residual mass that is included in the mass balance equations. The L1-norm
      of the residuals is minimized. Reactions in the checked set are assumed to
      have been manually checked and therefore have the residual fixed at zero.
      """
  
      # Create Flux 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=1)
  
      # Define residual mass variables and objective constriants
      z = prob.namespace(database.reactions, lower=0)
      r = prob.namespace(database.reactions)
  
      objective = z.expr((reaction_id, weights.get(reaction_id, 1))
                         for reaction_id in database.reactions)
      prob.set_objective(objective)
  
      rs = r.set(database.reactions)
      zs = z.set(database.reactions)
      prob.add_linear_constraints(zs >= rs, rs >= -zs)
  
      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:
              if reaction_id not in checked:
                  prob.add_linear_constraints(lhs + r(reaction_id) == 0)
              else:
                  prob.add_linear_constraints(lhs == 0)
  
      # Solve
      try:
          prob.solve(lp.ObjectiveSense.Minimize)
      except lp.SolverError as e:
          raise_from(
              MassConsistencyError('Failed to solve mass consistency: {}'.format(
                  e)), e)
  
      def iterate_reactions():
          for reaction_id in database.reactions:
              residual = r.value(reaction_id)
              yield reaction_id, residual
  
      def iterate_compounds():
          for compound in mass_compounds:
              yield compound, m.value(compound)
  
      return iterate_reactions(), iterate_compounds()