ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/models2010/modelbasic_varyC.py
Revision: 41
Committed: Tue Jan 18 00:35:23 2011 UTC (8 years, 9 months ago) by clausted
File size: 4130 byte(s)
Log Message:
Moved old data class "SPRdataclass" and accompanying surface interaction model modules to /models2010 subdirectory.  The plan is to implement these models for use with the "ba_class" and the modules in the parent directory.  

Should all the models be added to mdl_module or should they each go in their own module?  I am undecided.  
Line File contents
1 """
2 Provide the model for curvefitting.
3 This is a BASIC 1:1 REACTION model with time variable concentration.
4 The reaction equations are:
5 L + A ==ka=> LA
6 L + A <=kd== LA
7
8 Yuhang Wan, Rui Hou
9 Last modified on 100427 (yymmdd) by YW
10
11 Typical Pipeline:(work together with "modelclass.py" the father class):
12 >import modelbasic_varyC as mbvc
13 >m = mbvc.basicmodel_varyC()
14 >m.check_attribute()
15 >m.load_conc('example_conc.txt') or m.create_conc(t)
16 >m.wizard() or m.load_model('model.obj')
17 >m.simulate(t)
18 >m.save_model('newmodel')
19 >m.save_parainfo('parainfo.txt')
20 >m.save_sim_data('sim.dat')
21
22 """
23 __version__ = "100427"
24
25 import numpy as np
26 import pylab as plt
27 import copy
28 import os
29 import pickle
30 from modelclass import *
31
32
33 class basicmodel_varyC(modelclass):
34 """The data attributes of the class:
35 parainfo: the parameter information
36 sim_data: the simulated data
37 conc: the time variable concentration data
38 """
39
40 def __init__ (self, parainfo=[], sim_data=[], conc=[] ):
41 modelclass.__init__(self,parainfo,sim_data)
42 self.conc = conc
43
44 ##------------ not finished yet---------------
45 def load_conc(self, fname):
46 ##fp = file(fname, 'r+')
47 conc = np.loadtxt(fname)
48 self.conc = conc
49 return
50
51 def create_conc(self, t):
52 # create a step-function as concentration curve for debugging now
53 c = np.zeros(len(t), dtype=float)
54 conc = np.zeros((2,len(t)), dtype=float)
55 t = np.array(t)
56 t1 = input("start time:")
57 t2 = input("end time:")
58 ind = plt.find((t>t1-(t[1]-t[0]))&(t<t2-(t[1]-t[0])))
59 c[ind] = 1
60 conc[0] = t
61 conc[1] = c
62 self.conc = conc
63 plt.plot(t,c,'.')
64 plt.show()
65 return
66
67 ##----------------for debugging now------------
68
69 def wizard(self, ):
70 '''This function helps you to create a parameter information list.
71 This is a basic 1:1 reaction model, where there are 4 parameters:
72 rmax: Maximum analyte binding capacity(RU),
73 ka: Association rate constant(M-1S-1),
74 kd: Dissociation rate constant(S-1),
75 ca: Analyte concentration(M).'''
76 print ('')
77 pname = ['rmax','ka','kd','ca']
78 parainfo = []
79 for i in range(4):
80 print i
81 tmp={}
82 tmp['name'] = pname[i]
83 vstr = raw_input("input the value of parameter '%s', if a series, seperate with ',': " %pname[i])
84 vtmp = vstr.strip().split(',')
85 tmp['number'] = len(vtmp)
86 if len(vtmp) == 1:
87 vtmp = float(vtmp[0])
88 else:
89 vtmp = map(float, vtmp)
90 tmp['value'] = vtmp
91 tmp['fixed'] = bool(input("Is '%s' fixed? (1/0): " %pname[i]))
92
93 parainfo.append(tmp)
94 ##parainfo.append({'name':'', 'value':0., 'fixed':0, 'number':0})
95 self.parainfo = parainfo
96
97 return
98
99 def function(self, t, paralist):
100 '''This function calculates the theoretical curve through the
101 parameter list you give.
102 '''
103 # for the basic 1:1 model
104 ## Assign the parameters to calculate the curve
105 conc = copy.deepcopy(self.conc)
106 C = np.interp(t,conc[0],conc[1])
107 for p in paralist:
108 if p['name'] == 'rmax': rmax = p['value']
109 elif p['name'] == 'ka': ka = p['value']
110 elif p['name'] == 'kd': kd = p['value']
111 elif p['name'] == 'ca': ca = p['value']
112 else: print p['name'], p['value']
113 if type(ca) == list:
114 print "Error: This function can only generate data for a single concentration."
115 return
116
117 ## Must iterate through data, numerical integration.
118 g = np.zeros(len(C), dtype=float)
119
120 for i in range(1,len(C)):
121 dG = ka*ca*C[i-1]*(rmax-g[i-1]) - kd*g[i-1]
122 if (abs(g[i]) > 999999999): dG = 0
123
124 g[i] = g[i-1] + dG * (t[i] - t[i-1])
125
126 return g
127