File: ode.py

package info (click to toggle)
python-bumps 1.0.3-3
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 6,200 kB
  • sloc: python: 24,517; xml: 493; ansic: 373; makefile: 211; javascript: 99; sh: 94
file content (69 lines) | stat: -rw-r--r-- 1,764 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
# Fitting an ODE
# ==============
#
# Bumps can fit black-box functions, such as odeint from scipy.
#
# The following example is adapted from:
#
#     https://people.duke.edu/~ccc14/sta-663/CalibratingODEs.html.
#
# | Instructor: Cliburn Chan cliburn.chan@duke..edu
# | Instructor: Janice McCarthy janice.mccarthy@duke.edu

from bumps.names import *
import numpy as np
from scipy.integrate import odeint

# Define the ODE


def g(t, x0, a, b):
    """
    Solution to the ODE x'(t) = f(t,x,k) with initial condition x(0) = x0
    """
    return odeint(dfdt, x0, t, args=(a, b)).flatten()


def dfdt(x, t, a, b):
    """Receptor synthesis-internalization model."""
    return a - b * x


# Simulate some data.
#
# Note that the function :func:`bumps.util.push_seed` is to set the random
# number generator to a known state so that this function will create the
# same data every time the simulation is run.  If not, then you wouldn't
# be able to resume a fit since each time you resumed you would be fitting
# different data.


def simulate():
    from bumps.util import push_seed

    # Fake some data
    a = 2.0
    b = 0.5
    x0 = 10.0
    t = np.linspace(0, 10, 10)
    dy = 0.2 * np.ones_like(t)
    with push_seed(1):
        y = g(t, x0, a, b) + dy * np.random.normal(size=t.shape)
    # print(a, b, x0, t, dt, gt)
    return t, y, dy


t, y, dy = simulate()

# Define the fit problem.
#
# In this case :class:`bumps.curve.Curve` is initialized with *plot_x*
# as a vector of length 1000.  This is so that a smooth curve is drawn between
# the ten data points that were simulated in the fit.

M = Curve(g, t, y, dy, x0=1.0, a=1.0, b=1.0, plot_x=np.linspace(t[0], t[-1], 1000))
M.x0.range(0, 100)
M.a.range(0, 10)
M.b.range(0, 10)

problem = FitProblem(M)