from os import listdir
from os.path import join
import warnings
import re

import pandas
import sympy
import scipy.io

import cobra

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)

models = []
model_file = scipy.io.loadmat("m_model_collection/all_models.mat")
for possible_name in model_file.keys():
    if possible_name.startswith("_"):
        continue
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        model = cobra.io.mat.from_mat_struct(model_file[possible_name], model_id=possible_name)
    models.append(model)

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 Value": rational_solution.f,
                           "Decimal Value": float(rational_solution.f),
                           "Total Error": total_error(S, rational_v),
                           # Ensure the upper and lwoer bounds are satisfied.
                           "bounds_sat": 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)

exact_results.rename(columns={"bounds_sat": r"$ub \ge v \ge lb$",
                              "Total Error": r"$\sum\left|\mathbf S \cdot v\right|$"})

errors.T

results.T

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