Geek Answers Handbook

Efficient simulation of an arbitrary dynamic system in Python via numba?

For some time I wondered what is the most efficient way to model an arbitrary nonlinear (deterministic stochastic) dynamical system in Python. In the end, I do a lot of teaching or research. I am convinced that there must be a simple and effective way to do this.

In the pub this evening I came up with the following ...

def iterate(F, X, T, **params):
    """Iterate a non-linear map F starting from some initial condition X for T periods."""
    t = 0
    while t < T:
        yield X
        X = F(X, **params)
        t += 1

... test case using a Tinkerbell card ...

def tinker_bell_map(X, a, b, c, d):
    return [X[0]**2 - X[1]**2 + a * X[0] + b * X[1], 2 * X[0] * X[1] + c * X[0] + d * X[1]]

... gives ...

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 10, a=0.9, b=-0.6013, c=2.0, d=0.5)]
1 loops, best of 3: 26 µs per loop

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 100, a=0.9, b=-0.6013, c=2.0, d=0.5)]
1 loops, best of 3: 254 µs per loop

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 1000, a=0.9, b=-0.6013, c=2.0, d=0.5)]
1 loops, best of 3: 2.36 ms per loop

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 10000, a=0.9, b=-0.6013, c=2.0, d=0.5)]
1 loops, best of 3: 19.6 ms per loop

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 100000, a=0.9, b=-0.6013, c=2.0, d=0.5)]
1 loops, best of 3: 192 ms per loop

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 1000000, a=0.9, b=-0.6013, c=2.0, d=0.5)] 
1 loops, best of 3: 2.02 s per loop

%timeit -n 1 -r 3 [X for X in iterate(tinker_bell_map, [-0.72, -0.64], 10000000, a=0.9, b=-0.6013, c=2.0, d=0.5)]
1 loops, best of 3: 20.5 s per loop

... I tried several other test cases for deterministic and stochastic systems, and the above works are like charm. Although I think this is very good, I wonder if it can be done even faster using Numba ?

, ...

@njit
def tinker_bell_map(X, params):
    out = [X[0]**2 - X[1]**2 + params[0] * X[0] + params[1] * X[1],
           2 * X[0] * X[1] + params[2] * X[0] + params[3] * X[1]]
    return out

def simulator_factory(F):

    @njit
    def simulator(initial_condition, T, params):
        """Iterate a non-linear map starting from some X for T periods."""
        X = np.empty((initial_condition.shape[0], T + 1))
        X[:, 0] = initial_condition  # here is the offending line!
        for t in xrange(T):
            X[:, t+1] = F(X[:, t], params)
        return X

    return simulator


def iterator_factory(F):

    @njit
    def iterator(X, T, params):
        """Iterate a non-linear map starting from some X for T periods."""
        t = 0
        while t < T:
            yield X
            X = F(X, params)  # this is the offending line!
            t += 1

    return iterator

... , ...

In [8]: f(np.array([-0.72, -0.64]), 10, np.array([0.9, -0.6013, 2.0, 0.5]))
---------------------------------------------------------------------------
TypingError                               Traceback (most recent call last)
<ipython-input-8-d4c0195e7f4e> in <module>()
----> 1 f(np.array([-0.72, -0.64]), 10, np.array([0.9, -0.6013, 2.0, 0.5]))

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/dispatcher.pyc in _compile_for_args(self, *args, **kws)
    163         assert not kws
    164         sig = tuple([self.typeof_pyval(a) for a in args])
--> 165         return self.compile(sig)
    166 
    167     def inspect_llvm(self, signature=None):

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/dispatcher.pyc in compile(self, sig)
    301                                           self.py_func,
    302                                           args=args, return_type=return_type,
--> 303                                           flags=flags, locals=self.locals)
    304 
    305             # Check typing error if object mode is used

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in compile_extra(typingctx, targetctx, func, args, return_type, flags, locals, library)
    593     pipeline = Pipeline(typingctx, targetctx, library,
    594                         args, return_type, flags, locals)
--> 595     return pipeline.compile_extra(func)
    596 
    597 

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in compile_extra(self, func)
    316                 raise e
    317 
--> 318         return self.compile_bytecode(bc, func_attr=self.func_attr)
    319 
    320     def compile_bytecode(self, bc, lifted=(), lifted_from=None,

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in compile_bytecode(self, bc, lifted, lifted_from, func_attr)
    325         self.lifted_from = lifted_from
    326         self.func_attr = func_attr
--> 327         return self._compile_bytecode()
    328 
    329     def compile_internal(self, bc, func_attr=DEFAULT_FUNCTION_ATTRIBUTES):

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in _compile_bytecode(self)
    580 
    581         pm.finalize()
--> 582         return pm.run(self.status)
    583 
    584 

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in run(self, status)
    207                     # No more fallback pipelines?
    208                     if is_final_pipeline:
--> 209                         raise patched_exception
    210                     # Go to next fallback pipeline
    211                     else:

TypingError: Caused By:
Traceback (most recent call last):
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 201, in run
    res = stage()
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 415, in stage_nopython_frontend
    self.locals)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 710, in type_inference_stage
    infer.propagate()
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 408, in propagate
    self.constrains.propagate(self.context, self.typevars)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 113, in propagate
    loc=constrain.loc)
TypingError: Internal error at <numba.typeinfer.CallConstrain object at 0x10c5a7d50>:
Caused By:
Traceback (most recent call last):
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 201, in run
    res = stage()
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 415, in stage_nopython_frontend
    self.locals)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 709, in type_inference_stage
    infer.build_constrain()
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 395, in build_constrain
    self.constrain_statement(inst)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 519, in constrain_statement
    self.typeof_assign(inst)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 555, in typeof_assign
    self.typeof_expr(inst, inst.target, value)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 672, in typeof_expr
    raise NotImplementedError(type(expr), expr)
NotImplementedError: (<class 'numba.ir.Expr'>, build_list(items=[Var($0.27, /Users/drpugh/Research/python-dev/ramseyPy/sandbox.py (18)), Var($0.52, /Users/drpugh/Research/python-dev/ramseyPy/sandbox.py (19))]))

Failed at nopython (nopython frontend)
(<class 'numba.ir.Expr'>, build_list(items=[Var($0.27, /Users/drpugh/Research/python-dev/ramseyPy/sandbox.py (18)), Var($0.52, /Users/drpugh/Research/python-dev/ramseyPy/sandbox.py (19))]))
File "sandbox.py", line 45

Failed at nopython (nopython frontend)
Internal error at <numba.typeinfer.CallConstrain object at 0x10c5a7d50>:
Caused By:
Traceback (most recent call last):
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 201, in run
    res = stage()
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 415, in stage_nopython_frontend
    self.locals)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 709, in type_inference_stage
    infer.build_constrain()
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 395, in build_constrain
    self.constrain_statement(inst)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 519, in constrain_statement
    self.typeof_assign(inst)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 555, in typeof_assign
    self.typeof_expr(inst, inst.target, value)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 672, in typeof_expr
    raise NotImplementedError(type(expr), expr)
NotImplementedError: (<class 'numba.ir.Expr'>, build_list(items=[Var($0.27, /Users/drpugh/Research/python-dev/ramseyPy/sandbox.py (18)), Var($0.52, /Users/drpugh/Research/python-dev/ramseyPy/sandbox.py (19))]))

Failed at nopython (nopython frontend)
(<class 'numba.ir.Expr'>, build_list(items=[Var($0.27, /Users/drpugh/Research/python-dev/ramseyPy/sandbox.py (18)), Var($0.52, /Users/drpugh/Research/python-dev/ramseyPy/sandbox.py (19))]))
File "sandbox.py", line 45

... factory...

In [10]: s = simulator_factory(tinker_bell_map)

In [11]: s(np.array([-0.72, -0.64]), 10, np.array([0.9, -0.6013, 2.0, 0.5]))
---------------------------------------------------------------------------
TypingError                               Traceback (most recent call last)
<ipython-input-11-049d0797e27e> in <module>()
----> 1 s(np.array([-0.72, -0.64]), 10, np.array([0.9, -0.6013, 2.0, 0.5]))

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/dispatcher.pyc in _compile_for_args(self, *args, **kws)
    163         assert not kws
    164         sig = tuple([self.typeof_pyval(a) for a in args])
--> 165         return self.compile(sig)
    166 
    167     def inspect_llvm(self, signature=None):

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/dispatcher.pyc in compile(self, sig)
    301                                           self.py_func,
    302                                           args=args, return_type=return_type,
--> 303                                           flags=flags, locals=self.locals)
    304 
    305             # Check typing error if object mode is used

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in compile_extra(typingctx, targetctx, func, args, return_type, flags, locals, library)
    593     pipeline = Pipeline(typingctx, targetctx, library,
    594                         args, return_type, flags, locals)
--> 595     return pipeline.compile_extra(func)
    596 
    597 

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in compile_extra(self, func)
    316                 raise e
    317 
--> 318         return self.compile_bytecode(bc, func_attr=self.func_attr)
    319 
    320     def compile_bytecode(self, bc, lifted=(), lifted_from=None,

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in compile_bytecode(self, bc, lifted, lifted_from, func_attr)
    325         self.lifted_from = lifted_from
    326         self.func_attr = func_attr
--> 327         return self._compile_bytecode()
    328 
    329     def compile_internal(self, bc, func_attr=DEFAULT_FUNCTION_ATTRIBUTES):

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in _compile_bytecode(self)
    580 
    581         pm.finalize()
--> 582         return pm.run(self.status)
    583 
    584 

/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in run(self, status)
    207                     # No more fallback pipelines?
    208                     if is_final_pipeline:
--> 209                         raise patched_exception
    210                     # Go to next fallback pipeline
    211                     else:

TypingError: Caused By:
Traceback (most recent call last):
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 201, in run
    res = stage()
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 415, in stage_nopython_frontend
    self.locals)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 710, in type_inference_stage
    infer.propagate()
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 408, in propagate
    self.constrains.propagate(self.context, self.typevars)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 107, in propagate
    constrain(context, typevars)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/typeinfer.py", line 304, in __call__
    (ty, it, vt), loc=self.loc)
TypingError: Cannot resolve setitem: array(float64, 2d, C)[(slice3_type, int32)] = array(float64, 1d, C)
File "sandbox.py", line 29

Failed at nopython (nopython frontend)
Cannot resolve setitem: array(float64, 2d, C)[(slice3_type, int32)] = array(float64, 1d, C)
File "sandbox.py", line 29
```

.

Numba ...

+4
1

, -, njit, , nopython. , , njit. , njit -ed tinker_bell_map ( python), .

, , , , numba . (numpy 1.9.2 numba 0.14)

import numba
from numba import *
from numpy import *
import numpy as np

@njit
def simulator(initial_condition, params, X):
    a = params[0]
    b = params[1]
    c = params[2]
    d = params[3]
    X[0, 0] = initial_condition[0]
    X[1, 0] = initial_condition[1]
    for t in range(1, X.shape[1]):
        u = X[0, t-1]
        v = X[1, t-1]
        X[0, t] = u**2 - v**2 + a * u + b * v
        X[1, t] = 2 * u * v + c * u + d * v
    return X


x0 = np.array([-0.72, -0.64])
params = np.array([0.9, -0.6013, 2.0,0.5])

xs = np.zeros((2, 10000000 ))
%timeit -n 1 -r 3 simulator(x0, params, xs)
1 loops, best of 3: 70.7 ms per loop

xs = np.zeros((2, 100000000 ))
%timeit -n 1 -r 3 simulator(x0, params, xs)
1 loops, best of 3: 715 ms per loop


@njit
def tinker_bell_map(X, params, out):
    out[0] = X[0]**2 - X[1]**2 + params[0] * X[0] + params[1] * X[1]
    out[1] = 2 * X[0] * X[1] + params[2] * X[0] + params[3] * X[1]

def simulator_factory(f):
    def simulator(x0, params, x):
        for i in xrange(2):
            x[i,0] = x0[i]
        for t in xrange(1, x.shape[1]):
            f(x[:,t-1], params, x[:,t])
        return x
    return njit(simulator)

xs = np.zeros((2, 10))
sim = simulator_factory(tinker_bell_map)
print sim(x0, params, xs)

:

xs = np.zeros((2, 10000000 ))
%timeit -n 1 -r 3 sim(x0, params, xs)
1 loops, best of 3: 272 ms per loop

xs = np.zeros((2, 100000000 ))
%timeit -n 1 -r 3 sim(x0, params, xs)
1 loops, best of 3: 2.73 s per loop
+4

More articles:


All Articles