ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/modelbasic_varyC.py
Revision: 25
Committed: Wed Apr 28 20:22:24 2010 UTC (9 years, 5 months ago) by rjaynes
File size: 4130 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 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