ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/mdl_module.py
Revision: 69
Committed: Wed May 16 18:00:28 2012 UTC (7 years, 3 months ago) by clausted
File size: 26151 byte(s)
Log Message:
Change mdl_module.py range() to xrange() to increase speeds 5%.
Line File contents
1 """
2 mdl_module
3 -------------------------------------------------------------------------------
4
5 *Example model functions module for SPRI data.*
6 Just a few simple, common interaction models will go here.
7 In the future, each new model will go in its own module/file.
8 Each model is a function taking three parameters.
9
10 **Parameters:**
11 * *time* (numpy array) -- Time points, usually in seconds.
12 * *data* (numpy array) -- Initial SPR signal values.
13 * *params* (dictionary of dictionaries) -- Model parameter description.
14
15 **Returns:** numpy array of calculated SPR signal values.
16
17 The dictionary keys are the names of model parameters
18 (e.g. 'rmax' for maximal response, or 'kon' for kinetic on-rate).
19 Each model parameter is described by a subdictionary containing four entries.
20
21 * ``'value'`` The current value of the parameter.
22 * ``'min'`` The minimum allowable value.
23 * ``'max'`` The maximum allowable value.
24 * ``'fixed'`` Either 'float', 'fixed', or a reference to another ROI.
25
26 The ``min``, ``max``, and ``fixed`` keys are used during automatic
27 curve-fitting. A fixed parameter is not allowed to change, while a
28 float parameter is adjusted until the least-squares algorithm has
29 minimized the sum-squared error. The *fixed* parameter may also be an
30 integer, in which case it is fixed to the value of a parameter of the
31 same name in another ROI.
32
33 The model function returns an array of values obtained by numerical
34 integration. The model is represented by differential equations and
35 integrated using the rectangle rule or, preferentially, using the
36 trapezoidal rule.
37
38 Functions are also available that provide parameter dictionaries
39 containing all the necessary default values.
40
41 .. moduleauthor:: Christopher Lausted,
42 Institute for Systems Biology,
43 OSPRAI developers.
44
45 Examples::
46
47 >>> import mdl_module as mdl, numpy as np
48 >>> times, data = np.arange(300), np.zeros(300)
49 >>>
50 >>> param1 = {'rate': {'value':1.0, 'min':-10.0, 'max':10.0, 'fixed':True} }
51 >>> data1 = mdl.drift(times, data, param1)
52 >>> print round(data1[299], 1)
53 299.0
54 >>> param1 = dict(rate=dict(value=-1.0, min=-10.0, max=10.0, fixed=True))
55 >>> data1 = mdl.drift(times, data, param1)
56 >>> print round(data1[299], 1)
57 -299.0
58 >>>
59 >>> param2 = mdl.simple1to1_def_params()
60 >>> param2['rmax'] ['value'] = 1000.0
61 >>> param2['conc']['value'] = 1e-6
62 >>> param2['t1']['value'] = 100
63 >>> param2['kon']['value'] = 1e4
64 >>> param2['t2']['value'] = 200.0
65 >>> param2['koff']['value'] = 1e-2
66 >>> param2['t3']['value'] = 300.0
67 >>> data2 = mdl.simple1to1(times, data, param2)
68 >>> print round(data2[299], 1)
69 160.3
70
71 """
72 __version__ = "120515"
73
74
75 ## Import libraries
76 import sys
77 import numpy as np
78 from copy import deepcopy
79 from math import sqrt
80
81
82 ############ Classes needed to support the modeling. #################
83
84 class LinearRegression(object):
85 """
86 Linear Least-squares regression class.
87 Example::
88
89 >>> r = LinearRegression()
90 >>> r.data(0, 0)
91 >>> r.data(10, 23)
92 >>> r.data(20, 40)
93 >>> print r.slope(), r.intercept()
94 2.0 1.0
95 >>> print r.yvalues()
96 [1.0, 21.0, 41.0]
97
98 """
99 def __init__(self):
100 """Create a new regression object containin no data."""
101 self.sumx = 0.0 ## Sum of x values.
102 self.sumy = 0.0 ## Sum of x values.
103 self.sumxx = 0.0 ## Sum of x**2 values.
104 self.sumxy = 0.0 ## Sum of x*y values.
105 self.sumyy = 0.0 ## Sum of y**2 values.
106 self.xlist = [] ## Keep a list of x values.
107 return
108 def data(self, x, y):
109 """Add a data point by providing an x and y value."""
110 self.sumx += x
111 self.sumy += y
112 self.sumxx += x*x
113 self.sumxy += x*y
114 self.sumyy += y*y
115 self.xlist.append(float(x))
116 return
117 def slope(self):
118 """Return the slope from the linear regression."""
119 n = len(self.xlist)
120 if (n<2): raise ZeroDivisionError("Regression requires >2 points.")
121 det = n*self.sumxx - self.sumx**2
122 return ((n*self.sumxy - self.sumx*self.sumy) / det)
123 def intercept(self):
124 """Return the y-intercept from the linear regression."""
125 n = len(self.xlist)
126 if (n<2): raise ZeroDivisionError("Regression requires >2 points.")
127 det = n*self.sumxx - self.sumx**2
128 return ((self.sumxx*self.sumy - self.sumx*self.sumxy) / det)
129 def xvalues(self):
130 """Return the list of x values provided."""
131 return self.xlist
132 def yvalues(self):
133 """Return the calculated y values for the given x values."""
134 m,b = self.slope(), self.intercept()
135 return [(m*x+b) for x in self.xlist]
136 """End of LinearRegression class."""
137
138
139 ############ The simple example models follow below. #################
140
141 def noise(time, data, params):
142 """
143 This function simply adds gaussian noise.
144 It requires numpy arrays of times and starting data values,
145 It only requires one parameter in the params list,
146 ``params['rms']['value']``.
147 """
148 noise = float(params['rms']['value'])
149 return (noise * np.random.randn(len(data)) + data)
150
151
152 def drift(time, data, params):
153 """
154 This function simply models a constant signal drift in units/second.
155 It requires numpy arrays of times and starting data values.
156 It only requires one parameter in the params list,
157 ``params['rate']['value']``.
158 """
159 y = np.zeros(len(time), dtype=float)
160 try:
161 rate = params['rate']['value']
162 except KeyError:
163 print "Error: Rate parameter not supplied to drift() model."
164 return y
165
166 ## Numerical integration.
167 y[0] = data[0]
168 for i in xrange(1,len(time)):
169 dy = rate * (time[i]-time[i-1])
170 y[i] = y[i-1] + dy
171
172 return y
173 ## End of drift() function.
174
175 def drift_def_params():
176 """
177 Return dictionary of default parameters for drift model.
178 For example::
179
180 >>> from numpy import arange, zeros
181 >>> defaultparam0 = drift_def_params()
182 >>> data0 = drift(arange(100), zeros(100), defaultparam0)
183 >>> print round(data0[99], 1)
184 99.0
185 """
186 dpar = dict(rate=dict(value=1.0, min=-100.0, max=100.0, fixed=True))
187 return dpar
188
189
190 def simple1to1(time, data, params):
191 """
192 This function simply models a 1:1 interaction.
193 The model parameters are described in the table below.
194
195 Model::
196
197 [A] + [L] --kon--> [AL]
198 [A] + [L] <-koff-- [AL]
199
200 Derivation::
201
202 d[AL]/dt = kon*[A]*[L] - koff*[AL]
203 y = [AL]
204 (rmax-y) = [L]
205 conc = [A]
206 rmax = [AL] + [L]
207 dy/dt = conc*kon*(rmax-y) - koff*y
208
209 ======================= ============================================
210 Model Parameter Description
211 ======================= ============================================
212 params['t0']['value'] time of start baseline stabilization, (s)
213 params['t1']['value'] time of injection for binding, (s)
214 params['rmax']['value'] maximum response, (RIU)
215 params['conc']['value'] concentration of analyte [A], (M)
216 params['kon']['value'] on-rate of analyte, (1/Ms)
217 params['t2']['value'] time of end binding & begin washing, (s)
218 params['koff']['value'] off-rate of analyte, (1/s)
219 params['t3']['value'] time end of washing & data fitting, (s)
220 ----------------------- --------------------------------------------
221 *Optional*
222 params['cofa']['value'] concentration factor, (1/dilution factor)
223 ======================= ============================================
224 """
225
226 ## Skip parameter validation steps for now.
227 t0 = float(params['t0']['value'])
228 t1 = float(params['t1']['value'])
229 rmax = float(params['rmax']['value'])
230 conc = float(params['conc']['value'])
231 kon = float(params['kon']['value'])
232 t2 = float(params['t2']['value'])
233 koff = float(params['koff']['value'])
234 t3 = float(params['t3']['value'])
235 if ('cofa' in params.keys()):
236 conc *= float(params['cofa']['value'])
237
238 ## Initialize variables.
239 errlog = {} ## Error status dictionary.
240 y = np.zeros(len(time), dtype=float)
241 lreg = LinearRegression()
242 base = y[0] = data[0] ## Baseline SPR signal.
243
244 ## Must iterate through data, numerical integration.
245 for i in xrange(1,len(time)):
246 dt = time[i] - time[i-1]
247 ## Treat function as having three separate phases.
248
249 ## 1. If in pre-binding, baseline stabilization phase.
250 if (t0 < time[i] <= t1):
251 ## Append this data point to the linear regression.
252 lreg.data(time[i], data[i])
253 ## And if this is the last data point of pre-binding.
254 if (time[i+1] > t1):
255 ## Fit a line to the baseline data.
256 yreg = lreg.yvalues()
257 n = len(yreg)
258 ## Copy it in to the model and reset.
259 y[(1+i-n):(1+i)] = yreg
260 base = yreg[-1]
261 lreg = LinearRegression()
262
263 ## 2. If in binding phase, sample injection in progress.
264 elif (t1 < time[i] <= t2):
265 yb = y[i-1] - base
266 dy = conc*kon*(rmax-yb) - koff*yb
267 y[i] = y[i-1] + dy*dt
268 ## Check if we overshot the maximum response.
269 if (y[i] > (base+rmax)):
270 y[i] = base+rmax
271 errlog['Response overshot rmax.'] = True
272
273 ## 3. If in dissociation phase, conc=0.
274 elif (t2 < time[i] <= t3):
275 yb = y[i-1] - base
276 dy = 0 - koff*yb
277 y[i] = y[i-1] + dy*dt
278
279 ## End "for i" loop.
280
281 if any(errlog): print "Errors in simple1to1:", errlog.keys()
282
283 return y
284 """End of simple1to1() function"""
285
286 def simple1to1_def_params():
287 """
288 Return a dictionary of default parameters for simple1to1 model.
289 For example::
290
291 >>> from numpy import arange, zeros
292 >>> defaultparam0 = simple1to1_def_params()
293 >>> data0 = simple1to1(arange(300), zeros(300), defaultparam0)
294 >>> print round(data0[299], 1)
295 160.3
296 """
297 dpar = {}
298 dpar['t0'] = dict(value=0.0, min=0.0, max=0.0, fixed='fixed')
299 dpar['rmax'] = dict(value=1e3, min=1e1, max=1e4, fixed='fixed')
300 dpar['conc'] = dict(value=1e-6, min=1e-12, max=1e-1, fixed='fixed')
301 dpar['cofa'] = dict(value=1.0, min=1e-6, max=1e6, fixed='fixed')
302 dpar['t1'] = dict(value=100.0, min=100.0, max=100.0, fixed='fixed')
303 dpar['kon'] = dict(value=1e4, min=1e2, max=1e7, fixed='fixed')
304 dpar['t2'] = dict(value=200.0, min=200.0, max=200.0, fixed='fixed')
305 dpar['koff'] = dict(value=1e-2, min=1e-6, max=1e-1, fixed='fixed')
306 dpar['t3'] = dict(value=300.0, min=300.0, max=300.0, fixed='fixed')
307 return dpar
308
309
310 def simple1to1_mtl(time, data, params):
311 """
312 This function simply models a 1:1 interaction with mass transport limitation.
313 The model parameters are described in the table below.
314
315 Model::
316
317 [Abulk] --km-> [Asurf] + [L] --kon--> [AL]
318 [Abulk] <-km-- [Asurf] + [L] <-koff-- [AL]
319
320 Derivation::
321
322 d[AL]/dt = (kon*[A]*[L] - koff*[AL]) / (1 + kon*[L]/kmtl)
323 y = [AL]
324 (rmax-y) = [L]
325 conc = [Abulk]
326 rmax = [AL] + [L]
327 dy/dt = (kon*conc*(rmax-y) - koff*y) / (1 + kon*(rmax-y)/kmtl)
328
329 ======================= ==========================================
330 Model Parameter Description
331 ======================= ==========================================
332 params['t0']['value'] time of start baseline stabilization, (s)
333 params['t1']['value'] time of injection for binding, (s)
334 params['rmax']['value'] maximum response, (RIU)
335 params['conc']['value'] concentration of analyte [Abulk], (M)
336 params['kon']['value'] on-rate of analyte, (1/Ms)
337 params['t2']['value'] time of end binding & begin washing, (s)
338 params['koff']['value'] off-rate of analyte, (1/s)
339 params['t3']['value'] time end of washing & data fitting, (s)
340 params['kmtl']['value'] rate of diffusion, (RIU/Ms)
341 ----------------------- ------------------------------------------
342 *Optional*
343 params['cofa']['value'] concentration factor, (1/dilution factor)
344 ======================= ==========================================
345 """
346
347 ## Skip parameter validation steps for now.
348 t0 = params['t0']['value']
349 t1 = params['t1']['value']
350 rmax = params['rmax']['value']
351 conc = params['conc']['value']
352 kon = params['kon']['value']
353 t2 = params['t2']['value']
354 koff = params['koff']['value']
355 t3 = params['t3']['value']
356 kmtl = params['kmtl']['value']
357 if ('cofa' in params.keys()):
358 conc *= float(params['cofa']['value'])
359
360 ## Initialize variables.
361 errlog = {} ## Error status dictionary.
362 y = np.zeros(len(time), dtype=float)
363 lreg = LinearRegression()
364 base = y[0] = data[0] ## Baseline SPR signal.
365
366 ## Error checks.
367 if (kmtl<=0):
368 ## Avoid div/0, assume user wants no mass transport limitation.
369 stat['kmtl must be > 0'] = True
370 kmtl = 1e40 ## An arbitrary very large number.
371
372 ## Must iterate through data, numerical integration.
373 for i in xrange(1,len(time)):
374 dt = time[i] - time[i-1]
375 ## Treat function as having three separate phases.
376
377 ## 1. If in pre-binding, baseline stabilization phase.
378 if (t0 < time[i] <= t1):
379 ## Append data to linear regression
380 lreg.data(time[i], data[i])
381 ## And if this the the last data point pre-binding.
382 if (time[i+1] > t1):
383 ## Fit a line to the baseline data.
384 yreg = lreg.yvalues()
385 n = len(yreg)
386 ## Copy it in to the model and reset.
387 y[(1+i-n):(1+i)] = yreg
388 base = yreg[-1]
389 lreg = LinearRegression()
390
391 ## 2. If in binding phase, sample injection in progress.
392 elif (t1 < time[i] <= t2):
393 ## Binding phase
394 yb = y[i-1] - base
395 dy = (conc*kon*(rmax-yb) - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
396 y[i] = y[i-1] + dy*dt
397 ## Check if we overshot the maximum response.
398 if (y[i] > (base+rmax)):
399 y[i] = base+rmax
400 errlog['Response overshot rmax.'] = True
401
402 ## 3. If in dissociation phase, conc=0.
403 elif (t2 < time[i] <= t3):
404 ## Dissociation (conc=0)
405 yb = y[i-1] - base
406 dy = (0 - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
407 y[i] = y[i-1] + dy*dt
408
409 if any(errlog): print "Errors in simple1to1_mtl:", errlog.keys()
410
411 return y
412 """End of simple1to1_mtl() function"""
413
414 def simple1to1_mtl_def_params():
415 """
416 Return dictionary of default parameters for simple1to1_mtl model.
417 For example::
418
419 >>> from numpy import arange, zeros
420 >>> defaultparam0 = simple1to1_mtl_def_params()
421 >>> data0 = simple1to1_mtl(arange(300), zeros(300), defaultparam0)
422 >>> print round(data0[299], 1)
423 161.1
424 """
425 dpar = simple1to1_def_params() ## Based on simple1to1.
426 dpar['kmtl'] = dict(value=1e9, min=1e2, max=1e10, fixed='fixed')
427 return dpar
428
429
430 def simple1to1_series(time, data, params):
431 """
432 Model up to 8 simple 1:1 SPR injections of varying concentration.
433
434 ======================= ==========================================
435 Model Parameter Description
436 ======================= ==========================================
437 params['kon']['value'] on-rate of analyte, (1/Ms)
438 params['koff']['value'] off-rate of analyte, (1/s)
439 params['rmax']['value'] maximum response, (RIU)
440 params['tp']['value'] time prior to binding, (s)
441 params['ta']['value'] time/duration of association, (s)
442 params['td']['value'] time/duration of dissociation, (s)
443 ----------------------- ------------------------------------------
444 params['t0']['value'] timepoint of first (no.0) injection, (s)
445 ... ...
446 params['t7']['value'] timepoint of eighth (no.7) injection, (s)
447 ----------------------- ------------------------------------------
448 params['c0']['value'] concentration of first injection, (s)
449 ... ...
450 params['c7']['value'] concentration of eighth injection, (s)
451 ======================= ==========================================
452
453 Injections where tN=0 are ignored. For example, if you have only
454 six injections, set t6=0 and t7=0.
455 """
456
457 ## Parameters from params dictionary.
458 ## Skip parameter validation steps for now.
459 kon = float(params['kon']['value'])
460 koff = float(params['koff']['value'])
461 rmax = float(params['rmax']['value'])
462 ## Pre-association, Association duration, Dissociation duration.
463 tp = float(params['tp']['value']) ## Typical -100 s
464 ta = float(params['ta']['value']) ## Typical 300 s
465 td = float(params['td']['value']) ## Typical 300 s
466 td = td + ta ## After 300s binding and 300s washing, 600s elapse.
467 ## How many injections are there? Check valid keys and inj times.
468 tkey = ['t0','t1','t2','t3','t4','t5','t6','t7']
469 ckey = ['c0','c1','c2','c3','c4','c5','c6','c7']
470 for i in range(8):
471 if (tkey[i] in params) and (ckey[i] in params):
472 ## Valid inj time is after sufficient pre-association time.
473 if ((float(params[tkey[i]]['value']) + tp) > 0):
474 inj = i+1 ## Increase inj count.
475 ## Up to eight injection start times, concentrations.
476 ti = [float(params[k]['value']) for k in tkey[0:inj]]
477 c = [float(params[k]['value']) for k in ckey[0:inj]]
478
479 ## Begin.
480 y = np.zeros(len(time), dtype=float)
481 lreg = LinearRegression()
482 base = y[0] = data[0] ## Baseline SPR signal.
483
484 ## Must iterate through data, numerical integration.
485 for i in xrange(1,len(time)):
486 dt = time[i] - time[i-1]
487 #if (time[i] > (max(ti)+td)): break ## Speed things up.
488 inmodel = False
489
490 ## For-loop where j is each injection number.
491 ## Treat function as having three separate phases.
492 for j in xrange(inj):
493 ## 1. If in pre-binding, baseline stabilization phase.
494 if ((ti[j]+tp) <= time[i] <= (ti[j])):
495 ## Append this datapoint to linear regression.
496 lreg.data(time[i], data[i])
497 ## And if this is last data point of pre-binding phase.
498 if (time[i+1] > ti[j]):
499 ## Fit a line to the baseline data.
500 yreg = lreg.yvalues()
501 n = len(yreg)
502 ## Copy it in to the model and reset.
503 y[(1+i-n):(1+i)] = yreg
504 base = yreg[-1]
505 lreg = LinearRegression()
506
507 ## 2. If in binding phase, sample injection in progress.
508 elif ((ti[j]) <= time[i] <= (ti[j]+ta)):
509 yb = y[i-1] - base
510 dy = c[j]*kon*(rmax-yb) - koff*yb
511 y[i] = y[i-1] + dy*dt
512 if (y[i] > (base+rmax)): y[i] = (base+rmax) ## Avoid overshoot.
513
514 ## 3. If in dissociation phase, conc=0.
515 elif ((ti[j]+ta) <= time[i] <= (ti[j]+td)):
516 yb = y[i-1]-base
517 dy = 0 - koff*yb
518 y[i] = y[i-1] + dy*dt
519
520 ## Check if we are anywhere in the j-th model.
521 if ((ti[j]+tp) <= time[i] <= (ti[j]+td)):
522 inmodel = True
523
524 ## End "for j" loop.
525
526 ## Handle datapoints outside the model, e.g. regeneration.
527 if (inmodel==False):
528 ## Should I use zero or np.NaN or other? 0.
529 y[i] = 0
530
531 ## End "for i" loop.
532 return y
533
534 def simple1to1_series_def_params():
535 """
536 Return a dictionary of default parameters for simple1to1_series model.
537 While parameters for eight injection times are provided, only the
538 first injection time is nonzero. For example::
539
540 >>> from numpy import arange, zeros
541 >>> defaultparam0 = simple1to1_series_def_params()
542 >>> data0 = simple1to1_series(arange(300), zeros(300), defaultparam0)
543 >>> print round(data0[299], 1)
544 160.3
545 """
546 dpar = {}
547 dpar['kon'] = dict(value=1e4, min=1e2, max=1e7, fixed='fixed')
548 dpar['koff'] = dict(value=1e-2, min=1e-6, max=1e-1, fixed='fixed')
549 dpar['rmax'] = dict(value=1e3, min=1e1, max=1e4, fixed='fixed')
550 dpar['tp'] = dict(value=-50.0, min=-50.0, max=-50.0, fixed='fixed')
551 dpar['ta'] = dict(value=100.0, min=100.0, max=100.0, fixed='fixed')
552 dpar['td'] = dict(value=100.0, min=100.0, max=100.0, fixed='fixed')
553 ## Set all eight concentrations to 1uM.
554 for key in ['c0','c1','c2','c3','c4','c5','c6','c7']:
555 dpar[key] = dict(value=1e-6, min=1e-12, max=1e-1, fixed='fixed')
556 ## Just describe the first injection at t0=100s.
557 for key in ['t0','t1','t2','t3','t4','t5','t6','t7']:
558 dpar[key] = dict(value=0.0, min=0.0, max=0.0, fixed='fixed')
559 dpar['t0'] = dict(value=100.0, min=100.0, max=100.0, fixed='fixed')
560 return dpar
561
562
563 ###################### Code testing ##################################
564
565 def _demo():
566 """
567 Test this module by producing some graphical output.
568 Simulated data is provided by ba_class.
569 Output is provided by vu_module.
570 Use both simple1to1() and simple1to1_series().
571 Report how long it takes to run the model and fit curves.
572 """
573 import osprai_one as osp
574 from time import time
575
576 timer = time()
577 ## Create an example 300s sensorgram with ba_class default values.
578 ## Defaults: kon=1e4, koff=1e-2, conc=1e-6, rmax=1000, start at 100s.
579 ba1 = osp.BiosensorArray(1, 300)
580 roi = ba1.roi[0] ## We will use only one roi.
581 roi.simulate(bulkRI=0.0, kon=1e4, koff=1e-2, conc=1e-6)
582 ## Add gaussian noise of 5 RU and a jump of -100 at t=10s.
583 noise = 5.0 * np.random.randn(300)
584 roi.value += noise
585 roi.value[10:300] -= 50
586
587 ## Definine model & parameters for simple1to1.
588 roi.model = simple1to1
589 roi.params = simple1to1_def_params()
590 #roi.model = simple1to1_mtl
591 #roi.params = simple1to1_mtl_def_params()
592 ## Set the known or fixed parameters.
593 roi.params['t0']['value'] = 10
594 roi.params['t1']['value'] = 100
595 roi.params['t2']['value'] = 200
596 roi.params['t3']['value'] = 299
597 roi.params['conc']['value'] = 1e-6
598 roi.params['rmax']['value'] = 1000
599 ## Setup initial guesses for floating paratmers.
600 roi.params['kon']['fixed'] = 'float'
601 roi.params['koff']['fixed'] = 'float'
602 roi.params['kon']['value'] = 1e5
603 roi.params['koff']['value'] = 1e-3
604 ## Perform curve-fitting and display.
605 osp.mclma([roi])
606 timer = time() - timer
607 print "Elapsed time was {0:0.1f} milliseconds".format(1000*timer)
608 osp.dualgraph(ba1, ba1, "Data fitted by simple1to1() model.")
609
610 timer = time()
611 ## Definine model & parameters for simple1to1_series.
612 roi.model = simple1to1_series
613 roi.params = simple1to1_series_def_params()
614 ## Set the known or fixed parameters.
615 roi.params['t0']['value'] = 100
616 roi.params['tp']['value'] = -90
617 roi.params['ta']['value'] = 100
618 roi.params['td']['value'] = 200
619 roi.params['c0']['value'] = 1e-6
620 roi.params['rmax']['value'] = 1000
621 ## Setup initial guesses for floating paratmers.
622 roi.params['kon']['fixed'] = 'float'
623 roi.params['koff']['fixed'] = 'float'
624 roi.params['kon']['value'] = 1e5
625 roi.params['koff']['value'] = 1e-3
626 ## Perform curve-fitting and display.
627 osp.mclma([roi])
628 timer = time() - timer
629 print "Elapsed time was {0:0.1f} milliseconds".format(1000*timer)
630 osp.dualgraph(ba1, ba1, "Data fitted by simple1to1_series() model.")
631
632 return
633
634 def _test():
635 """
636 Automatic Code testing with doctest.
637 Enable use of ellipses as wildcard for returned text.
638 Doctest runs the example code in the docstrings as well as the following.
639 We test that diff eq models here match the simulation in BioSensorArray.
640 Here are the tests::
641
642 >>> from numpy import arange, zeros
643 >>> dpar = simple1to1_series_def_params()
644 >>> values = simple1to1_series(arange(300), zeros(300), dpar)
645 >>> print round(values[298], -1) ## Round to 10s place.
646 160.0
647 >>> import ba_class as ba
648 >>> ba1 = ba.BiosensorArray(1, 300)
649 >>> ba1.roi[0].simulate(bulkRI=0.0, kon=dpar['kon']['value'], \
650 koff=dpar['koff']['value'], conc=dpar['c0']['value'])
651 >>> print round(ba1.roi[0].value[299], -1)
652 160.0
653 """
654 import doctest
655 doctest.testmod(optionflags=doctest.ELLIPSIS)
656 print "Code testing with doctest is complete."
657 return
658
659 if __name__ == "__main__":
660 """
661 Use 'python mdl_module.py' for code testing.
662 Use 'python mdl_module.py demo' to also show a graphic demo.
663 """
664 _test() ## Automatic Code testing with doctest.
665 if ("demo" in sys.argv):
666 _demo() ## A more interactive, graphical test.
667
668
669 ##################### End of module ##################################