ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/modelcompeting_varyC.py
Revision: 25
Committed: Wed Apr 28 20:22:24 2010 UTC (9 years ago) by rjaynes
File size: 5205 byte(s)
Log Message:
Add py and obj files to allow modeling of more SPR experiments with converter and curvefitting modules.  This is the work of Yuhang Wan and Rui Hou.

1. In "converter.py": 
      Add the saving and reading function for the sprclass data object.
      Also add function "keyfile_read_fake" to provide default information for SPRit and ICM formats in case of the bug when do background_subtract.
      Fix the bugs in "background_subtract".
      Tested by DAM and ICM formats.
2. In model modules:
      "modelclass.py" is the parent class for all the other model classes that performs the theoretical simulating, loading and saving of the parameter or simulated data. Rui and I also add some other model modules like competing model, twostate model, parallel model, and the time variable concentrated models, where the simulated result is compared with Clamp's simulation to make sure the equations are correct. 
       The basicmodel and basicmodel_varyC class are tested. 
3. In "curvefitting.py":
      Add typical pipeline for operation. The examples are packed with the file. 
      Add function to show the Elapsed time for each fitting.
Line File contents
1 """
2 Provide the model for curvefitting.
3 This is a COMPETING REACTION model with time variable concentration.
4 The reaction equations are:
5 L + A1 ==ka1=> LA1
6 L + A1 <=kd1== LA1
7 L + A2 ==ka2=> LA2
8 L + A2 <=kd2== LA2
9
10 Rui Hou, Yuhang Wan
11 Last modified on 100427 (yymmdd) by YW
12
13 Typical Pipeline:(work together with "modelclass.py" the father class):
14 >import modelcompeting_varyC as mcvc
15 >m = mcvc.competingmodel_varyC()
16 >...
17
18 """
19 __version__ = "100427"
20
21 import numpy as np
22 import pylab as plt
23 import copy
24 import os
25 import pickle
26 from modelclass import *
27
28 class competingmodel_varyC(modelclass):
29 """The data attributes of the class:
30 parainfo: the parameter information
31 sim_data: the simulated data
32 conc: the time variable concentration data
33 """
34
35 def __init__ (self, parainfo=[], sim_data=[], conc=[] ):
36 modelclass.__init__(self,parainfo,sim_data)
37 self.conc = conc
38
39 ##------------ not finished yet---------------
40 def load_conc(self, fname):
41 ##fp = file(fname, 'r+')
42 conc = np.loadtxt(fname)
43 self.conc = conc
44 return
45
46 def create_conc(self, t):
47 # create a step-function as concentration curve for debugging now
48 c = np.zeros(len(t), dtype=float)
49 conc = np.zeros((2,len(t)), dtype=float)
50 t = np.array(t)
51 t1 = input("start time:")
52 t2 = input("end time:")
53 ind = plt.find((t>t1-(t[1]-t[0]))&(t<t2-(t[1]-t[0])))
54 c[ind] = 1
55 conc[0] = t
56 conc[1] = c
57 self.conc = conc
58 plt.plot(t,c,'.')
59 plt.show()
60 return
61
62 ##----------------for debugging now------------
63
64
65 def wizard(self, ):
66 '''This function helps you to create a parameter information list.
67 This is a competing reactions model, where there are 9 parameters:
68 rmax: Maximum analyte binding capacity(RU),
69 mw1: (Relative) molecular weight of A1,
70 ka1: Association rate constant for L+A1=LA1(M-1S-1),
71 kd1: Dissociation rate constant for LA1=L+A1(S-1),
72 mw2: (Relative) molecular weight of A2,
73 ka2: Association rate constant for L+A2=LA2(M-1S-1),
74 kd2: Dissociation rate constant for LA2=L+A2(S-1),
75 ca1: Analyte concentration for A1(M),
76 ca2: Analyte concentration for A2(M).'''
77 print ('')
78 pname = ['rmax','mw1','ka1','kd1','mw2','ka2','kd2','ca1','ca2']
79 parainfo = []
80 for i in range(9):
81 print i
82 tmp={}
83 tmp['name'] = pname[i]
84 vstr = raw_input("input the value of parameter '%s', if a series, seperate with ',': " %pname[i])
85 vtmp = vstr.strip().split(',')
86 tmp['number'] = len(vtmp)
87 if len(vtmp) == 1:
88 vtmp = float(vtmp[0])
89 else:
90 vtmp = map(float, vtmp)
91 tmp['value'] = vtmp
92 tmp['fixed'] = bool(input("Is '%s' fixed? (1/0): " %pname[i]))
93
94 parainfo.append(tmp)
95 ##parainfo.append({'name':'', 'value':0., 'fixed':0, 'number':0})
96 self.parainfo = parainfo
97
98
99 def function(self, t, paralist):
100 '''This function calculates the theoretical curve through the
101 parameter list you give.
102 '''
103 # for competing reactions 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'] == 'mw1': mw1 = p['value']
110 elif p['name'] == 'ka1': ka1 = p['value']
111 elif p['name'] == 'kd1': kd1 = p['value']
112 elif p['name'] == 'mw2': mw2 = p['value']
113 elif p['name'] == 'ka2': ka2 = p['value']
114 elif p['name'] == 'kd2': kd2 = p['value']
115 elif p['name'] == 'ca1': ca1 = p['value']
116 elif p['name'] == 'ca2': ca2 = p['value']
117 else: print p['name'], p['value']
118 if type(ca1) == list or type(ca2) == list:
119 print "Error: This function can only generate data for a single concentration."
120 return
121
122 rau = float(mw1/mw2)
123
124 ## Must iterate through data, numerical integration.
125 g = np.zeros(len(C), dtype=float)
126 g1 = np.zeros(len(C), dtype=float)
127 g2 = np.zeros(len(C), dtype=float)
128
129 for i in range(1,len(C)):
130 dG1 = ka1*ca1*C[i-1]*(rmax-g1[i-1]-g2[i-1]*rau) - kd1*g1[i-1]
131 dG2 = ka2*ca2*C[i-1]*(rmax/rau-g1[i-1]/rau-g2[i-1]) - kd2*g2[i-1]
132
133 if (abs(g1[i]) > 999999999): dG1 = 0
134 if (abs(g2[i]) > 999999999): dG2 = 0
135
136 g1[i] = g1[i-1] + dG1 * (t[i] - t[i-1])
137 g2[i] = g2[i-1] + dG2 * (t[i] - t[i-1])
138 g[i] = g1[i] + g2[i]
139
140 return g
141
142
143 ##### End of time variable concentrated competing model class definition.