osprai/trunk/modeltwostate_varyC.py

"""
Provide the model for curvefitting.
This is a TWO STATE REACTION model with time variable concentration..
The reaction equations are:
L + A ==ka1=> LA
L + A <=kd1== LA
LA ==ka2=> LA*
LA <=kd2== LA*

Rui Hou, Yuhang Wan
Last modified on 100427 (yymmdd) by YW

Typical Pipeline:(work together with "modelclass.py" the father class):
>import modeltwostate_varyC as mtvc
>m = mtvc.twostatemodel_varyC()
>...

"""
__version__ = "100427"

import numpy as np
import pylab as plt
import copy
import os
import pickle
from modelclass import *


class twostatemodel_varyC(modelclass):
    """The data attributes of the class:
    parainfo:           the parameter information
    sim_data:           the simulated data    
    conc:               the time variable concentration data
    """
    def __init__ (self, parainfo=[],  sim_data=[], conc=[] ):
        modelclass.__init__(self,parainfo,sim_data)
        self.conc = conc

##------------ not finished yet---------------        
    def load_conc(self, fname):
        ##fp = file(fname, 'r+')
        conc = np.loadtxt(fname)
        self.conc = conc
        return
        
    def create_conc(self, t):
        # create a step-function as concentration curve for debugging now
        c = np.zeros(len(t), dtype=float)
        conc = np.zeros((2,len(t)), dtype=float)
        t = np.array(t)
        t1 = input("start time:")
        t2 = input("end time:")
        ind = plt.find((t>t1-(t[1]-t[0]))&(t<t2-(t[1]-t[0])))
        c[ind] = 1
        conc[0] = t
        conc[1] = c
        self.conc = conc
        plt.plot(t,c,'.')
        plt.show()
        return
        
##----------------for debugging now------------

    def wizard(self, ):
        '''This function helps you to create a parameter information list.
        This is a competing reactions model, where there are 6 parameters: 
        rmax:   Maximum analyte binding capacity(RU),
        ka1:    Association rate constant for L+A=LA(M-1S-1),
        kd1:    Dissociation rate constant for LA=L+A(S-1),
        ka2:    Forward rate constant for LA=LA*(S-1),
        kd2:    Backward rate constant for LA*=LA(S-1),
        ca:     Analyte concentration(M).'''
        print ('')
        pname = ['rmax','ka1','kd1','ka2','kd2','ca']
        parainfo = []
        for i in range(6):
            print i
            tmp={}
            tmp['name'] = pname[i]
            vstr = raw_input("input the value of parameter '%s', if a series, seperate with ',': " %pname[i])
            vtmp = vstr.strip().split(',')
            tmp['number'] = len(vtmp)
            if len(vtmp) == 1:  
                vtmp = float(vtmp[0])
            else:
                vtmp = map(float, vtmp)
            tmp['value'] = vtmp
            tmp['fixed'] = bool(input("Is '%s' fixed? (1/0): " %pname[i]))

            parainfo.append(tmp)
            ##parainfo.append({'name':'', 'value':0., 'fixed':0, 'number':0}) 
        self.parainfo = parainfo
        
        
    def function(self, t, paralist):
        '''This function calculates the theoretical curve through the 
        parameter list you give.
        '''
        # for two state reaction model
        conc = copy.deepcopy(self.conc)        
        C = np.interp(t,conc[0],conc[1])
        ## Assign the parameters to calculate the curve
        for p in paralist:
            if p['name'] == 'rmax':    rmax = p['value']
            elif p['name'] == 'ka1':    ka1 = p['value']
            elif p['name'] == 'kd1':    kd1 = p['value']
            elif p['name'] == 'ka2':    ka2 = p['value']
            elif p['name'] == 'kd2':    kd2 = p['value']
            elif p['name'] == 'ca':     ca = p['value']
            else: print p['name'], p['value']
        if type(ca1) == list or type(ca2) == list:
            print "Error: This function can only generate data for a single concentration."
            return            
    
        ## Must iterate through data, numerical integration.
        g = np.zeros(len(C), dtype=float)
        g1 = np.zeros(len(C), dtype=float)
        g2 = np.zeros(len(C), dtype=float)    

        for i in range(1,len(C)):
            dG1 = (ka1*ca*C[i-1]*(rmax-g1[i-1]-g2[i-1]) - kd1*g1[i-1])-(ka2*g1[i-1]- kd2*g2[i-1])
            dG2 = ka2*g1[i-1]- kd2*g2[i-1]

            if (abs(g1[i]) > 999999999): dG1 = 0
            if (abs(g2[i]) > 999999999): dG2 = 0
            
            g1[i] = g1[i-1] + dG1 * (t[i] - t[i-1])
            g2[i] = g2[i-1] + dG2 * (t[i] - t[i-1])
            g[i] = g1[i] + g2[i]
        
        return g

##### End of time variable concentrated two state model class definition. 
Revision:	25
Committed:	Wed Apr 28 20:22:24 2010 UTC (9 years, 3 months ago) by rjaynes
File size:	4751 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 TWO STATE REACTION model with time variable concentration..
4	The reaction equations are:
5	L + A ==ka1=> LA
6	L + A <=kd1== LA
7	LA ==ka2=> LA*
8	LA <=kd2== LA*
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 modeltwostate_varyC as mtvc
15	>m = mtvc.twostatemodel_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
29	class twostatemodel_varyC(modelclass):
30	"""The data attributes of the class:
31	parainfo: the parameter information
32	sim_data: the simulated data
33	conc: the time variable concentration data
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	def wizard(self, ):
65	'''This function helps you to create a parameter information list.
66	This is a competing reactions model, where there are 6 parameters:
67	rmax: Maximum analyte binding capacity(RU),
68	ka1: Association rate constant for L+A=LA(M-1S-1),
69	kd1: Dissociation rate constant for LA=L+A(S-1),
70	ka2: Forward rate constant for LA=LA*(S-1),
71	kd2: Backward rate constant for LA*=LA(S-1),
72	ca: Analyte concentration(M).'''
73	print ('')
74	pname = ['rmax','ka1','kd1','ka2','kd2','ca']
75	parainfo = []
76	for i in range(6):
77	print i
78	tmp={}
79	tmp['name'] = pname[i]
80	vstr = raw_input("input the value of parameter '%s', if a series, seperate with ',': " %pname[i])
81	vtmp = vstr.strip().split(',')
82	tmp['number'] = len(vtmp)
83	if len(vtmp) == 1:
84	vtmp = float(vtmp[0])
85	else:
86	vtmp = map(float, vtmp)
87	tmp['value'] = vtmp
88	tmp['fixed'] = bool(input("Is '%s' fixed? (1/0): " %pname[i]))
89
90	parainfo.append(tmp)
91	##parainfo.append({'name':'', 'value':0., 'fixed':0, 'number':0})
92	self.parainfo = parainfo
93
94
95	def function(self, t, paralist):
96	'''This function calculates the theoretical curve through the
97	parameter list you give.
98	'''
99	# for two state reaction model
100	conc = copy.deepcopy(self.conc)
101	C = np.interp(t,conc[0],conc[1])
102	## Assign the parameters to calculate the curve
103	for p in paralist:
104	if p['name'] == 'rmax': rmax = p['value']
105	elif p['name'] == 'ka1': ka1 = p['value']
106	elif p['name'] == 'kd1': kd1 = p['value']
107	elif p['name'] == 'ka2': ka2 = p['value']
108	elif p['name'] == 'kd2': kd2 = p['value']
109	elif p['name'] == 'ca': ca = p['value']
110	else: print p['name'], p['value']
111	if type(ca1) == list or type(ca2) == list:
112	print "Error: This function can only generate data for a single concentration."
113	return
114
115	## Must iterate through data, numerical integration.
116	g = np.zeros(len(C), dtype=float)
117	g1 = np.zeros(len(C), dtype=float)
118	g2 = np.zeros(len(C), dtype=float)
119
120	for i in range(1,len(C)):
121	dG1 = (ka1caC[i-1](rmax-g1[i-1]-g2[i-1]) - kd1g1[i-1])-(ka2g1[i-1]- kd2g2[i-1])
122	dG2 = ka2g1[i-1]- kd2g2[i-1]
123
124	if (abs(g1[i]) > 999999999): dG1 = 0
125	if (abs(g2[i]) > 999999999): dG2 = 0
126
127	g1[i] = g1[i-1] + dG1 * (t[i] - t[i-1])
128	g2[i] = g2[i-1] + dG2 * (t[i] - t[i-1])
129	g[i] = g1[i] + g2[i]
130
131	return g
132
133	##### End of time variable concentrated two state model class definition.