Run reaction production check method.

def run_reaction_production_check(self, model, solver, threshold,
                                      implicit_sinks=True):
        """Run reaction production check method."""
        prob = solver.create_problem()

        # Create flux variables
        v = prob.namespace()
        for reaction_id in model.reactions:
            lower, upper = model.limits[reaction_id]
            v.define([reaction_id], lower=lower, upper=upper)

        # Build mass balance constraints
        massbalance_lhs = {compound: 0 for compound in model.compounds}
        for spec, value in iteritems(model.matrix):
            compound, reaction_id = spec
            massbalance_lhs[compound] += v(reaction_id) * value

        # Create production variables and apply constraints
        for compound, lhs in iteritems(massbalance_lhs):
            if implicit_sinks:
                # The constraint is merely >0 meaning that we have implicit
                # sinks for all compounds.
                prob.add_linear_constraints(lhs >= 0)
            else:
                prob.add_linear_constraints(lhs == 0)

        confirmed_production = set()
        for reaction in model.reactions:
            if all(c in confirmed_production for c, _ in
                   model.get_reaction_values(reaction)):
                continue

            prob.set_objective(v(reaction))
            for sense in (lp.ObjectiveSense.Maximize,
                          lp.ObjectiveSense.Minimize):
                try:
                    result = prob.solve(sense)
                except lp.SolverError as e:
                    self.fail(
                        'Failed to solve for compound, reaction: {}, {}:'
                        ' {}'.format(compound, reaction, e))

                flux = result.get_value(v(reaction))
                for compound, value in model.get_reaction_values(reaction):
                    if compound in confirmed_production:
                        continue

                    production = 0
                    if sense == lp.ObjectiveSense.Maximize and flux > 0:
                        production = float(value) * flux
                    elif sense == lp.ObjectiveSense.Minimize and flux < 0:
                        production = float(value) * flux

                    if production >= threshold:
                        confirmed_production.add(compound)

        for compound in sorted(model.compounds):
            if compound not in confirmed_production:
                yield compound