diff --git a/cvode_with_sympy.py b/cvode_with_sympy.py index 9ae2a4b..0b8b1f1 100644 --- a/cvode_with_sympy.py +++ b/cvode_with_sympy.py @@ -3,7 +3,7 @@ from helper_funs import * -def parseSym(model_dict): +def parseSym(model_dict): ''' parse model variables and parameters to symbolic variables and calculate derivatives @@ -16,14 +16,14 @@ def parseSym(model_dict): # get model name model_name = cleanModelName(model_dict) - + # define x names - ode_species = [species for species in model_dict['odes']] + ode_species = list(model_dict['odes']) ode_alg_species = model_dict['vars'] # read algebraic equations if 'alg_eqs' in model_dict: - alg_eqs_species = [species for species in model_dict['alg_eqs']] + alg_eqs_species = list(model_dict['alg_eqs']) # define rhs of alg. eqs alg_eqs = [model_dict['alg_eqs'][species] for species in alg_eqs_species] else: @@ -32,7 +32,7 @@ def parseSym(model_dict): # define p names parameters = model_dict['pars'] - + # define t name t_name = 't' @@ -105,27 +105,23 @@ def writeInitSundials(model_dict,xs,ps,fs,xs_alg,gs,atol=1e-6,rtol=1e-6,hmin=0.0 if not os.path.exists('includes'): os.mkdir('./includes') - fname_def = model_name+'_define.c' - fid = open('./includes/'+fname_def, 'w') - define_str = "#define NPARS %d\t\t/* number of parameters */\n" % (len(ps)) - define_str = define_str+"#define NEQ %d\t\t/* number of equations */\n" % len(fs) - fid.write(define_str) - fid.close() - - fname_init = model_name+'_initialize.c' - fid = open('./includes/'+fname_init, 'w') - if (type(atol)!=list): - atol = [atol for x in range(len(fs))] - - init_str = "hmin = RCONST(%e);\t\t/* minimal stepsize */\n" % hmin - init_str = init_str + "hmax = RCONST(%e);\t\t/* maximal stepsize */\n" % hmax - init_str = init_str + "mxsteps = RCONST(%e);\t\t/* maximal number of steps */\n" % mxsteps - init_str = init_str + "reltol = RCONST(%e);\t\t/* scalar relative tolerance */\n" % rtol - for i in range(len(fs)): - init_str = init_str + "Ith(abstol,%d) = RCONST(%e);\t\t/* vector absolute tolerance components */\n" % (i+1, atol[i]) - fid.write(init_str) - fid.close() - + fname_def = f'{model_name}_define.c' + with open(f'./includes/{fname_def}', 'w') as fid: + define_str = "#define NPARS %d\t\t/* number of parameters */\n" % (len(ps)) + define_str = define_str+"#define NEQ %d\t\t/* number of equations */\n" % len(fs) + fid.write(define_str) + fname_init = f'{model_name}_initialize.c' + with open(f'./includes/{fname_init}', 'w') as fid: + if (type(atol)!=list): + atol = [atol for _ in range(len(fs))] + + init_str = "hmin = RCONST(%e);\t\t/* minimal stepsize */\n" % hmin + init_str = init_str + "hmax = RCONST(%e);\t\t/* maximal stepsize */\n" % hmax + init_str = init_str + "mxsteps = RCONST(%e);\t\t/* maximal number of steps */\n" % mxsteps + init_str = init_str + "reltol = RCONST(%e);\t\t/* scalar relative tolerance */\n" % rtol + for i in range(len(fs)): + init_str = init_str + "Ith(abstol,%d) = RCONST(%e);\t\t/* vector absolute tolerance components */\n" % (i+1, atol[i]) + fid.write(init_str) return 1 @@ -134,10 +130,7 @@ def writeOdeSundials(model_dict,xs,ps,fs,xs_alg,gs,checknegative): write ODEs to c-file for sundials ''' - if checknegative: - checknegative_c = "TRUE" - else: - checknegative_c = "FALSE" + checknegative_c = "TRUE" if checknegative else "FALSE" # get model name model_name = cleanModelName(model_dict) @@ -146,82 +139,83 @@ def writeOdeSundials(model_dict,xs,ps,fs,xs_alg,gs,checknegative): gs_c = convToCstr(gs) alg_dict = dict(zip(xs_alg,gs_c)) - fname_ode = model_name+'_ode_f.c' - fid = open('./includes/'+fname_ode, 'w') - - # write header - header = "static int f(realtype t, N_Vector x, N_Vector xdot, void *user_data)\n" - header = header+"{\n\n" - header = header+"\tint i;" - header = header+"\tUserData data;\n" - header = header+"\tbooleantype check_negative = %s;\n" % (checknegative_c) - fid.write(header) - - # write parameter definitions - par_def = "\t/* Extract needed constants from data */\n" - par_def = par_def+"\tdata = (UserData) user_data;\n" - par_def = par_def+"\tdouble p[NPARS];\n" - par_def = par_def+"\tfor (i=0; ip[i];\n\n" - fid.write(par_def) - - for i, p in enumerate(ps): - tmp_str = '\trealtype %s = p[%d];\n' % (p, i) - fid.write(tmp_str) - fid.write('\n') - - # tmp_str = "\tprintf(\"\\n\\n\");\n" - # fid.write(tmp_str) - - # write variables definitions - for i, x in enumerate(xs): - tmp_str = '\trealtype %s = Ith(x,%d);\n' % (x, i+1) - fid.write(tmp_str) - # tmp_str = "\tprintf(\"%%f\\n\", %s);\n" % (x) - # fid.write(tmp_str) - fid.write('\n') + fname_ode = f'{model_name}_ode_f.c' + with open(f'./includes/{fname_ode}', 'w') as fid: + header = ( + "static int f(realtype t, N_Vector x, N_Vector xdot, void *user_data)\n" + + "{\n\n" + ) + + header += "\tint i;" + header += "\tUserData data;\n" + header += "\tbooleantype check_negative = %s;\n" % (checknegative_c) + fid.write(header) + + par_def = ( + "\t/* Extract needed constants from data */\n" + + "\tdata = (UserData) user_data;\n" + ) + + par_def += "\tdouble p[NPARS];\n" + par_def += "\tfor (i=0; ip[i];\n\n" + fid.write(par_def) + + for i, p in enumerate(ps): + tmp_str = '\trealtype %s = p[%d];\n' % (p, i) + fid.write(tmp_str) + fid.write('\n') + # tmp_str = "\tprintf(\"\\n\\n\");\n" + # fid.write(tmp_str) - # write alg eqs definitions - for alg in alg_dict: - tmp_str = '\t%s = %s;\n' % (alg,alg_dict[alg]) - fid.write(tmp_str) - # tmp_str = "\tprintf(\"%%f\\n\", %s);\n" % (alg) - # fid.write(tmp_str) - fid.write('\n') - - # tmp_str = "\tprintf(\"\\n\\n\");\n" - # fid.write(tmp_str) - - # # check for neg species # FIXME! - # non_neg_spec = [sp.origins[model.name] for sp in model.sp_obj if not sp.type == 'thermodynamic'] - # non_neg_spec_cond = ["%s<-1e-5" % (spec) for spec in non_neg_spec] - - # fid.write("\t%s%s%s" % ("if (check_negative && (", " || ".join(non_neg_spec_cond),")) {\n")) - # # fid.write("\t\tprintf(\"%s\\n\", \"At least one species turned negative!\");\n") - # fid.write("\t\treturn(1);\n") - # fid.write("\t}\n") - - # write odes - for i, f in enumerate(fs_c): - if not (f=='0' or f=='0.0' or f==0): - tmp_str = '\tIth(xdot,%d) = %s;\n' % (i+1, f) + # write variables definitions + for i, x in enumerate(xs): + tmp_str = '\trealtype %s = Ith(x,%d);\n' % (x, i+1) fid.write(tmp_str) + # tmp_str = "\tprintf(\"%%f\\n\", %s);\n" % (x) + # fid.write(tmp_str) fid.write('\n') - # write variables definitions - for i, x in enumerate(xs): - if x in xs_alg: - tmp_str = '\tIth(x,%d) = %s;\n' % (i+1, x) + + # write alg eqs definitions + for alg in alg_dict: + tmp_str = '\t%s = %s;\n' % (alg,alg_dict[alg]) fid.write(tmp_str) - # tmp_str = "\tprintf(\"%%f\\n\", %s);\n" % (x) + # tmp_str = "\tprintf(\"%%f\\n\", %s);\n" % (alg) # fid.write(tmp_str) fid.write('\n') - # write footer - footer = '\treturn(0);\n}\n' - fid.write(footer) - fid.close() + # tmp_str = "\tprintf(\"\\n\\n\");\n" + # fid.write(tmp_str) + + # # check for neg species # FIXME! + # non_neg_spec = [sp.origins[model.name] for sp in model.sp_obj if not sp.type == 'thermodynamic'] + # non_neg_spec_cond = ["%s<-1e-5" % (spec) for spec in non_neg_spec] + + # fid.write("\t%s%s%s" % ("if (check_negative && (", " || ".join(non_neg_spec_cond),")) {\n")) + # # fid.write("\t\tprintf(\"%s\\n\", \"At least one species turned negative!\");\n") + # fid.write("\t\treturn(1);\n") + # fid.write("\t}\n") + + # write odes + for i, f in enumerate(fs_c): + if f not in ['0', '0.0', 0]: + tmp_str = '\tIth(xdot,%d) = %s;\n' % (i+1, f) + fid.write(tmp_str) + fid.write('\n') + + # write variables definitions + for i, x in enumerate(xs): + if x in xs_alg: + tmp_str = '\tIth(x,%d) = %s;\n' % (i+1, x) + fid.write(tmp_str) + # tmp_str = "\tprintf(\"%%f\\n\", %s);\n" % (x) + # fid.write(tmp_str) + fid.write('\n') + # write footer + footer = '\treturn(0);\n}\n' + fid.write(footer) return 1 def writeJacSundials(model_dict,xs,ps,fs,xs_alg,gs,dfdx): @@ -235,63 +229,64 @@ def writeJacSundials(model_dict,xs,ps,fs,xs_alg,gs,dfdx): gs_c = convToCstr(gs) alg_dict = dict(zip(xs_alg,gs_c)) - fname_jac = model_name+'_ode_jac.c' - fid = open('./includes/'+fname_jac, 'w') - - # write header - header = "static int Jac(long int N, realtype t, N_Vector x, N_Vector fx, DlsMat J, void *user_data,\nN_Vector tmp1, N_Vector tmp2, N_Vector tmp3)\n" - header = header+"{\n\n" - header = header+"\tint i, j;\n" - header = header+"\tUserData data;\n" - fid.write(header) - - # write parameter definitions - par_def = "\t/* Extract needed constants from data */\n" - par_def = par_def+"\tdata = (UserData) user_data;\n" - par_def = par_def+"\tdouble p[NPARS];\n" - par_def = par_def+"\tfor (i=0; ip[i];\n\n" - fid.write(par_def) - - # write parameter definitions - for i, p in enumerate(ps): - tmp_str = '\trealtype %s = p[%d];\n' % (p, i) - fid.write(tmp_str) - fid.write('\n') - - #write variables definitions - for i, x in enumerate(xs): - tmp_str = '\trealtype %s = Ith(x,%d);\n' % (x, i+1) - fid.write(tmp_str) - fid.write('\n') + fname_jac = f'{model_name}_ode_jac.c' + with open(f'./includes/{fname_jac}', 'w') as fid: + header = ( + "static int Jac(long int N, realtype t, N_Vector x, N_Vector fx, DlsMat J, void *user_data,\nN_Vector tmp1, N_Vector tmp2, N_Vector tmp3)\n" + + "{\n\n" + ) + + header += "\tint i, j;\n" + header += "\tUserData data;\n" + fid.write(header) + + par_def = ( + "\t/* Extract needed constants from data */\n" + + "\tdata = (UserData) user_data;\n" + ) + + par_def += "\tdouble p[NPARS];\n" + par_def += "\tfor (i=0; ip[i];\n\n" + fid.write(par_def) + + # write parameter definitions + for i, p in enumerate(ps): + tmp_str = '\trealtype %s = p[%d];\n' % (p, i) + fid.write(tmp_str) + fid.write('\n') - #write alg eqs - for alg in alg_dict: - tmp_str = '\t%s = %s;\n' % (alg,alg_dict[alg]) - fid.write(tmp_str) - fid.write('\n') - - # write jacobian - for i in range(dfdx.shape[0]): - for j in range(dfdx.shape[1]): - if not dfdx[i,j]==0: - # convert formula to c-style - dfdx_c = convToCstr([dfdx[i,j]]) - tmp_str = '\tIJth(J,%d,%d) = %s;\n' % (i+1,j+1, dfdx_c[0]) - fid.write(tmp_str) - fid.write('\n') + #write variables definitions + for i, x in enumerate(xs): + tmp_str = '\trealtype %s = Ith(x,%d);\n' % (x, i+1) + fid.write(tmp_str) + fid.write('\n') - # write variables definitions - for i, x in enumerate(xs): - if x in xs_alg: - tmp_str = '\tIth(x,%d) = %s;\n' % (i+1, x) + #write alg eqs + for alg in alg_dict: + tmp_str = '\t%s = %s;\n' % (alg,alg_dict[alg]) fid.write(tmp_str) fid.write('\n') - # write footer - footer = '\treturn(0);\n}\n' - fid.write(footer) - fid.close() + # write jacobian + for i in range(dfdx.shape[0]): + for j in range(dfdx.shape[1]): + if dfdx[i, j] != 0: + # convert formula to c-style + dfdx_c = convToCstr([dfdx[i,j]]) + tmp_str = '\tIJth(J,%d,%d) = %s;\n' % (i+1,j+1, dfdx_c[0]) + fid.write(tmp_str) + fid.write('\n') + + # write variables definitions + for i, x in enumerate(xs): + if x in xs_alg: + tmp_str = '\tIth(x,%d) = %s;\n' % (i+1, x) + fid.write(tmp_str) + fid.write('\n') + # write footer + footer = '\treturn(0);\n}\n' + fid.write(footer) return 1 @@ -452,7 +447,7 @@ def writeModelFiles(model_dict,force=False,atol=1e-6,rtol=1e-6,hmin=0.0,hmax=0.0 if force or (len(bin_exists) != 1 or bin_exists[0]==False): # remove old bin files for bin in bin_files: - rm_cmd = "rm -rf %s" % os.path.join('./bin/',bin) + rm_cmd = f"rm -rf {os.path.join('./bin/', bin)}" os.system(rm_cmd) # parse model (xs,ps,fs,xs_alg,gs,dfdx) = parseSym(model_dict) @@ -474,9 +469,7 @@ def objectiveFunction(model_dict,initvars,initpars,data,tExp): model_dict['initvars'] = initvars t,x = integrateSundials(model_dict,tSim=tExp) - chi2 = np.sum((x-data['x'])**2/data['sd']**2) - - return chi2 + return np.sum((x-data['x'])**2/data['sd']**2) def convertToD2D(model_dict,savepath='./D2D'): diff --git a/helper_funs.py b/helper_funs.py index cb17744..5e23124 100644 --- a/helper_funs.py +++ b/helper_funs.py @@ -14,16 +14,16 @@ def removeModelFiles(): # check if bin directory exists if os.path.exists('bin'): for binfile in os.listdir('./bin'): - os.remove('./bin/%s' % binfile) + os.remove(f'./bin/{binfile}') os.rmdir('bin') if os.path.exists('includes'): for includesfile in os.listdir('./includes'): - os.remove('./includes/%s' % includesfile) + os.remove(f'./includes/{includesfile}') os.rmdir('includes') if os.path.exists('src'): for cfile in os.listdir('./src'): - if not cfile in ['integrate_cvode.c','integrate_idas.c']: - os.remove('./src/%s' % cfile) + if cfile not in ['integrate_cvode.c', 'integrate_idas.c']: + os.remove(f'./src/{cfile}') return 1 @@ -46,7 +46,7 @@ def cleanModelName(model_dict): short_name = '' # limit short name to 20 characters for part in model_name.split('_'): - short_name = short_name+'_'+part + short_name = f'{short_name}_{part}' if len(short_name) > 20: break @@ -80,15 +80,13 @@ def cleanFormula(formula): """ # spaces around special characters for s in ['(', ')', '*', '/', '+', '-', ',']: - formula = formula.replace(s, ' %s ' % s) + formula = formula.replace(s, f' {s} ') # reduce spaces ' ' -> ' ' - for i in range(3): + for _ in range(3): formula = formula.replace(" ", " ") - # substitutes xe-xx for xxe - xx - match = re.search(r'\de - ', formula) - if match: - my_digit = match.group(0)[0] - formula = formula.replace(match.group(0), my_digit + 'e-') + if match := re.search(r'\de - ', formula): + my_digit = match[0][0] + formula = formula.replace(match[0], f'{my_digit}e-') # convert string to float with decimal point and back ('1' to '1.0') split_formula = formula.split(' ') n_split_formula = [] @@ -126,12 +124,11 @@ def modelHash(model_dict): """ generate checksum for model dict """ - + model_structure = [model_dict[m] for m in model_dict if m not in ['initpars','initvars']] - checksum = abs(hash(repr(model_structure))) #checksum = abs(reduce(lambda x,y : x^y, [hash(str(item)) for item in model_dict.items()])) - return checksum + return abs(hash(repr(model_structure))) def SBMLwritePythonModule(module_dict,doc_txt=''): @@ -452,16 +449,15 @@ def convToDict(model_name,x_names,p_names,dxdt,x0=[],p0=[]): """ converts lists of variables, parameters and odes to model dictionary """ - model_dict = {} - model_dict['name'] = model_name - model_dict['odes'] = dict(zip(x_names,dxdt)) - model_dict['vars'] = x_names - model_dict['pars'] = p_names - - if x0 == []: - model_dict['initvars'] = dict(zip(x_names, 0.1*np.ones([len(x_names)]))) - else: - model_dict['initvars'] = dict(zip(x_names, x0)) + model_dict = { + 'name': model_name, + 'odes': dict(zip(x_names, dxdt)), + 'vars': x_names, + 'pars': p_names, + 'initvars': dict(zip(x_names, 0.1 * np.ones([len(x_names)]))) + if x0 == [] + else dict(zip(x_names, x0)), + } if p0 == []: model_dict['initpars'] = dict(zip(p_names, 0.1*np.ones([len(p_names)]))) diff --git a/model_def.py b/model_def.py index e5d5a60..7267686 100644 --- a/model_def.py +++ b/model_def.py @@ -15,14 +15,15 @@ def EpoEpoR(): t_name = 't' # define rhs of ODEs - dxdt = [] + dxdt = [ + '- kon * Epo * EpoR + koff * EpoEpoR + kex * EpoEpoRi', + '- kon * Epo * EpoR + koff * EpoEpoR + kt * EpoR0 - kt * EpoR + kex * EpoEpoRi', + 'kon * Epo * EpoR - koff * EpoEpoR - ke * EpoEpoR', + 'ke * EpoEpoR - kex * EpoEpoRi - kdi * EpoEpoRi - kde * EpoEpoRi', + 'kdi * EpoEpoRi', + 'kde * EpoEpoRi', + ] - dxdt.append('- kon * Epo * EpoR + koff * EpoEpoR + kex * EpoEpoRi') - dxdt.append('- kon * Epo * EpoR + koff * EpoEpoR + kt * EpoR0 - kt * EpoR + kex * EpoEpoRi') - dxdt.append('kon * Epo * EpoR - koff * EpoEpoR - ke * EpoEpoR') - dxdt.append('ke * EpoEpoR - kex * EpoEpoRi - kdi * EpoEpoRi - kde * EpoEpoRi') - dxdt.append('kdi * EpoEpoRi') - dxdt.append('kde * EpoEpoRi') return (t_name, x_names, p_names, dxdt) @@ -38,11 +39,12 @@ def ToyModel(): t_name = 't' # define rhs of ODEs - dxdt = [] + dxdt = [ + '- k1 * A + k2 * B - k3 * A * B + k4 * C', + '+ k1 * A - k2 * B - k3 * A * B + k4 * C', + '+ k3 * A * B - k4 * C', + ] - dxdt.append('- k1 * A + k2 * B - k3 * A * B + k4 * C') - dxdt.append('+ k1 * A - k2 * B - k3 * A * B + k4 * C') - dxdt.append('+ k3 * A * B - k4 * C') return (t_name, x_names, p_names, dxdt) @@ -58,11 +60,12 @@ def ToyModel2(): t_name = 't' # define rhs of ODEs - dxdt = [] + dxdt = [ + '- k1 * log ( A ) + k2 * B ** k3 * log ( A * B ) + C**k4', + '+ k1 * A - k2 * B - k3 * A * B + k4 * C', + '+ k3 * A * B - k4 * C', + ] - dxdt.append('- k1 * log ( A ) + k2 * B ** k3 * log ( A * B ) + C**k4') - dxdt.append('+ k1 * A - k2 * B - k3 * A * B + k4 * C') - dxdt.append('+ k3 * A * B - k4 * C') return (t_name, x_names, p_names, dxdt) @@ -94,13 +97,16 @@ def MAPK(): "k10*MAPKp/(KK10+MAPKp)"] # define rhs of ODEs - dxdt = [v[2]+"-"+v[1], # MKKK - v[1]+"-"+v[2], # MKKKp - v[6]+"-"+v[3], # MKK - v[3]+"+"+v[5]+"-"+v[4]+"-"+v[6], # MKKp - v[4]+"-"+v[5], # MKKpp - v[10]+"-"+v[7], # MAPK - v[7]+"+"+v[9]+"-"+v[8]+"-"+v[10], # MAPKp - v[8]+"-"+v[9]] # MAPKpp + dxdt = [ + f"{v[2]}-{v[1]}", + f"{v[1]}-{v[2]}", + f"{v[6]}-{v[3]}", + f"{v[3]}+{v[5]}-{v[4]}-{v[6]}", + f"{v[4]}-{v[5]}", + f"{v[10]}-{v[7]}", + f"{v[7]}+{v[9]}-{v[8]}-{v[10]}", + f"{v[8]}-{v[9]}", + ] + return (t_name, x_names, p_names, dxdt)