ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/test.py
Revision: 65
Committed: Fri May 20 16:33:02 2011 UTC (8 years, 3 months ago) by clausted
File size: 8958 byte(s)
Log Message:
Added much testing to all modules using docstrings and doctest.  To help this, I also added two optional parameters to mclma(), verbose and maxiter.  For example, use "mclma(rois, verbose=False, maxiter=99)" to run the regression up to 99 iterations without echoing progress to the standard output.

Line File contents
1 """
2 Christopher Lausted, Institute for Systems Biology
3 Last modified on 110519 (yymmdd)
4 Scripts to test various OSPRAI components.
5 """
6
7 ## Set the PYTHONPATH to your preferred OSPRAI installation.
8 import sys
9 sys.path.append("/home/clausted/Dropbox/ISB/subversion/osprai/trunk")
10
11 ## Import libaries
12 import os
13 import numpy as np
14 ## Import Osprai libraries.
15 import ba_class
16 import vu_module
17 import io_module
18 import cal_module
19 import mdl_module
20 import fit_module
21 import osprai_one
22 ## Reload Osprai libraries in case code has changed.
23 reload(ba_class)
24 reload (vu_module)
25 reload(io_module)
26 reload(cal_module)
27 reload(mdl_module)
28 reload(fit_module)
29 reload(osprai_one)
30 ## Or just import this next one...
31 from osprai_one import *
32
33
34 def test_io_module():
35 print "Loading SPR data files..."
36 os.chdir("exampledata")
37 ba1 = readbiosensor("example-biosensor.txt")
38 ba1 = readicmtxt("example-icm.txt")
39 ba1 = readsprit("example-sprit.txt")
40 ba1 = applykey(ba1, "example-key.tsv")
41 print "Writing SPR data files..."
42 writesprit(ba1,"testwritesprit.txt")
43 writeclamp(ba1,"testwriteclamp.txt")
44 writebiosensor(ba1, "testwritebiosensor.txt")
45 os.chdir("..")
46 print "Done with test_io_module."
47 return
48
49 def test_cal_module_vu_module():
50 print "Loading an SPR data file..."
51 ba1 = readicmtxt("./exampledata/example-icm.txt")
52
53 print "Plotting the raw data file..."
54 dotgraph(ba1)
55
56 print "Plotting the calibrated data..."
57 ba2 = calibrate(ba1, ba1, 2850, 3120)
58 linegraph(ba2, "Calibrated SPR Data")
59 dualgraph(ba1, ba2, "Uncalibrated and Calibrated SPR Data")
60
61 print "Use ROI #16 as background, subtract, and plot..."
62 bgset(ba2, 16)
63 ba3 = bgsubt(ba2)
64 linegraph(ba3, "Calibrated Referenced SPR Data")
65
66 print "Test of scatterplot using time points 500, 1000, 1500, 2000"
67 scatterplot(ba3, 3600, 4200, 3600, 4700, "Scatterplot")
68
69 print "Cutting out an interval..."
70 ba3 = copyinterval(ba3, 2800, 6200)
71 linegraph(ba3)
72 print "Done with test test_cal_module_vu_module."
73 return
74
75
76 def test_mdl_module_fit_module():
77 print "Test the model module..."
78 ## Set some default parameters.
79 dpar = {'t1': {'value':30.0, 'min':30.0, 'max':30.0, 'fixed':'fixed'} }
80 dpar['rmax'] = {'value':100.0, 'min':1.0, 'max':1000.0, 'fixed':'fixed'}
81 dpar['conc'] = {'value':1e-5, 'min':1e-12, 'max':1e1, 'fixed':'fixed'}
82 dpar['cofa'] = {'value':1.0, 'min':1e-10, 'max':1e10, 'fixed':'fixed'}
83 dpar['kon'] = {'value':2e4, 'min':1e1, 'max':1e8, 'fixed':'fixed'}
84 dpar['t2'] = {'value':150.0, 'min':10.0, 'max':1000.0, 'fixed':'fixed'}
85 dpar['koff'] = {'value':1.0e-2, 'min':1e-8, 'max':1e-1, 'fixed':'fixed'}
86 dpar['t3'] = {'value':270.0, 'min':270.0, 'max':270.0, 'fixed':'fixed'}
87 ## Create ba with three ROIs and 300 datapoints.
88 ba0 = BiosensorArray(3,300)
89 for roi in ba0.roi:
90 roi.time = arange(300, dtype=float) ## Samples every 1 second.
91 roi.value = zeros(300, dtype=float) + 20 ## Baseline signal is 20 units.
92 roi.params = deepcopy(dpar)
93 roi.model = simple1to1
94 ## Set the dilution factors
95 ba0.roi[0].params['cofa']['value'] = 1.0 / 27
96 ba0.roi[1].params['cofa']['value'] = 1.0 / 9
97 ba0.roi[2].params['cofa']['value'] = 1.0 / 3
98 ## Simulate the data and plot it.
99 for roi in ba0.roi:
100 roi.value = roi.model(roi.time, roi.value, roi.params)
101 linegraph(ba0, "Simulation")
102
103 ## Constrain rmax belowtrue value to show what happens.
104 print "Testing the fitting..."
105 ba1 = deepcopy(ba0)
106 ## Set one kon, koff, rmax to float.
107 ba1.roi[0].params['koff']['fixed'] = 'float'
108 ba1.roi[0].params['kon']['fixed'] = 'float'
109 ba1.roi[0].params['rmax'] = {'value':10.0, 'min':1.0, 'max':80.0, 'fixed':'float'}
110 ba1.roi[1].params['koff']['fixed'] = 0
111 ba1.roi[1].params['kon']['fixed'] = 0
112 ba1.roi[1].params['rmax']['fixed'] = 0
113 ba1.roi[2].params['koff']['fixed'] = 0
114 ba1.roi[2].params['kon']['fixed'] = 0
115 ba1.roi[2].params['rmax']['fixed'] = 0
116 ## Fit and then plot.
117 mclma(ba1.roi)
118 for roi in ba1.roi:
119 roi.value =roi.model(roi.time, roi.value, roi.params)
120 dualgraph(ba0, ba1, "Fitted Data with an incorrect rmax constraint")
121 print "Done with test test_mdl_module_fit_module."
122 return
123
124
125 def demo_mtl():
126 print "Test the model module..."
127 ba0 = BiosensorArray(1,300) ## Create ba1 with one ROI and 300 datapoints.
128 roi0 = ba0.roi[0]
129 roi0.time = arange(300, dtype=float) ## Samples every 1 second.
130 roi0.value = zeros(300, dtype=float) + 20 ## Baseline signal is 20 units.
131 roi0.params = {'t1': {'value':30.0, 'min':30.0, 'max':30.0, 'fixed':True} }
132 roi0.params['rmax'] = {'value': 100.0}
133 roi0.params['conc'] = {'value': 1e-6}
134 roi0.params['kon'] = {'value': 5e4}
135 roi0.params['t2'] = {'value': 150.0}
136 roi0.params['koff'] = {'value': 5e-3}
137 roi0.params['t3'] = {'value': 270.0}
138 roi0.params['kmtl'] = {'value': 2e6}
139 roi0.name = "kmtl=2e6"
140 roi0.model = simple1to1_mtl
141 roi0.value =roi0.model(roi0.time, roi0.value, roi0.params)
142
143 ba1 = deepcopy(ba0)
144 roi0 = ba1.roi[0]
145 roi0.name = "No MTL"
146 roi0.model = simple1to1
147 roi0.value =roi0.model(roi0.time, roi0.value, roi0.params)
148 dualgraph(ba0, ba1, "With and without MTL")
149 print "Done with test demo_mtl."
150
151
152 def test_mdl_module_series():
153 print "Test the model module..."
154
155 ## Multiple injection model (ba0)
156 ## Create ba with two ROIs and 950 datapoints.
157 ba0 = BiosensorArray(2,950)
158 for roi in ba0.roi:
159 roi.time = arange(950, dtype=float) ## Samples every 1 second.
160 roi.value = zeros(950, dtype=float) + 20 ## Baseline signal is 20.
161 roi.params = simple1to1_series_def_params()
162 roi.model = simple1to1_series
163 ## Define ROI0
164 roi = ba0.roi[0]
165 roi.params['kon']['value'] = 1e4
166 roi.params['t0']['value'] = 100
167 roi.params['c0']['value'] = 1e-6
168 roi.params['t1']['value'] = 400
169 roi.params['c1']['value'] = 2e-6
170 roi.params['t2']['value'] = 700
171 roi.params['c2']['value'] = 4e-6
172 roi.value = roi.model(roi.time, roi.value, roi.params)
173 ## Define ROI1
174 roi = ba0.roi[1]
175 roi.params['kon']['value'] = 2e3
176 roi.params['t0']['value'] = 100
177 roi.params['c0']['value'] = 1e-6
178 roi.params['t1']['value'] = 400
179 roi.params['c1']['value'] = 2e-6
180 roi.params['t2']['value'] = 700
181 roi.params['c2']['value'] = 4e-6
182 roi.value = roi.model(roi.time, roi.value, roi.params)
183
184 ## Single-injection model (ba1)
185 ba1 = BiosensorArray(2,950)
186 for roi in ba1.roi:
187 roi.time = arange(950, dtype=float) ## Samples every 1 second.
188 roi.value = zeros(950, dtype=float) + 15 ## Baseline signal is 19.
189 roi.params = simple1to1_def_params()
190 roi.model = simple1to1
191 ## Define ROI0
192 roi = ba1.roi[0]
193 roi.params['kon']['value'] = 1e4
194 roi.params['conc']['value'] = 1e-6
195 roi.value = roi.model(roi.time, roi.value, roi.params)
196 ## Define ROI1
197 roi = ba1.roi[1]
198 roi.params['kon']['value'] = 2e3
199 roi.params['conc']['value'] = 1e-6
200 roi.value = roi.model(roi.time, roi.value, roi.params)
201
202 ## Display graph.
203 dualgraph(ba0, ba1, "test_mdl_module_series")
204 return
205
206
207 def test_noisy_fit():
208 print "test_noisy_fit"
209 ## Create a ba1 with a single ROI.
210 ba1 = BiosensorArray(1,300)
211 roi = ba1.roi[0]
212 ## Simple 1to1 binding curve.
213 ## Note KD = 2.2e-2 / 4.4e-4 = 5e-7, log10(5e-7) = -6.3
214 roi.simulate(bulkRI=0.0, kon=4.4e4, koff=2.2e-2, conc=7.7e-7)
215 ## Add 5 units of noise.
216 params = dict(rms=dict(value=5))
217 roi.value = noise(roi.time, roi.value, params)
218 ## Create a copy ba2 and use for fitting three parameters.
219 ba2 = deepcopy(ba1)
220 roi = ba2.roi[0]
221 roi.model = simple1to1
222 roi.params = simple1to1_def_params()
223 roi.params['rmax']['fixed'] = 'float'
224 roi.params['kon']['fixed'] = 'float'
225 roi.params['koff']['fixed'] = 'float'
226 roi.params['conc']['value'] = 7.7e-7
227 print roi.params
228 ## Fit.
229 #mclma(ba2.roi, maxiter=15, verbose=False)
230 mclma(ba2.roi)
231 roi.value = zeros(len(roi.value))
232 roi.value = roi.model(roi.time, roi.value, roi.params)
233 ## View
234 dualgraph(ba1, ba2, "Fit to noisy data")
235 kD = float(roi.params['koff']['value']) / float(roi.params['kon']['value'])
236 print "The expected value is -6.3. The fitted value is:"
237 print round((np.log10(kD)), 1)
238 return
239
240
241 ## Run selected test or tests.
242
243 #test_io_module()
244 #test_cal_module_vu_module()
245
246 #demo_mtl()
247 #test_mdl_module_fit_module()
248 #test_mdl_module_series()
249 test_noisy_fit()
250
251