docs: def predict_compound_pairs_iterated(
          reactions, formulas, prior=(1, 43), max_iterations=None,
          element_weight=element_weight):
      """Predict reaction pairs using iterated method.
  
      Returns a tuple containing a dictionary of predictions keyed by the
      reaction IDs, and the final number of iterations. Each reaction prediction
      entry contains a tuple with a dictionary of transfers and a dictionary of
      unbalanced compounds. The dictionary of unbalanced compounds is empty only
      if the reaction is balanced.
  
      Args:
          reactions: Dictionary or pair-iterable of (id, equation) pairs.
              IDs must be any hashable reaction identifier (e.g. string) and
              equation must be :class:`psamm.reaction.Reaction` objects.
          formulas: Dictionary mapping compound IDs to
              :class:`psamm.formula.Formula`. Formulas must be flattened.
          prior: Tuple of (alpha, beta) parameters for the MAP inference.
              If not provided, the default parameters will be used: (1, 43).
          max_iterations: Maximum iterations to run before stopping. If the
              stopping condition is reached before this number of iterations,
              the procedure also stops. If None, the procedure only stops when
              the stopping condition is reached.
          element_weight: A function providing returning weight value for the
              given :class:`psamm.formula.Atom` or
              :class:`psamm.formula.Radical`. If not provided, the default weight
              will be used (H=0, C=1, *=0.82)
      """
      prior_alpha, prior_beta = prior
      reactions = dict(reactions)
  
      pair_reactions = {}
      possible_pairs = Counter()
      for reaction_id, equation in iteritems(reactions):
          for (c1, _), (c2, _) in product(equation.left, equation.right):
              spair = tuple(sorted([c1.name, c2.name]))
              possible_pairs[spair] += 1
              pair_reactions.setdefault(spair, set()).add(reaction_id)
  
      next_reactions = set(reactions)
      pairs_predicted = None
      prediction = {}
      weights = {}
      iteration = 0
      while len(next_reactions) > 0:
          iteration += 1
          if max_iterations is not None and iteration > max_iterations:
              break
  
          logger.info('Iteration {}: {} reactions...'.format(
              iteration, len(next_reactions)))
  
          for reaction_id in next_reactions:
              result = predict_compound_pairs(
                  reactions[reaction_id], formulas, weights, element_weight)
              if result is None:
                  continue
  
              transfer, balance = result
  
              rpairs = {}
              for ((c1, _), (c2, _)), form in iteritems(transfer):
                  rpairs.setdefault((c1, c2), []).append(form)
  
              prediction[reaction_id] = rpairs, balance
  
          pairs_predicted = Counter()
          for reaction_id, (rpairs, _) in iteritems(prediction):
              for c1, c2 in rpairs:
                  spair = tuple(sorted([c1.name, c2.name]))
                  pairs_predicted[spair] += 1
  
          next_reactions = set()
          for spair, total in sorted(iteritems(possible_pairs)):
              pred = pairs_predicted[spair]
  
              # The weight is set to the maximum a posteriori (MAP) estimate
              # of the primary pair probability distribution.
              posterior_alpha = prior_alpha + pred
              posterior_beta = prior_beta + total - pred
              pair_weight = ((posterior_alpha - 1) /
                             (posterior_alpha + posterior_beta - 2))
  
              if (spair not in weights or
                      abs(pair_weight - weights[spair]) > 1e-5):
                  next_reactions.update(pair_reactions[spair])
  
              c1, c2 = spair
              weights[c1, c2] = pair_weight
              weights[c2, c1] = pair_weight
  
      return prediction, iteration