Wrapper class 'integrate.complex_ode' does not accept function parameters (Trac #1451)

[scipy-gitbot]

Original ticket http://projects.scipy.org/scipy/ticket/1451 on 2011-06-01 by trac user gideon, assigned to unknown.

Hi,
I have tried to use 'integrate.complex_ode' in conjunction with integrator.set_f_params(*args). I consistently get the message:
"'float' object is not subscriptable" (or 'int' or 'complex' depending on the parameters passed). This happens whether *args is one parameter or a list. Error seems to be in :

/usr/lib/python2.7/dist-packages/scipy/integrate/ode.py in: 

_wrap(self, t, y, *f_args)

    370     def _wrap(self, t, y, *f_args):
--> 371         f = self.cf(*((t, y[::2] + 1j*y[1::2]) + f_args))
    372         self.tmp[::2] = real(f)
    373         self.tmp[1::2] = imag(f)

The wrapper complex_ode seems to work fine without trying to pass any additional function parameters.

Code to show problem:

from scipy import integrate
from pylab import *

def f(t, y, arg):
    f = -1j*arg*exp(1j*pi)*y
    return f


ic = [1.0]
t0 = -.015
t1 = .015
dt = 0.1

solver = integrate.complex_ode(f)
solver.set_integrator('vode',method='bdf',rtol=1.e-12)
solver.set_f_params(arg)

t = []
sol = []
solver.set_initial_value(ic, t0)
while solver.successful() and solver.t < t1:
    solver.integrate(t1,step=1)
    t.append(solver.t)
    sol.append(solver.y)
print solver.y
t = array(t)
sol = array(sol)
plot(t,abs(sol))

I am running:[[BR]]

• Linux 2.6.38-8-generic Allow loading from file-like in scipy.io.arff #42-Ubuntu SMP Mon Apr 11 03:31:50 UTC 2011 i686 i686 i386 GNU/Linux
• Python 2.7.1+ (r271:86832, Apr 11 2011, 18:05:24)

*IPython 0.10.1

Thanks[[BR]]

-Gideon

[scipy-gitbot]

trac user borishim wrote on 2012-03-20

After a quick reading of the source code, I think that to fix this, a major lifting in the f2py interface is required, and all integrators needs to be fixed separately. For example, for the dopri5 integrator, additional arguments will not work even for real valued problems. For the workaround, functools.partial or lambda can be used in place of the generator, although performance hits might be experienced:

import functools
f_with_args = functools.partial(f, args=1)
...
solver = integrate.complex_ode(f_with_args)
...

I actually am starting to think that the entire rewritting of set_f_params() method in terms of functools.partial could be a better design decision, considering that the target function is a python function after all.

[WarrenWeckesser]

Another report, on stackoverflow: http://stackoverflow.com/questions/34577870/using-scipy-integrate-complex-ode-instead-of-scipy-integrate-ode

[pathfinder49]

It seems like this issue is still not fixed! Why is the wrapper implemented in such a complicated way? All that should be needed here is numpy.ndarray.view(dtype=np.floatxx), where np.floatxx should be chosen based on the complex dtype. This results in an interleaved vector of floats and is suitable for direct use with scipy.integrate.ode.

This approach also avoids unnessecary copying of input arrays.

Error messeges could also be more helpful. In two different use-cases I got:

• TypeError: unhashable type: 'slice'
• TypeError: can't multiply sequence by non-int of type 'complex'

Edit: Minimally complex_ode should be clearly documented as BROKEN.

[pathfinder49]

This behaviour seems like it could be integrated into scipy.integrate.ode. However, I appreciate that changing the behaviour of ode may not be worth it.

If someone more experienced with scipy development could give their thoughts, I'm happy to make a PR. Depending on feedback this would be:

• Refactoring complex_ode to be a thin wrapper using numpy.ndarray.view or
• Adding complex support for all solvers in ode

[weinbe58]

I am pretty sure the ode and complex_ode interfaces are deprecated. Because of this, it is not clear if any of the Scipy developers are going to be willing to do this update.

Since ode accepts arguments, it is simple enough to write a wrapper function that uses the view of an array to get around this:

def real_func(t,y,complex_func,args):
    return complex_func(t,y.view(np.complex128),*args).view(np.float64)

I forget the order of the arguments, but this hopefully is clear. Also, do not forget to switch the access to the solution of the solver from solver.y to solver.y.view(np.complex128). It is a little annoying, but it is a pretty minor change to any code.

[pv]

Looking at this, it actually seems to be a f2py bug. It does not appear to correctly handle the `f(t, y, *args)` variable-argument construction in callbacks with extra arguments, as you can see also with `ode`: ``` from scipy import integrate def f(t, y, *args): return 0*y solver = integrate.ode(f) solver.set_integrator('vode') solver.set_f_params(1) # <- works, if you comment this out solver.set_initial_value([0, 0], 0) while solver.successful() and solver.t < 1: solver.integrate(1, step=1) ``` It's possible to work around this in complex_ode by hooking `complex_ode.set_f_params` to set an attribute on `self`, and not use `*args` in the wrapper routines. The f2py bug itself is a bit nastier to fix. Adding the workaround would be fine. Regarding using `np.asarray(x, np.complex128).view(np.float64)` for unraveling, and more extensive refactoring --- the former is simple change that can be made. The latter probably changes the API, and can't be made.

[pathfinder49]

Since ode accepts arguments, it is simple enough to write a wrapper function that uses the view of an array to get around this

Agreed complex_ode doesn't really do much and isn't really essential. However, unfamiliar users may well decide to use complex_ode as it seems like the intended way of integrating complex ODEs. (They may also be unaware of view.) Unfortunatly, the errors from this issue don't point the user to complex_ode being broken. I also wouldn't expect inexperienced users to find this issue.

I am pretty sure the ode and complex_ode interfaces are deprecated.

They are? Where can I see this?

Regarding using np.asarray(x, np.complex128).view(np.float64) for unraveling, and more extensive refactoring --- the former is simple change that can be made.

I'll submit a PR for np.asarray(x, np.complex128).view(np.float64) when I have a moment.