#!/usr/bin/env python
# coding: utf-8

# # Solving models with rational and floating point.

# Using [cobrapy](http://opencobra.github.io/cobrapy/) (version 0.4.0b1 or later), we have interfaces to serveral floating point solvers ([gurobi](http://www.gurobi.com), [MOSEK](http://www.mosek.com/), [CPLEX](http://www-01.ibm.com/software/commerce/optimization/cplex-optimizer/), [Clp](https://projects.coin-or.org/Clp), and two different versions of [GLPK](https://www.gnu.org/software/glpk/)). Additionally, we can use the rational solving capabilities of GLPK and esolver (also known as [QSopt_ex](http://www.dii.uchile.cl/~daespino/ESolver_doc/main.html)). We will use all of these solvers on all the models in the collection.

# In[1]:


import os

from numpy import array
import pandas
import sympy

import cobra

pandas.set_option("display.max_rows", 100)


# We will verify all of our solutions by writing out the Stoichiometrix matrix $\mathbf S$ using the [sympy](http://sympy.org/) symbolic math library. The total error will be equivalent to 
# $$ \sum \left|\mathbf S \cdot v \right| $$

# In[2]:


def convert_to_rational(value):
    return sympy.Rational("%.15g" % value)

def construct_exact_S(model):
    # intialize to 0
    S = sympy.SparseMatrix(len(model.metabolites), len(model.reactions), 0)
    # populate with stoichiometry
    for i, r in enumerate(model.reactions):
        for met, value in r._metabolites.iteritems():
            S[model.metabolites.index(met), i] = convert_to_rational(value)
    return S

def total_error(S, v):
    return sum(abs(i) for i in S * v)


# Some of these models were exported from [BiGG](http://bigg.ucsd.edu/bigg/export.pl) as SBML, and others were downloaded from their respective publications. All of the models are available in a [git repository](https://github.com/opencobra/m_model_collection). They were parsed into a single mat file by an [included script](load_models.ipynb).

# In[3]:


models = []
for filename in sorted(os.listdir("sbml3")):
    if not filename.endswith(".xml"):
        continue
    models.append(cobra.io.read_sbml_model(os.path.join("sbml3", filename)))


# These models all compute rationally with esolver and also with the other floating point solvers, as shown below. They also compute with floating point solvers in the [COBRA toolbox](http://dx.doi.org/doi:10.1038/nprot.2011.308), as shown [here](http://nbviewer.ipython.org/github/opencobra/m_model_collection/blob/master/MATLAB.ipynb).

# In[4]:


results = {}
exact_results = {}
errors = {}
for m in models:
    S = construct_exact_S(m)
    S_float = m.to_array_based_model().S
    rational_solution = m.optimize(solver="esolver", rational_solution=True)
    rational_v = sympy.Matrix(rational_solution.x)
    exact_results[m.id] = {"Rational": rational_solution.f,
                           "Decimal": float(rational_solution.f),
                           "Error": total_error(S, rational_v),
                           # Ensure the upper and lwoer bounds are satisfied.
                           "Bounds": all([r.upper_bound >= value >= r.lower_bound
                                              for r, value in zip(m.reactions, rational_v)])}
    # solve this model with all the solvers
    solutions = {solver: m.optimize(solver=solver)
                 for solver in cobra.solvers.solver_dict}
    solutions["cglpk_exact"] = m.optimize(solver="cglpk", exact=True)
    # store the objective value and errors
    results[m.id] = {k: v.f for k, v in solutions.iteritems()}
    errors[m.id] = {k: sum(abs(S_float * array(v.x))) for k, v in solutions.iteritems()}
# format the results as pandas dataframes
exact_results = pandas.DataFrame.from_dict(exact_results).T
results = pandas.DataFrame.from_dict(results)
errors = pandas.DataFrame.from_dict(errors)


# For all of these models, we can demonstrate they satisfy the model constraints using exact operations.

# $$\sum\left|\mathbf S \cdot v\right| = 0$$

# In[5]:


abs(exact_results.Error).max()


# $$ub \ge v \ge lb$$

# In[6]:


exact_results.Bounds.all()


# Here are objective values of the rational results provided by esolver:

# In[7]:


exact_results[["Decimal", "Rational"]]


# Here is the $\sum\left|\mathbf S \cdot v\right|$ error computed using floating point operations for every solver. When computed rationally with esolver above, this value was exactly 0. Howver, when rounding the fractional values to floating point, there is a very small amount of resulting error, so even esolver does not give 0 error for this computation.

# In[8]:


errors.T


# Here are all the computed biomass flux rates for each solver and model.

# In[9]:


results.T


# The objectives computed for these solvers are effectively the same as those computed by esolver.

# In[10]:


differences = (results - results.ix["esolver"]).T
differences.pop("esolver")
differences


# In[11]:


abs(differences).max().max()