ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/modelparallel.py
Revision: 25
Committed: Wed Apr 28 20:22:24 2010 UTC (9 years, 6 months ago) by rjaynes
File size: 4753 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 User Rev File contents
1 rjaynes 25 """
2     Provide the model for curvefitting.
3     This is a PARALLEL REACTION model.
4     The reaction equations:
5     L1 + A ==ka1=> L1A
6     L1 + A <=kd1== L1A
7     L2 + A ==ka2=> L2A
8     L2 + A <=kd2== L2A
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 modelparallel as mp
15     >m = mp.parallelmodel()
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    
29     class parallelmodel(modelclass):
30     """The data attributes of the class:
31     parainfo: the parameter information
32     sim_data: the simulated data
33     """
34    
35     def __init__ (self, parainfo=[], sim_data=[] ):
36     modelclass.__init__(self,parainfo,sim_data)
37    
38     def wizard(self, ):
39     '''This function helps you to create a parameter information list.
40     This is a competing reactions model, where there are 10 parameters:
41     rmax1: Maximum analyte binding capacity for L1(RU),
42     ka1: Association rate constant for L1+A=L1A(M-1S-1),
43     kd1: Dissociation rate constant for L1A=L1+A(S-1),
44     rmax2: Maximum analyte binding capacity for L2(RU),
45     ka2: Association rate constant for L2+A=L2A(M-1S-1),
46     kd2: Dissociation rate constant for L2A=L2+A(S-1),
47     ca: Analyte concentration(M),
48     ton: Starting time for sample injection(s),
49     toff: Ending time for sample injection(s),
50     tout: Total time(s).'''
51     print ('')
52     pname = ['rmax1','ka1','kd1','rmax2','ka2','kd2','ca','ton','toff','tout']
53     parainfo = []
54     for i in range(10):
55     print i
56     tmp={}
57     tmp['name'] = pname[i]
58     vstr = raw_input("input the value of parameter '%s', if a series, seperate with ',': " %pname[i])
59     vtmp = vstr.strip().split(',')
60     tmp['number'] = len(vtmp)
61     if len(vtmp) == 1:
62     vtmp = float(vtmp[0])
63     else:
64     vtmp = map(float, vtmp)
65     tmp['value'] = vtmp
66     tmp['fixed'] = bool(input("Is '%s' fixed? (1/0): " %pname[i]))
67    
68     parainfo.append(tmp)
69     ##parainfo.append({'name':'', 'value':0., 'fixed':0, 'number':0})
70     self.parainfo = parainfo
71    
72    
73     def function(self, t, paralist):
74     '''This function calculates the theoretical curve through the
75     parameter list you give.
76     '''
77     # for parallel reactions model
78    
79     ## Assign the parameters to calculate the curve
80     for p in paralist:
81     if p['name'] == 'rmax1': rmax1 = p['value']
82     elif p['name'] == 'ka1': ka1 = p['value']
83     elif p['name'] == 'kd1': kd1 = p['value']
84     elif p['name'] == 'rmax2': rmax2 = p['value']
85     elif p['name'] == 'ka2': ka2 = p['value']
86     elif p['name'] == 'kd2': kd2 = p['value']
87     elif p['name'] == 'ca': ca = p['value']
88     elif p['name'] == 'ton': ton = p['value']
89     elif p['name'] == 'toff': toff = p['value']
90     elif p['name'] == 'tout': tout = p['value']
91     else: print p['name'], p['value']
92     if type(ca1) == list or type(ca2) == list:
93     print "Error: This function can only generate data for a single concentration."
94     return
95    
96     t1 = ton ##100 ##234 # Initial injection
97     t2 = toff ##178 ##570 # Wash, begin of the dissociation
98     t3 = tout ##400 ##840 # End wash
99    
100     ## Must iterate through data, numerical integration.
101     g = np.zeros(len(t), dtype=float)
102     g1 = np.zeros(len(t), dtype=float)
103     g2 = np.zeros(len(t), dtype=float)
104    
105     for i in range(1,len(t)):
106     if (t[i] > t3): break #Speed things up.
107     if (t1 < t[i] < t2):
108     ## Association
109     dG1 = ka1*ca*(rmax1-g1[i-1]) - kd1*g1[i-1]
110     dG2 = ka2*ca*(rmax2-g2[i-1]) - kd2*g2[i-1]
111     elif (t2 < t[i] < t3):
112     ## Dissociation
113     dG1 = 0 - kd1*g1[i-1]
114     dG2 = 0 - kd2*g2[i-1]
115     else:
116     dG1 = 0
117     dG2 = 0
118     if (abs(g1[i]) > 999999999): dG1 = 0
119     if (abs(g2[i]) > 999999999): dG2 = 0
120    
121     g1[i] = g1[i-1] + dG1 * (t[i] - t[i-1])
122     g2[i] = g2[i-1] + dG2 * (t[i] - t[i-1])
123     g[i] = g1[i] + g2[i]
124    
125     return g
126    
127    
128     ##### End of parallel model class definition.