def predict_compound_pairs(reaction, compound_formula, solver, epsilon=<NUM_LIT>,<EOL>alt_elements=None, weight_func=default_weight):
elements = set()<EOL>for compound, value in reaction.compounds:<EOL><INDENT>if compound.name not in compound_formula:<EOL><INDENT>return<EOL><DEDENT>f = compound_formula[compound.name]<EOL>elements.update(e for e, _, _ in _weighted_formula(f, weight_func))<EOL><DEDENT>if len(elements) == <NUM_LIT:0>:<EOL><INDENT>return<EOL><DEDENT>p = solver.create_problem()<EOL>x = p.namespace(name='<STR_LIT:x>')<EOL>y = p.namespace(name='<STR_LIT:y>')<EOL>m = p.namespace(name='<STR_LIT:m>')<EOL>q = p.namespace(name='<STR_LIT:q>')<EOL>omega = p.namespace(name='<STR_LIT>')<EOL>gamma = p.namespace(name='<STR_LIT>')<EOL>delta = p.namespace(name='<STR_LIT>')<EOL>left, right = _reaction_to_dicts(reaction)<EOL>objective = <NUM_LIT:0><EOL>for c1, c2 in product(left, right):<EOL><INDENT>m.define([(c1, c2)], lower=<NUM_LIT:0>)<EOL>q.define([(c1, c2)], lower=<NUM_LIT:0>)<EOL>omega.define([(c1, c2)], types=lp.VariableType.Binary)<EOL>objective += <NUM_LIT:100> * m[c1, c2] + <NUM_LIT> * q[c1, c2] - <NUM_LIT:0.1> * omega[c1, c2]<EOL>gamma.define(((c1, c2, e) for e in elements),<EOL>types=lp.VariableType.Binary)<EOL>delta.define(((c1, c2, e) for e in elements), lower=<NUM_LIT:0>)<EOL>x.define([(c1, c2)], types=lp.VariableType.Binary)<EOL>y.define([(c1, c2)], types=lp.VariableType.Binary)<EOL>p.add_linear_constraints(y[c1, c2] <= <NUM_LIT:1> - x[c1, c2])<EOL>delta_wsum = delta.expr(<EOL>((c1, c2, e), weight_func(e)) for e in elements)<EOL>p.add_linear_constraints(<EOL>m[c1, c2] <= delta_wsum, q[c1, c2] <= delta_wsum)<EOL>objective -= gamma.sum((c1, c2, e) for e in elements)<EOL><DEDENT>p.set_objective(objective)<EOL>zs = {}<EOL>for c1, v in iteritems(left):<EOL><INDENT>f = compound_formula[c1.name]<EOL>for e in elements:<EOL><INDENT>mf = f.get(e, <NUM_LIT:0>)<EOL>delta_sum = delta.sum((c1, c2, e) for c2 in right)<EOL>try:<EOL><INDENT>p.add_linear_constraints(delta_sum == v * mf)<EOL><DEDENT>except ValueError:<EOL><INDENT>raise UnbalancedReactionError('<STR_LIT>')<EOL><DEDENT><DEDENT>x_sum = x.sum((c1, c2) for c2 in right)<EOL>y_sum = y.sum((c1, c2) for c2 in right)<EOL>p.add_linear_constraints(x_sum <= <NUM_LIT:1>, y_sum <= <NUM_LIT:1>)<EOL>zs[c1] = <NUM_LIT:0><EOL>for e, mf, w in _weighted_formula(f, weight_func):<EOL><INDENT>zs[c1] += w * mf<EOL><DEDENT><DEDENT>for c2, v in iteritems(right):<EOL><INDENT>f = compound_formula[c2.name]<EOL>for e in elements:<EOL><INDENT>mf = f.get(e, <NUM_LIT:0>)<EOL>delta_sum = delta.sum((c1, c2, e) for c1 in left)<EOL>try:<EOL><INDENT>p.add_linear_constraints(delta_sum == v * mf)<EOL><DEDENT>except ValueError:<EOL><INDENT>raise UnbalancedReactionError('<STR_LIT>')<EOL><DEDENT><DEDENT><DEDENT>for c1, v1 in iteritems(left):<EOL><INDENT>for c2 in right:<EOL><INDENT>f1 = compound_formula[c1.name]<EOL>f2 = compound_formula[c2.name]<EOL>for e in elements:<EOL><INDENT>mf = f1.get(e, <NUM_LIT:0>)<EOL>p.add_linear_constraints(<EOL>delta[c1, c2, e] <= float(v1) * mf * omega[c1, c2])<EOL>if e in f2:<EOL><INDENT>p.add_linear_constraints(<EOL>delta[c1, c2, e] <= float(v1) * mf * gamma[c1, c2, e])<EOL><DEDENT><DEDENT>p.add_linear_constraints(<EOL>m[c1, c2] <= float(v1) * zs[c1] * x[c1, c2],<EOL>q[c1, c2] <= float(v1) * zs[c1] * y[c1, c2])<EOL><DEDENT><DEDENT>try:<EOL><INDENT>result = p.solve(lp.ObjectiveSense.Maximize)<EOL><DEDENT>except lp.SolverError:<EOL><INDENT>raise UnbalancedReactionError('<STR_LIT>')<EOL><DEDENT>first_objective = result.get_value(objective)<EOL>yield dict(_transfers(reaction, delta, elements, result, epsilon))<EOL>if alt_elements is None:<EOL><INDENT>alt_elements = elements<EOL><DEDENT>else:<EOL><INDENT>alt_elements = elements.intersection(alt_elements)<EOL><DEDENT>if len(alt_elements) == <NUM_LIT:0>:<EOL><INDENT>return<EOL><DEDENT>add_objective_constraint = True<EOL>while True:<EOL><INDENT>elem_vars = <NUM_LIT:0><EOL>elem_count = <NUM_LIT:0><EOL>for c1, c2 in product(left, right):<EOL><INDENT>for e in alt_elements:<EOL><INDENT>if result.get_value(gamma[c1, c2, e]) > <NUM_LIT:0.5>:<EOL><INDENT>elem_count += <NUM_LIT:1><EOL>elem_vars += gamma[c1, c2, e]<EOL><DEDENT><DEDENT><DEDENT>if elem_count == <NUM_LIT:0>:<EOL><INDENT>break<EOL><DEDENT>if add_objective_constraint:<EOL><INDENT>p.add_linear_constraints(objective >= first_objective)<EOL>add_objective_constraint = False<EOL><DEDENT>p.add_linear_constraints(elem_vars <= elem_count - <NUM_LIT:1>)<EOL>try:<EOL><INDENT>result = p.solve(lp.ObjectiveSense.Maximize)<EOL><DEDENT>except lp.SolverError:<EOL><INDENT>break<EOL><DEDENT>yield dict(_transfers(reaction, delta, elements, result, epsilon))<EOL><DEDENT>
Predict compound pairs of reaction using MapMaker.

    Yields all solutions as dictionaries with compound pairs as keys and
    formula objects as values.

    Args:
        reaction: :class:`psamm.reaction.Reaction` object.
        compound_formula: Dictionary mapping compound IDs to formulas. Formulas
            must be flattened.
        solver: LP solver (MILP).
        epsilon: Threshold for rounding floating-point values to integers in
            the predicted transfers.
        alt_elements: Iterable of elements to consider for alternative
              solutions. Only alternate solutions that have different transfers
              of these elements will be returned (default=all elements).
        weight_func: Optional function that returns a weight for a formula
            element (should handle specific Atom and Radical objects). By
            default, the standard MapMaker weights will be used
            (H=0, R=40, N=0.4, O=0.4, P=0.4, *=1).