It is possible to have random seed input as an argument
(e.g. itoEuler(f, G, y0, tspan, random_state=1)) for reproducibility purpose?

Thanks

Hi guys, i'm citing you, and I believe you need to add a citation format to the .md file
I got this far...

  @software{sdeint,
    author = {Matthew Aburn and Yoav Ram},
    url = {https://github.com/mattja/sdeint},
  }

thanks a lot for implementing all this!
BTW, we noticed that the algorithm for stochastic Ito integration you have is order 1.0, however Rossler paper is titled "second order RK...", and mathematica has a routine that implements order 3/2 (in dt). Is this correct ? Maybe i'm misleading what second order means here (probabily \sqrt{dt}^2), but are you planning to implement higher order algorithms?
thanks again

sdeint needs some tighter input validation, if people try to use a state vector y with shape (n,1) instead of (n,)
e.g. mattja/nsim#2

I see that scipy.integrate.odeint also refuses this kind of input. But it gives a proper error message:

RuntimeError: The array return by func must be one-dimensional, but got ndim=2.

odeint has args=[] for this purpose. In case of sdeint we'd possibly need two separate parameters, argsf and argsG.

Hi,

thanks for the cool package. However, there seems to be a small issue with the return value of the input function

if I return lists in function f and function G, the integrator seems to be confused and then allocate a huge amount of memory. After manually converting my return values into numpy.array, the problem is solved. It took me an hour to figure this out.

Therefore, it would be easier if the integrator automatically checks the return values of the input function, to avoid confusing the reader.

Thanks!

Great package, works very well.
It could be useful to add some examples to the README.
Here's one I made for the logistic growth model:
https://gist.github.com/yoavram/5550c3f5daf1772e8352

Hi! Great library.

I wonder how can I get repeated simulations with differente noise realizations

I was tring this:

import numpy as np
import sdeint

M = 10

b = 1.0
tau = 0.8
tspan = np.linspace(0.0, 5.0, 5001)
x0 = 0.1*ones(M)

def f(x, t):
    return -x/tau

def g(x, t):
    return b*np.ones(x.shape)

result = sdeint.itoSRI2(f, [g]*M, x0, tspan)

But the outputs are all the same.

Thanks in advance.

I have a mechanical system described as d/dt (q) = f(q, t, F), where
q …generalized coordinates,
F…external forces to be modeled as white noise
(Where I used sympy.physics.mechanics to generate these equations of motion)

What I did in order to apply SDEINT, I linearized about the forces F:
B = f(q, t, F).diff(F).subs({F: 0.})
F1 = f(q, t, F).subs({F: 0})

Then I calculated as
Result = itoint(F1, B, x0, tspan)

Is this the right approach, or what would be the correct one?

Thanks for any help!

I hope it's not a too stupid question since I have no background in SDEs in general but I need to use the package to generate sample paths for our machine learning project.

From my understanding, suppose I have a set of SDEs with d states, if I assume independent additive constant diffusion noise, then the G function should return a dxd diagonal matrix right? Should the entries of the matrix contain the "variance or std" of the noise? I am not even sure if the terms variance and std are proper here.

Thanks a lot for helping.

Although scipy.integrate.odeint does not support that sort of thing.
But hey, why not?

This would then be an improved fix for mattja/nsim#2 , since it's reasonable for people to
supply a state of shape (n,1) interpreted as a column vector and expect it to work.

Implement with numpy.ravel on input and numpy.reshape on output but make sure data is not copied during reshape.

Nice project. I cannot see boundary conditions implemented (please correct me if I am wrong) which would make a easy and very powerful enhancement!

I tried your 1D example given in the introduction, and it seemed to work fine.
Then I tried the 2D example, and it did not seem to work:
No matter what I entered, the graphs of both solutions dropped from the starting point straight to zero and remained there.

where is my mistake?
Thanks for any help!

There's a bug that I suppose stems from the type interpretation of initial conditions. The following code which simulates a simple damped harmonic oscillator (with no noise at all) illustrates the error:

import numpy as np
import matplotlib.pyplot as plt
import sdeint


w0 = np.pi  # arbitrary
beta = 0.5  # arbitrary
tspan = np.linspace(0.0, 5.0, 5001)

x0 = 1, 1    # types interpreted as int (wrong solution) 
x1 = 1, 1.   # types interpreted as int (wrong solution)
x2 = 1., 1   # types interpreted as float (correct solution)

def f(x, t):
    return np.array([x[1], -w0*x[0]-beta*x[1]])

def g(x, t):
    return np.diag([0, 0])

result0 = sdeint.itoint(f, g, x0, tspan)
result1 = sdeint.itoint(f, g, x1, tspan)
result2 = sdeint.itoint(f, g, x2, tspan)

# vis the results
for i in range(3):
    l = plt.plot(tspan, eval('result'+str(i))[:,0], label=str(I), linewidth=5-2*i)
    plt.plot(tspan, eval('result'+str(i))[:,1], '--', color=l[0].get_color(), linewidth=5-2*i)
plt.legend()
plt.show()

Might be related to #15 and #16.

I'd like to use sdeint to generate a couple of figures. :)

How should I cite sdeint for that?

I'd like to brighten the possibility of integrating stochastic delay ODEs, if possible.

Can the approach used in, https://github.com/Zulko/ddeint/blob/master/ddeint/ddeint.py, be adapted for sdeint? I'm quite happy to write the code, if there's interest. Wanted to see if that looks like a viable approach. My understanding of numerical methods for stochastic ODEs is not good enough to be sure this is a sane route.

Hi,
I've been using this package for my master's thesis. I ran into problems with high-dimensional systems, where an unreasonable amount of memory was required. I rewrote the main loop for the itoint solver to generate the Brownian increments to be generated on the fly. This may be slightly lower but allowed me to handle much larger systems.
I would like to share this with you since this package helped me a lot.
I could either incorporate this into the current solver or write an additional option for larger systems. Please let me know your preferences.

I am using sdeint as it is shown in the 2nd example. So multi-dimensional independent driving Wiener processes. Is it possible to save the noise for every dimension/timestep?

I love to play around (that is all I am doing, I am semi retired) with SDEINT.
It gives possible realizations of the stochastic system. ( I usually use sympy.physics.mechanics to get equations of motion for some mechanical system, and then ‚exite it‘ with white noise forces).

It would be nice, if in addition, one could get the expectation and the variance of the solution.
Something like:

R, E, V = itoint(f, G, times, all=True) giving everything
R = itoint(f, G, times, all=False) giving only a realisation

Just a thought.

 pip install sdeint
Collecting sdeint
  Downloading sdeint-0.2.0.tar.gz
    Complete output from command python setup.py egg_info:
    Traceback (most recent call last):
      File "<string>", line 1, in <module>
      File "C:\Users\yoavram\AppData\Local\Temp\pip-build-ptlgr1og\sdeint\setup.py", line 48, in <module>
        long_description=read('README.rst'),
      File "C:\Users\yoavram\AppData\Local\Temp\pip-build-ptlgr1og\sdeint\setup.py", line 11, in read
        return codecs.open(os.path.join(here, *parts), 'r').read()
      File "D:\Anaconda3\lib\encodings\cp1255.py", line 23, in decode
        return codecs.charmap_decode(input,self.errors,decoding_table)[0]
    UnicodeDecodeError: 'charmap' codec can't decode byte 0x9f in position 1592: character maps to <undefined>

    ----------------------------------------
Command "python setup.py egg_info" failed with error code 1 in C:\Users\yoavram\AppData\Local\Temp\pip-build-ptlgr1og\sdeint

Same thing happens when installing directly from the code downloaded from github.

Python 3.4, Windows 7 64-bit.

Hello, big man. I want to ask a question. I have made a two-dimensional project according to your method. I always can't get the result I want. Can you explain a two-dimensional example?

First of all, thanks a lot for the package, it's very useful!

I was wondering whether the following application could be a relevant addition for sdeint:
In several fields, including neuroscience and decision making, it's often of interest to calculate the first passage time of (simple) stochastic processes reaching a certain bound or bounds (usually constant for simplicity although it doesn't need to be). You might have heard about the DDM model and variants.

Your functions such as itoSRI2() are very useful to calculate the trajectories of this process, but one still needs to evaluate when the trajectories reach the bound for the first time. To do this, most commonly one integrates trajectories until a fixed maximum (long) time Tmax, and then finds when the trajectories reach the bound. This is not very efficient, as many trajectories will reach the bound much earlier than Tmax. The other alternative is to evaluate the bound crossing at every time step of the trajectory, but since this is outside of the itoSRI2() function, it's also probably not efficient (multiple calls to the function with subsequent overhead, etc.)

I was wondering if it would be worthwhile to implement a method in which you pass already the bound function(s), and it returns directly the first passage time, and which bound was reached first (in case of multiple), all from within the itoSRI2() akin method. I think it would be valuable for a lot of researchers that currently are using the first technique described above, or implementing their own (and likely inefficient) Euler-Maruyama simulations.

Hello,

This is a follow up on #14 given numpy's latest recommendation regarding random number generation, which can be found here.

In short, using np.random.XXX is now legacy, and the recommended usage is:

rng = np.random.default_rng(seed=0)
rng.normal()

And instead of setting the global seed, one should pass the Generator instance to the place where it is needed.
AFAIK, sdeint does not allow passing a Generator instance.

It would be nice to incorporate this feature for compatibility with newer code, e.g., with sdeint using a global Generator instance that could be set by the user.

Thanks

Hi, thanks for your really great library! it really saves my day.

I'm not expert in maths, but am currently interested in simulating sde include a poisson process. is it possible or would you incorporate a poisson process dNt into your library in the future? it could be great if you can suggest some ways to implement it. Thanks for your help in advance!

I am playing around with sdeint.itoint a lot - and enjoy it!
It seems to me that when I integrate a given interval, say, (0, 10), the more 'steps' I give the more 'accurate' it gets - and the longer it calculates.
With scipy.odeint this does not seem to be the case.
Of course, only my observation, I do not really know how either program works internally.

Is my observation correct, or am I doing something wrong?

Thanks!