@@ -49,16 +49,12 @@ def augmented_system(Y, t, p):
49
49
# I used CAS to get these derivatives
50
50
y0_sensitivity = np .exp (- a * t )
51
51
a_sensitivity = (
52
- - (np .exp (t * (a - 1 )) - 1 + (a - 1 ) * (y0 * a - y0 - 1 ) * t )
53
- * np .exp (- a * t )
54
- / (a - 1 ) ** 2
52
+ - (np .exp (t * (a - 1 )) - 1 + (a - 1 ) * (y0 * a - y0 - 1 ) * t ) * np .exp (- a * t ) / (a - 1 ) ** 2
55
53
)
56
54
57
55
sensitivity = np .c_ [y0_sensitivity , a_sensitivity ]
58
56
59
- integrated_solutions = odeint (
60
- func = augmented_system , y0 = [y0 , 1 , 0 ], t = t .ravel (), args = (p ,)
61
- )
57
+ integrated_solutions = odeint (func = augmented_system , y0 = [y0 , 1 , 0 ], t = t .ravel (), args = (p ,))
62
58
simulated_sensitivity = integrated_solutions [:, 1 :]
63
59
64
60
np .testing .assert_allclose (sensitivity , simulated_sensitivity , rtol = 1e-5 )
@@ -78,14 +74,12 @@ def ode_func(y, t, p):
78
74
y = 1.0 / (a - 1 ) * (np .exp (- t ) - np .exp (- a * t ))
79
75
80
76
# Instantiate ODE model
81
- ode_model = DifferentialEquation (
82
- func = ode_func , t0 = 0 , times = t , n_states = 1 , n_theta = 1
83
- )
77
+ ode_model = DifferentialEquation (func = ode_func , t0 = 0 , times = t , n_states = 1 , n_theta = 1 )
84
78
85
79
simulated_y , sens = ode_model ._simulate ([y0 ], [a ])
86
-
87
- assert simulated_y .shape == (len (t ),1 )
88
- assert sens .shape == (len (t ), 1 , 1 + 1 )
80
+
81
+ assert simulated_y .shape == (len (t ), 1 )
82
+ assert sens .shape == (len (t ), 1 , 1 + 1 )
89
83
np .testing .assert_allclose (y , simulated_y , rtol = 1e-5 )
90
84
91
85
@@ -102,9 +96,7 @@ def ode_func_1(y, t, p):
102
96
103
97
# Instantiate ODE model
104
98
# Instantiate ODE model
105
- model1 = DifferentialEquation (
106
- func = ode_func_1 , t0 = 0 , times = self .t , n_states = 1 , n_theta = 1
107
- )
99
+ model1 = DifferentialEquation (func = ode_func_1 , t0 = 0 , times = self .t , n_states = 1 , n_theta = 1 )
108
100
109
101
# Sensitivity initial condition for this model should be 1 by 2
110
102
model1_sens_ic = np .array ([1 , 0 ])
@@ -117,9 +109,7 @@ def ode_func_2(y, t, p):
117
109
return p [0 ] * np .exp (- p [0 ] * t ) - p [1 ] * y [0 ]
118
110
119
111
# Instantiate ODE model
120
- model2 = DifferentialEquation (
121
- func = ode_func_2 , t0 = 0 , times = self .t , n_states = 1 , n_theta = 2
122
- )
112
+ model2 = DifferentialEquation (func = ode_func_2 , t0 = 0 , times = self .t , n_states = 1 , n_theta = 2 )
123
113
124
114
model2_sens_ic = np .array ([1 , 0 , 0 ])
125
115
@@ -134,14 +124,9 @@ def ode_func_3(y, t, p):
134
124
return [ds , di ]
135
125
136
126
# Instantiate ODE model
137
- model3 = DifferentialEquation (
138
- func = ode_func_3 , t0 = 0 , times = self .t , n_states = 2 , n_theta = 1
139
- )
127
+ model3 = DifferentialEquation (func = ode_func_3 , t0 = 0 , times = self .t , n_states = 2 , n_theta = 1 )
140
128
141
- model3_sens_ic = np .array ([
142
- 1 , 0 , 0 ,
143
- 0 , 1 , 0
144
- ])
129
+ model3_sens_ic = np .array ([1 , 0 , 0 , 0 , 1 , 0 ])
145
130
146
131
np .testing .assert_array_equal (model3_sens_ic , model3 ._sens_ic )
147
132
@@ -154,14 +139,9 @@ def ode_func_4(y, t, p):
154
139
return [ds , di ]
155
140
156
141
# Instantiate ODE model
157
- model4 = DifferentialEquation (
158
- func = ode_func_4 , t0 = 0 , times = self .t , n_states = 2 , n_theta = 2
159
- )
142
+ model4 = DifferentialEquation (func = ode_func_4 , t0 = 0 , times = self .t , n_states = 2 , n_theta = 2 )
160
143
161
- model4_sens_ic = np .array ([
162
- 1 , 0 , 0 , 0 ,
163
- 0 , 1 , 0 , 0
164
- ])
144
+ model4_sens_ic = np .array ([1 , 0 , 0 , 0 , 0 , 1 , 0 , 0 ])
165
145
166
146
np .testing .assert_array_equal (model4_sens_ic , model4 ._sens_ic )
167
147
@@ -175,18 +155,12 @@ def ode_func_5(y, t, p):
175
155
return [dx , ds , dz ]
176
156
177
157
# Instantiate ODE model
178
- model5 = DifferentialEquation (
179
- func = ode_func_5 , t0 = 0 , times = self .t , n_states = 3 , n_theta = 3
180
- )
158
+ model5 = DifferentialEquation (func = ode_func_5 , t0 = 0 , times = self .t , n_states = 3 , n_theta = 3 )
181
159
182
160
# First three columns are derivatives with respect to ode parameters
183
161
# Last three coluimns are derivatives with repsect to initial condition
184
162
# So identity matrix should appear in last 3 columns
185
- model5_sens_ic = np .array ([
186
- [1 , 0 , 0 , 0 , 0 , 0 ],
187
- [0 , 1 , 0 , 0 , 0 , 0 ],
188
- [0 , 0 , 1 , 0 , 0 , 0 ]
189
- ])
163
+ model5_sens_ic = np .array ([[1 , 0 , 0 , 0 , 0 , 0 ], [0 , 1 , 0 , 0 , 0 , 0 ], [0 , 0 , 1 , 0 , 0 , 0 ]])
190
164
191
165
np .testing .assert_array_equal (np .ravel (model5_sens_ic ), model5 ._sens_ic )
192
166
@@ -203,25 +177,18 @@ def system_1(y, t, p):
203
177
y0 = 0.0
204
178
times = np .arange (0.5 , 8 , 0.5 )
205
179
206
- yobs = np .array ([
207
- 0.30 , 0.56 , 0.51 , 0.55 , 0.47 , 0.42 , 0.38 , 0.30 ,
208
- 0.26 , 0.21 , 0.22 , 0.13 , 0.13 , 0.09 , 0.09
209
- ])[:,np .newaxis ]
180
+ yobs = np .array (
181
+ [0.30 , 0.56 , 0.51 , 0.55 , 0.47 , 0.42 , 0.38 , 0.30 , 0.26 , 0.21 , 0.22 , 0.13 , 0.13 , 0.09 , 0.09 ]
182
+ )[:, np .newaxis ]
210
183
211
- ode_model = DifferentialEquation (
212
- func = system_1 , t0 = 0 , times = times , n_theta = 1 , n_states = 1
213
- )
184
+ ode_model = DifferentialEquation (func = system_1 , t0 = 0 , times = times , n_theta = 1 , n_states = 1 )
214
185
215
186
integrated_solution , * _ = ode_model ._simulate ([y0 ], [alpha ])
216
187
217
188
assert integrated_solution .shape == yobs .shape
218
189
219
190
# compare automatic and manual logp values
220
- manual_logp = norm .logpdf (
221
- x = np .ravel (yobs ),
222
- loc = np .ravel (integrated_solution ),
223
- scale = 1
224
- ).sum ()
191
+ manual_logp = norm .logpdf (x = np .ravel (yobs ), loc = np .ravel (integrated_solution ), scale = 1 ).sum ()
225
192
with pm .Model () as model_1 :
226
193
forward = ode_model (theta = [alpha ], y0 = [y0 ])
227
194
y = pm .Normal ("y" , mu = forward , sd = 1 , observed = yobs )
@@ -238,9 +205,7 @@ def system(y, t, p):
238
205
239
206
times = np .arange (0 , 9 )
240
207
241
- ode_model = DifferentialEquation (
242
- func = system , t0 = 0 , times = times , n_states = 1 , n_theta = 1
243
- )
208
+ ode_model = DifferentialEquation (func = system , t0 = 0 , times = times , n_states = 1 , n_theta = 1 )
244
209
245
210
@pytest .mark .xfail (condition = (theano .config .floatX == "float32" ), reason = "Fails on float32" )
246
211
def test_too_many_params (self ):
@@ -264,21 +229,15 @@ def test_too_few_y0(self):
264
229
265
230
def test_func_callable (self ):
266
231
with pytest .raises (ValueError ):
267
- DifferentialEquation (
268
- func = 1 , t0 = 0 , times = self .times , n_states = 1 , n_theta = 1
269
- )
232
+ DifferentialEquation (func = 1 , t0 = 0 , times = self .times , n_states = 1 , n_theta = 1 )
270
233
271
234
def test_number_of_states (self ):
272
235
with pytest .raises (ValueError ):
273
- DifferentialEquation (
274
- func = self .system , t0 = 0 , times = self .times , n_states = 0 , n_theta = 1
275
- )
236
+ DifferentialEquation (func = self .system , t0 = 0 , times = self .times , n_states = 0 , n_theta = 1 )
276
237
277
238
def test_number_of_params (self ):
278
239
with pytest .raises (ValueError ):
279
- DifferentialEquation (
280
- func = self .system , t0 = 0 , times = self .times , n_states = 1 , n_theta = 0
281
- )
240
+ DifferentialEquation (func = self .system , t0 = 0 , times = self .times , n_states = 1 , n_theta = 0 )
282
241
283
242
284
243
class TestDiffEqModel :
@@ -294,7 +253,9 @@ def ode_func(y, t, p):
294
253
# Instantiate two Ops
295
254
op_1 = DifferentialEquation (func = ode_func , t0 = 0 , times = t , n_states = 1 , n_theta = 1 )
296
255
op_2 = DifferentialEquation (func = ode_func , t0 = 0 , times = t , n_states = 1 , n_theta = 1 )
297
- op_other = DifferentialEquation (func = ode_func , t0 = 0 , times = np .linspace (0 , 2 , 16 ), n_states = 1 , n_theta = 1 )
256
+ op_other = DifferentialEquation (
257
+ func = ode_func , t0 = 0 , times = np .linspace (0 , 2 , 16 ), n_states = 1 , n_theta = 1
258
+ )
298
259
299
260
assert op_1 == op_2
300
261
assert op_1 != op_other
@@ -310,14 +271,11 @@ def system(y, t, p):
310
271
[0.5 , 1.0 , 1.5 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 , 4.5 , 5.0 , 5.5 , 6.0 , 6.5 , 7.0 , 7.5 ]
311
272
)
312
273
313
- yobs = np .array ([
314
- 0.31 , 0.57 , 0.51 , 0.55 , 0.47 , 0.42 , 0.38 , 0.3 ,
315
- 0.26 , 0.22 , 0.22 , 0.14 , 0.14 , 0.09 , 0.1
316
- ])[:,np .newaxis ]
274
+ yobs = np .array (
275
+ [0.31 , 0.57 , 0.51 , 0.55 , 0.47 , 0.42 , 0.38 , 0.3 , 0.26 , 0.22 , 0.22 , 0.14 , 0.14 , 0.09 , 0.1 ]
276
+ )[:, np .newaxis ]
317
277
318
- ode_model = DifferentialEquation (
319
- func = system , t0 = 0 , times = times , n_states = 1 , n_theta = 1
320
- )
278
+ ode_model = DifferentialEquation (func = system , t0 = 0 , times = times , n_states = 1 , n_theta = 1 )
321
279
322
280
with pm .Model () as model :
323
281
alpha = pm .HalfCauchy ("alpha" , 1 )
@@ -341,14 +299,11 @@ def system(y, t, p):
341
299
[0.5 , 1.0 , 1.5 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 , 4.5 , 5.0 , 5.5 , 6.0 , 6.5 , 7.0 , 7.5 ]
342
300
)
343
301
344
- yobs = np .array ([
345
- 0.31 , 0.57 , 0.51 , 0.55 , 0.47 , 0.42 , 0.38 , 0.3 ,
346
- 0.26 , 0.22 , 0.22 , 0.14 , 0.14 , 0.09 , 0.1
347
- ])[:,np .newaxis ]
302
+ yobs = np .array (
303
+ [0.31 , 0.57 , 0.51 , 0.55 , 0.47 , 0.42 , 0.38 , 0.3 , 0.26 , 0.22 , 0.22 , 0.14 , 0.14 , 0.09 , 0.1 ]
304
+ )[:, np .newaxis ]
348
305
349
- ode_model = DifferentialEquation (
350
- func = system , t0 = 0 , times = times , n_states = 1 , n_theta = 2
351
- )
306
+ ode_model = DifferentialEquation (func = system , t0 = 0 , times = times , n_states = 1 , n_theta = 2 )
352
307
353
308
with pm .Model () as model :
354
309
alpha = pm .HalfCauchy ("alpha" , 1 )
@@ -375,24 +330,24 @@ def system(y, t, p):
375
330
376
331
times = np .array ([0.0 , 0.8 , 1.6 , 2.4 , 3.2 , 4.0 , 4.8 , 5.6 , 6.4 , 7.2 , 8.0 ])
377
332
378
- yobs = np .array ([
379
- [1.02 , 0.02 ],
380
- [0.86 , 0.12 ],
381
- [0.43 , 0.37 ],
382
- [0.14 , 0.42 ],
383
- [0.05 , 0.43 ],
384
- [0.03 , 0.14 ],
385
- [0.02 , 0.08 ],
386
- [0.02 , 0.04 ],
387
- [0.02 , 0.01 ],
388
- [0.02 , 0.01 ],
389
- [0.02 , 0.01 ],
390
- ])
391
-
392
- ode_model = DifferentialEquation (
393
- func = system , t0 = 0 , times = times , n_states = 2 , n_theta = 1
333
+ yobs = np .array (
334
+ [
335
+ [1.02 , 0.02 ],
336
+ [0.86 , 0.12 ],
337
+ [0.43 , 0.37 ],
338
+ [0.14 , 0.42 ],
339
+ [0.05 , 0.43 ],
340
+ [0.03 , 0.14 ],
341
+ [0.02 , 0.08 ],
342
+ [0.02 , 0.04 ],
343
+ [0.02 , 0.01 ],
344
+ [0.02 , 0.01 ],
345
+ [0.02 , 0.01 ],
346
+ ]
394
347
)
395
348
349
+ ode_model = DifferentialEquation (func = system , t0 = 0 , times = times , n_states = 2 , n_theta = 1 )
350
+
396
351
with pm .Model () as model :
397
352
R = pm .Lognormal ("R" , 1 , 5 )
398
353
sigma = pm .HalfCauchy ("sigma" , 1 , shape = 2 )
@@ -413,25 +368,25 @@ def system(y, t, p):
413
368
return [ds , di ]
414
369
415
370
times = np .array ([0.0 , 0.8 , 1.6 , 2.4 , 3.2 , 4.0 , 4.8 , 5.6 , 6.4 , 7.2 , 8.0 ])
416
-
417
- yobs = np .array ([
418
- [1.02 , 0.02 ],
419
- [0.86 , 0.12 ],
420
- [0.43 , 0.37 ],
421
- [0.14 , 0.42 ],
422
- [0.05 , 0.43 ],
423
- [0.03 , 0.14 ],
424
- [0.02 , 0.08 ],
425
- [0.02 , 0.04 ],
426
- [0.02 , 0.01 ],
427
- [0.02 , 0.01 ],
428
- [0.02 , 0.01 ],
429
- ])
430
-
431
- ode_model = DifferentialEquation (
432
- func = system , t0 = 0 , times = times , n_states = 2 , n_theta = 2
371
+
372
+ yobs = np .array (
373
+ [
374
+ [1.02 , 0.02 ],
375
+ [0.86 , 0.12 ],
376
+ [0.43 , 0.37 ],
377
+ [0.14 , 0.42 ],
378
+ [0.05 , 0.43 ],
379
+ [0.03 , 0.14 ],
380
+ [0.02 , 0.08 ],
381
+ [0.02 , 0.04 ],
382
+ [0.02 , 0.01 ],
383
+ [0.02 , 0.01 ],
384
+ [0.02 , 0.01 ],
385
+ ]
433
386
)
434
387
388
+ ode_model = DifferentialEquation (func = system , t0 = 0 , times = times , n_states = 2 , n_theta = 2 )
389
+
435
390
with pm .Model () as model :
436
391
beta = pm .HalfCauchy ("beta" , 1 )
437
392
gamma = pm .HalfCauchy ("gamma" , 1 )
0 commit comments