QuantEcon.py icon indicating copy to clipboard operation
QuantEcon.py copied to clipboard

Published 20 hours ago • verified publisher icon QuantEcon

Efficient iteration of arbitrary non-linear dynamical system (deterministic or stochastic)

Open davidrpugh opened this issue 9 years ago • 6 comments

For a while I have been wondering what the most efficient way to simulate an arbitrary non-linear (deterministic of stochastic) dynamical system in Python. I end up doing this a lot either teaching or research. I am convinced that there must be a simple and efficient way of doing this.

At 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

...a test case using the Tinkerbell Map...

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]]

...yields...

%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 have tried several other test cases for deterministic and stochastic systems and the above works like a charm. While I think the above is pretty good, I wonder if it could be made even faster using Numba? I tried doing

from numba import jit

@jit
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

But this didn't lead to any real speedups. Thoughts on how to speed this up in Numba? Thoughts on a better way to accomplish the above? Thoughts on whether something like this would be useful for Quantecon?

davidrpugh avatar Jul 23 '15 22:07 davidrpugh

Apparently numba support for generators in nopython mode in 0.18.1 (http://numba.pydata.org/numba-doc/0.18.1/release-notes.html) so in principle it should work. However, passing keyword arguments or unpacking them is not supported, which is probably the reason the code remains slow. Maybe if you pass parameters as a vector instead of a dict it would work. I am not sure though, I would need to go to the pub to get a clearer view ;-)

albop avatar Jul 23 '15 22:07 albop

@albop

The naive Numba interpretation...

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


@jit
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

...fails with the following...

In [24]: iterate([-0.72, -0.64], tinker_bell_map, 10, [0.9, -0.6013, 2.0, 0.5])
---------------------------------------------------------------------------
LoweringError                             Traceback (most recent call last)
<ipython-input-24-2e4b4320a906> in <module>()
----> 1 iterate([-0.72, -0.64], tinker_bell_map, 10, [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:

LoweringError: 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 515, in stage_objectmode_backend
    res = self._backend(lowerfn, objectmode=True)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 493, in _backend
    lowered = lowerfn()
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 471, in backend_object_mode
    self.flags)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/compiler.py", line 750, in py_lowering_stage
    lower.lower()
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/lowering.py", line 88, in lower
    self.genlower.lower_next_func(self)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/generators.py", line 128, in lower_next_func
    entry_block_tail = lower.lower_function_body()
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/lowering.py", line 140, in lower_function_body
    self.lower_block(block)
  File "/Users/drpugh/anaconda/lib/python2.7/site-packages/numba/lowering.py", line 158, in lower_block
    raise LoweringError(msg, inst.loc)
LoweringError: Internal error:
TypeError: decref() takes exactly 2 arguments (3 given)
File "<ipython-input-22-c3f19d5a7249>", line 6

Failed at object (object mode backend)
Internal error:
TypeError: decref() takes exactly 2 arguments (3 given)
File "<ipython-input-22-c3f19d5a7249>", line 6

...I feel like this is some obvious mistake that only a Numba newbie would make. Google failed to shed light on this stack trace...I am hoping that you have some idea...

davidrpugh avatar Jul 27 '15 14:07 davidrpugh

...hum, I honnestly have no idea what the error means. But there are two potential problems:

  • returning a list doesn't work in nopython mode. You can try with a tuple or write preallocated arrays.
  • I should have spotted that first but you are not allowed to pass functions as arguments, at least for now. I have had this problem many times (cf https://groups.google.com/a/continuum.io/forum/#!searchin/numba-users/function$20argument/numba-users/xiEQHOOe-O4/30klEMS6m4EJ) and have found some workaround, but it's not perfect so I don' t know if you will want to go down that route (example: https://github.com/albop/numba_experiments)

albop avatar Jul 27 '15 15:07 albop

@albop

Thanks for the links. I did notice earlier today that passing a function as an argument is problematic. Here are two tentative solutions that I have been playing with...

@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

...neither works unfortunately...

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

...and for the simulator 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

Here the issue seems to be with my trying to assign an array to a slice.

Finding Numba a bit frustrating at the moment...

davidrpugh avatar Jul 27 '15 20:07 davidrpugh

I have a scalar version of this working, but am having issues getting the slicing to work in array version. I thought I'd throw this out there in case it helps someone else tidy up what they're doing:

from numba import jit
import numpy as np
import timeit


def littlefunc(X):
    out = .95 * X
    return out

def simulator_factory_scalar(F):

    FJIT = jit(F, nopython=True)

    @jit(nopython=True)
    def simulator(initial_condition, T):
        """Iterate a non-linear map starting from some X for T periods."""
        X = np.empty(T + 1)
        X[0] = initial_condition
        for t in range(T):
            X[t+1] = FJIT(X[t])
        return X

    return simulator

test_scalar = simulator_factory_scalar(littlefunc)
test_scalar(50., 1000)

cc7768 avatar Jul 28 '15 16:07 cc7768

Partially solved with help from SO...

@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):

    @njit
    def simulator(X0, params, out):
        """Iterate a map F starting from some X0 for T periods."""
        for i in xrange(out.shape[0]):
            out[i, 0] = X0[i]
        for t in xrange(1, out.shape[1]):
            F(out[:, t-1], params, out[:, t])

    return simulator

This works and yields a factor of 60x speed up!

In [78]: xs = np.empty((2, 10000000))

In [79]: s = simulator_factory(tinker_bell_map)

In [80]: %timeit -n 1 -r 3 s(np.array([-0.72, -0.64]), np.array([0.9, -0.6013, 2.0, 0.5]), xs)
1 loops, best of 3: 372 ms per loop

davidrpugh avatar Aug 17 '15 14:08 davidrpugh