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).

def predict_compound_pairs(reaction, compound_formula, solver, epsilon=1e-5,
                           alt_elements=None, weight_func=default_weight):
    """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).
    """
    elements = set()
    for compound, value in reaction.compounds:
        if compound.name not in compound_formula:
            return

        f = compound_formula[compound.name]
        elements.update(e for e, _, _ in _weighted_formula(f, weight_func))

    if len(elements) == 0:
        return

    p = solver.create_problem()

    x = p.namespace(name='x')
    y = p.namespace(name='y')
    m = p.namespace(name='m')
    q = p.namespace(name='q')
    omega = p.namespace(name='omega')
    gamma = p.namespace(name='gamma')
    delta = p.namespace(name='delta')

    left, right = _reaction_to_dicts(reaction)

    objective = 0
    for c1, c2 in product(left, right):
        m.define([(c1, c2)], lower=0)
        q.define([(c1, c2)], lower=0)
        omega.define([(c1, c2)], types=lp.VariableType.Binary)

        objective += 100 * m[c1, c2] + 90 * q[c1, c2] - 0.1 * omega[c1, c2]

        gamma.define(((c1, c2, e) for e in elements),
                     types=lp.VariableType.Binary)
        delta.define(((c1, c2, e) for e in elements), lower=0)

        # Eq 12
        x.define([(c1, c2)], types=lp.VariableType.Binary)
        y.define([(c1, c2)], types=lp.VariableType.Binary)
        p.add_linear_constraints(y[c1, c2] <= 1 - x[c1, c2])

        # Eq 6, 9
        delta_wsum = delta.expr(
            ((c1, c2, e), weight_func(e)) for e in elements)
        p.add_linear_constraints(
            m[c1, c2] <= delta_wsum, q[c1, c2] <= delta_wsum)

        objective -= gamma.sum((c1, c2, e) for e in elements)

    p.set_objective(objective)

    # Eq 3
    zs = {}
    for c1, v in iteritems(left):
        f = compound_formula[c1.name]
        for e in elements:
            mf = f.get(e, 0)
            delta_sum = delta.sum((c1, c2, e) for c2 in right)
            try:
                p.add_linear_constraints(delta_sum == v * mf)
            except ValueError:
                raise UnbalancedReactionError('Unable to add constraint')

        # Eq 8, 11
        x_sum = x.sum((c1, c2) for c2 in right)
        y_sum = y.sum((c1, c2) for c2 in right)
        p.add_linear_constraints(x_sum <= 1, y_sum <= 1)

        # Eq 13
        zs[c1] = 0
        for e, mf, w in _weighted_formula(f, weight_func):
            zs[c1] += w * mf

    # Eq 2
    for c2, v in iteritems(right):
        f = compound_formula[c2.name]
        for e in elements:
            mf = f.get(e, 0)
            delta_sum = delta.sum((c1, c2, e) for c1 in left)
            try:
                p.add_linear_constraints(delta_sum == v * mf)
            except ValueError:
                raise UnbalancedReactionError('Unable to add constraint')

    for c1, v1 in iteritems(left):
        for c2 in right:
            f1 = compound_formula[c1.name]
            f2 = compound_formula[c2.name]
            for e in elements:
                mf = f1.get(e, 0)

                # Eq 4
                p.add_linear_constraints(
                    delta[c1, c2, e] <= float(v1) * mf * omega[c1, c2])

                # Eq 5
                if e in f2:
                    p.add_linear_constraints(
                        delta[c1, c2, e] <= float(v1) * mf * gamma[c1, c2, e])

            # Eq 7, 10
            p.add_linear_constraints(
                m[c1, c2] <= float(v1) * zs[c1] * x[c1, c2],
                q[c1, c2] <= float(v1) * zs[c1] * y[c1, c2])

    try:
        result = p.solve(lp.ObjectiveSense.Maximize)
    except lp.SolverError:
        raise UnbalancedReactionError('Unable to solve')

    first_objective = result.get_value(objective)

    # Yield the first result
    yield dict(_transfers(reaction, delta, elements, result, epsilon))

    # Find alternative solutions
    if alt_elements is None:
        alt_elements = elements
    else:
        alt_elements = elements.intersection(alt_elements)

    if len(alt_elements) == 0:
        return

    # Add a constraint that disallows the current transfer solution until no
    # alternative solutions are left.
    add_objective_constraint = True
    while True:
        elem_vars = 0
        elem_count = 0
        for c1, c2 in product(left, right):
            for e in alt_elements:
                if result.get_value(gamma[c1, c2, e]) > 0.5:
                    elem_count += 1
                    elem_vars += gamma[c1, c2, e]

        if elem_count == 0:
            break

        if add_objective_constraint:
            p.add_linear_constraints(objective >= first_objective)
            add_objective_constraint = False

        p.add_linear_constraints(elem_vars <= elem_count - 1)
        try:
            result = p.solve(lp.ObjectiveSense.Maximize)
        except lp.SolverError:
            break

        yield dict(_transfers(reaction, delta, elements, result, epsilon))