Scripts as supporting information to Palm, Rohwer, Hofmeyr (2014)

fit_robust.py

import os
import pod
import random
import csv
import itertools

from time import strftime, sleep
from timeit import default_timer as clock
from operator import attrgetter

begin = clock()

print 'Setting up...'

independent_allos = '(1 + {G6P})**4.0*(1 + {ATP})**4.0'.format(
    G6P='G6P/K_G6P_{state}',
    ATP='ATP/K_ATP_{state}'
)
dependent_allos = '(1 + {G6P} + {ATP})**4.0'.format(
    G6P='G6P/K_G6P_{state}',
    ATP='ATP/K_ATP_{state}'
)

num = {
    'Dep2'       : '{{L}}{V}*4.0*{S}*(1 + {S} + {X})**(3.0)*{allos}'.format(
        S='UDPG/K_UDPG_{state}',
        P='UDP/K_UDP_{state}',
        X='ATP/Kc_ATP_{state}',
        allos=dependent_allos,
        V='{{V}}_{state}'
    ),
    'Dep3' : '{{L}}{V}*4.0*{S}*(1 + {S} + {X})**(3.0)*{allos}'.format(
        S='UDPG/K_UDPG_r',      
        P='UDP/K_UDP_{state}',
        X='ATP/Kc_ATP_{state}',
        allos=dependent_allos,
        V='{{V}}_{state}'
    ),
    'Dep1'           : '{{L}}{V}*4.0*{S}*(1 + {S})**(3.0)*{allos}'.format(
        S='UDPG/K_UDPG_{state}',
        P='UDP/K_UDP_{state}',
        allos=dependent_allos,
        V='{{V}}_{state}'
    ),
    'Indep1'   : '{{L}}{V}*4.0*{S}*(1 + {S})**(3.0)*{allos}'.format(
        S='UDPG/K_UDPG_{state}',
        P='UDP/K_UDP_{state}',
        allos=independent_allos,
        V='{{V}}_{state}'
    ),
    'Indep2'   : '{{L}}{V}*4.0*{S}*(1 + {S})**(3.0)*{allos}'.format(
        S='UDPG/K_UDPG_{state}',
        P='UDP/K_UDP_{state}',
        allos=independent_allos,
        V='{{V}}_{state}'
    ),
}

den = {
    'Dep2'       : '{{L}}(1 + {S} + {X})**4.0*{allos}'.format(
        S='UDPG/K_UDPG_{state}',
        P='UDP/K_UDP_{state}',
        X='ATP/Kc_ATP_{state}',
        allos=dependent_allos
    ),
    'Dep3' : '{{L}}(1 + {S} + {X})**4.0*{allos}'.format(
        S='UDPG/K_UDPG_r',      
        P='UDP/K_UDP_{state}',
        X='ATP/Kc_ATP_{state}',
        allos=dependent_allos
    ),
    'Dep1'           : '{{L}}(1 + {S})**4.0*{allos}'.format(
        S='UDPG/K_UDPG_{state}',
        P='UDP/K_UDP_{state}',
        allos=dependent_allos
    ),
    'Indep1'   : '{{L}}(1 + {S})**4.0*{allos}'.format(
        S='UDPG/K_UDPG_{state}',
        P='UDP/K_UDP_{state}',
        allos=independent_allos
    ),
    'Indep2'   : '{{L}}(1 + {S})**4.0*{allos}'.format(
        S='UDPG/K_UDPG_{state}',
        P='UDP/K_UDP_{state}',
        allos=independent_allos
    ),
}


def buildEquation(id, states, only_r_active=False):
    equation = '(' + num[id].format(L='', state='r')

    if not only_r_active:
        for s in states:
            L = '{L}_' + s + '*'
            equation = equation + ' + ' + num[id].format(L=L, state=s)

    equation = equation + ')/(' + den[id].format(L='', state='r')

    for s in states:
        L = '{L}_' + s + '*'
        equation = equation + ' + ' + den[id].format(L=L, state=s)

    equation = equation + ')'

    if id == 'Indep2':
        equation = equation.replace('ATP/K_ATP_r', '0')

    return equation


equationIds = ['Dep1', 'Dep2', 'Dep3', 'Indep1', 'Indep2']

equations = {}

for i in equationIds:
    equations[i + '_2'] = buildEquation(i, ['t'])
    equations[i + '_2_only_r_active'] = buildEquation(i, ['t'], only_r_active=True)


def get_equation(equation, V, L):
    return equation.format(V=V, L=L)

working_dir = ''
# working_dir = os.getenv('HOME') + '/Documents/Engauge/'

I_UDPG_datasource = pod.CSV_ScanDataSource(
    working_dir + 'PirRotCab1968_I_UDPG.csv',
    ['UDPG', 'ATP', 'G6P'],
    'v'
)
I_G6P_datasource = pod.CSV_ScanDataSource(
    working_dir + 'PirRotCab1968_I_G6P.csv',
    ['G6P', 'ATP', 'UDPG'],
    'v'
)
D_UDPG_datasource = pod.CSV_ScanDataSource(
    working_dir + 'PirRotCab1968_D_UDPG.csv',
    ['UDPG', 'ATP', 'G6P'],
    'v'
)
D_G6P_datasource = pod.CSV_ScanDataSource(
    working_dir + 'PirRotCab1968_D_G6P.csv',
    ['G6P', 'ATP', 'UDPG'],
    'v'
)

from multiprocessing import Pool


def solve_all(fitter):
    r = fitter.solve()
    r._prepareForPickle()
    return r

pool = Pool(processes=4)

tasks = []

robust = None

count = 0

for eq in equations:
    equation_I_UDPG = get_equation(equations[eq], V='V_I', L='L_I')
    equation_D_UDPG = get_equation(equations[eq], V='V_D', L='L_D')
    equation_I_G6P = get_equation(equations[eq], V='V_I', L='L_I')
    equation_D_G6P = get_equation(equations[eq], V='V_D', L='L_D')

    I_G6P_model = pod.Equation_Model(equation_I_G6P, ['G6P', 'ATP', 'UDPG'])
    I_UDPG_model = pod.Equation_Model(equation_I_UDPG, ['UDPG', 'ATP', 'G6P'])
    D_G6P_model = pod.Equation_Model(equation_D_G6P, ['G6P', 'ATP', 'UDPG'])
    D_UDPG_model = pod.Equation_Model(equation_D_UDPG, ['UDPG', 'ATP', 'G6P'])

    for i in range(1000):
        f = pod.Fitter('PirRotCab1968_{0}'.format(eq))

        pod.ScanDataSet('I_G6P', f, I_G6P_datasource, I_G6P_model)
        pod.ScanDataSet('I_UDPG', f, I_UDPG_datasource, I_UDPG_model)
        pod.ScanDataSet('D_G6P', f, D_G6P_datasource, D_G6P_model)
        pod.ScanDataSet('D_UDPG', f, D_UDPG_datasource, D_UDPG_model)

        # Specify parameters
        f.addParameter('K_UDPG_r', init=random.uniform(0.1, 1.0), min=1e-25)
        if not 'Dep3' in eq:
            f.addParameter('K_UDPG_t', init=random.uniform(0.1, 1.0), min=1e-25)
        f.addParameter('K_G6P_r', init=random.uniform(0.001, 0.01), min=1e-25)
        f.addParameter('K_G6P_t', init=random.uniform(0.01, 0.1), min=1e-25)
        if not 'Indep2' in eq:
            f.addParameter('K_ATP_r', init=random.uniform(0.1, 5.0), min=1e-25)
        f.addParameter('K_ATP_t', init=random.uniform(0.1, 5.0), min=1e-25)
        if 'Dep2' in eq:
            f.addParameter('Kc_ATP_r', init=random.uniform(0.1, 5.0), min=1e-25)
            f.addParameter('Kc_ATP_t', init=random.uniform(0.1, 5.0), min=1e-25)
        f.addParameter('L_I_t', init=random.uniform(0.1, 10.0), min=1e-25)
        f.addParameter('L_D_t', init=random.uniform(10.0, 1000.0), min=1e-25)
        f.addParameter('V_I_r', init=random.uniform(0.1, 0.5), min=1e-25)
        f.addParameter('V_D_r', init=random.uniform(0.1, 0.5), min=1e-25)
        if not 'only_r_active' in eq:
            f.addParameter('V_I_t', init=random.uniform(0.0, 1e-25), min=1e-25)
            f.addParameter('V_D_t', init=random.uniform(0.0, 1e-25), min=1e-25)

        alg = pod.robust_biweight()

        f.setAlgorithm(alg)

        tasks.append(f)

        count += 1

print 'Submitted tasks: {0}'.format(count)

done = pool.map(solve_all, tasks)  # Parallel map

resdict = {}

def process(r):
    if r.success:
        key = r.name
        if key in resdict:
            res = resdict[key]
            if (r.fitted_parameters > 0).all():
                res.append(r)
            res = sorted(res, key=attrgetter('r2'), reverse=True)
            if len(res) >= 50:
                res = res[:50]
            resdict[key] = res
        else:
            resdict[key] = [r]

map(process, done)

print 'Processing...'

reslist = []

for k in resdict:
    # Save the results for the best fit of each equation
    res = resdict[k]
    r = res[0]
    r._initAfterPickle()
    print '{0} : {1}'.format(k, r.r2)
    r.addOutputClass(pod.FullDataTableOutput)
    r.addOutputClass(pod.ResidualPlotOutput)
    r.addOutputClass(pod.PgfPlots3DOutput, reverse_x=True, reverse_y=True)
    r.addOutputClass(pod.ParametersTxtOutput)
    r.writeOutput()
    reslist.append(r)

    # Write all the results to CSV
    writeHeader = True
    for res in resdict[k]:
        if not r.success:
            continue

        csv_file = open('pod_results/{0}.csv'.format('PirRotCab1968_robust_{0}'.format(k)), 'a')
        csv_writer = csv.writer(csv_file)

        if writeHeader:
            row = []
            row.extend(['id', 'R2', 'R2_adj', 'numFuncEval', 'SS_err', 'SS_total', 'n', 'p'])

            for p in r.parameters:
                row.append(p.name)
            csv_writer.writerow(row)
            writeHeader = False

        row = []
        row.extend([
            res.name,
            res.r2,
            res.r2_adj,
            res.num_func_evals,
            res.ss_err,
            res.ss_tot,
            len(res.residuals),
            len(res.parameters)
        ])

        for p in res.fitted_parameters:
            row.append(p)

        csv_writer.writerow(row)
        csv_file.close()

end = clock()
print 'Done, took {0} seconds'.format(end-begin)

ftestpairs = itertools.combinations(reslist, 2)


def ftest(r):
    f = pod.FTester(0.99)
    f.one_sided(r[0], r[1])
    return f

fres = map(ftest, ftestpairs)

for fr in fres:
    print '{} vs {}:'.format(fr.r1.name, fr.r2.name)
    print fr.winner.name
    print '{} : {}'.format(fr.n1, fr.n2)
    print '{} : {}'.format(fr.p1, fr.p2)
    print '{} : {}'.format(fr.sse1, fr.sse2)
    print 'f-stat: {}'.format(fr.f_stat)
    print 'f-crit: {}'.format(fr.f_crit)
    print '\n'

fv_vs_dop.py

import numpy
import scipy
from openopt import GLP
import itertools
import multiprocessing

RANDOM_18 = 0
SEQUENTIAL_18 = 1

method = SEQUENTIAL_18

s = 0.2
s_std = 0.2
atp = 0.0
g6p_max = 10.0

KusUDPGusr = 0.421297274949
KusGsixPusr = 0.0592919176332
KusGsixPust = 0.390755000664
KusATPusr = 4.66652071316
KusATPust = 2.95378015821
KcusATPusr = 17.9751827042
KcusATPust = 2.86243403008

previous_best = numpy.zeros(18)
previous_best[0] = 1.0
previous_best_min = numpy.zeros(18)
previous_best_min[17] = 1.0

alpha = {}
alpha[''] = 0.250198199
alpha['2'] = 4.8915760016
alpha['2a'] = 4.0035073056
alpha['5'] = 1.0
alpha['4'] = 1.0
alpha['3c'] = 1.0
alpha['3b'] = 2.6673715616
alpha['3a'] = 13.5717499383
alpha['1a'] = 1.0
alpha['1b'] = 1.0

seq2 = ['', '2']
seq543 = ['', '5', '4', '3c', '3b', '3a']
seq3a = ['', '3a']
seq3b = ['', '3b']
seq3c = ['', '3c']
seq4 = ['', '4']
seq5 = ['', '5']
seq1a = ['', '1a']
seq1b = ['', '1b']
seq2a = ['', '2a']
seq22a = ['', '2', '2a']
sites = 7

if method == SEQUENTIAL_18:
    sequences = [seq22a, seq543]
elif method == RANDOM_18:
    sequences = [seq2, seq2a, seq3a, seq3b, seq3c, seq4, seq5]

def state(seq):
    state = []
    for i in range(len(seq)):
        inner = []
        for j in range(i+1):
            inner.append(seq[j])
        state.append(inner)

    return state

cluster_states = map(state, sequences)
states = [i for i in itertools.product(*cluster_states)]
statesx = [set(itertools.chain(*i)) for i in states]
raw_input('Pause')


def L_func(state):
    innerL = 1.0
    for i in state:
        innerL = innerL * alpha[i]
    return innerL

L_values = map(L_func, statesx)


def pdeg_num(c, s):
    return c * (len(s)-1)


def pdeg_func(gs_concentrations):
    num = map(pdeg_num, gs_concentrations, statesx)
    num = numpy.sum(num)

    den = numpy.sum(gs_concentrations)

    ratio = num/den

    return ratio


def v_func(a, E, L, s):
    return (1.0*E*4.0*s/KusUDPGusr*(1 + s/KusUDPGusr + atp/KcusATPusr)**(3.0)*(1 + a/KusGsixPusr + atp/KusATPusr)**4.0)/((1 + s/KusUDPGusr + atp/KcusATPusr)**4.0*(1 + a/KusGsixPusr + atp/KusATPusr)**4.0 + L*(1 + s/KusUDPGusr + atp/KcusATPust)**4.0*(1 + a/KusGsixPust + atp/KusATPust)**4.0)


def fv_func(gs, g6p_min):
    v_num = 0.0
    v_den = 0.0

    for g, l in itertools.izip(gs, L_values):
        v_num += v_func(g6p_min, g, l, s)
        v_den += v_func(g6p_max, g, l, s_std)

    v = v_num/v_den

    return v


def pdeg_fv_pair(g6p_min, GS):
    return pdeg_func(GS), fv_func(GS, g6p_min)


def maximize_fv(goal_pdeg):
    parlength = len(statesx)

    f = lambda x: func(x, goal_pdeg)

    problem = GLP(f, numpy.zeros(parlength), maxNonSuccess=15, lb=numpy.zeros(parlength), ub=numpy.ones(parlength)*50, maxFunEvals=1e7, maxIter=700, storeIterPoints=False)
    ooresult = problem.solve('de')

    result = ooresult.xf
    del problem
    del ooresult
    return result


def func(GS, goal_pdeg):
    pdeg = pdeg_func(GS)
    return numpy.abs(pdeg - goal_pdeg)


master_GS = []

p_degrees = scipy.linspace(0.0, sites, 20, endpoint=True)

pool = multiprocessing.Pool()

master_GS = pool.map(maximize_fv, p_degrees)

k = 0

G6P_concentrations = [0.0, 0.1, 0.25, 0.5]

if method == SEQUENTIAL_18:
    f = open('/home/danie/Documents/LatexProjects/gs_rate_equation/scripts/guisalmas_sequential.dat', 'w')
elif method == RANDOM_18:
    f = open('/home/danie/Documents/LatexProjects/gs_rate_equation/scripts/guisalmas_random.dat', 'w')
f.write('#{0}\t{1}'.format('pdeg', 'fv'))
f.write('\n')

for i in G6P_concentrations:
    master_Pdeg = []
    master_FV = []

    for j in master_GS:
        pdeg, fv = pdeg_fv_pair(i, j)
        master_Pdeg.append(pdeg)
        master_FV.append(fv)

    for x, y in zip(master_Pdeg, master_FV):
        f.write('{0}\t{1}'.format(x, y))
        f.write('\n')
    f.write('\n')
    f.write('\n')

f.close()

PirRotCab1968_D_G6P.csv

PirRotCab1968_D_UDPG.csv

PirRotCab1968_I_G6P.csv

PirRotCab1968_I_UDPG.csv