real and complex noises?

[AshtonSBradley]

I am working on setting up some quantum phase space simulations. The issue is that I would like to be able to pass f and g to SDEProblem. However, when I pass f and go to SDEProblem, although complex variables are supported, the noise is called with the integrated type. When prob is passed to solve, it triggers the error

MethodError: no method matching randn(::MersenneTwister, ::Type{Complex{Float64}})
Closest candidates are:
randn{T}(::AbstractRNG, ::Type{T}, ::Tuple{Vararg{Int64,N}}) at random.jl:1216
randn{T}(::AbstractRNG, ::Type{T}, ::Integer, ::Integer...) at random.jl:12

It seems that this is triggered by array problems: if I try to solve the single variable Kubo problem (which requires a real noise), everything goes well, problem solved.

So there are two issues

1. can noises be more directly called and manipulated to construct a new noise in general? I would like to be able to, say, take real Gaussian noises, construct an object in momentum space, then Fourier transform (or matrix transform depending on the basis) to position space, and use the new spatially coloured noise in a position space SDE. This would allow for a broader class of SDEs as, for example, the spatial basis would presumably be free for the user to construct.

2. is there future scope to have the option of passing a complex variable problem to solve, but having access to real or complex noises, rather than triggering noise from the integrated type? Access to the real noises would be great, allowing construction of complex noises if desired:

from the standard Wiener process dw1, dw2, a unit-normalised complex noise is constructed as

dW = (dw1 + i*dw2)/sqrt(2)

It would be great to be able to define a vector drift and factorized diffusion matrix for complex variables using @ode_def_bare, and pass this to SDEProblem, with a some choice as to whether complex or real noises are assumed in the product B*dW, or whether a more direct noise construction is needed.

[ChrisRackauckas]

This error is because randn cannot create complex numbers. So, to answer point 1, yes, but it's not flexible enough yet. You can see here:

https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/problems/sde_problems.jl#L19

That you can pass in a noise process, but

https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/noise_process.jl#L9

it's just a function. It should be a function (t,u), and for inplace also take in the rand_vec. That will address part 1.

It will even go a little further. One way is that it will allow many forms of temporal noise. One way these are implemented is by generating them from an SDE. This will allow one to use a system of equations, and make say the first equation have white-noise via randn(), and then use the first equation's solution as the noise process for another.

I am thinking of actually expanding it all the way to (rand_vec,t,u,W), where W is the history vector of the noise. It's not always saved, but when it isn't it's just a [] so it'll never error. This will allow someone to define non-Markovian noise which relies of previous values of the noise. But the standard noise process would still be (rand_vec,t,u,W) -> randn!(rand_vec) which just does an inplace update to the random vector with white noise, and ignores all else.

For point 2, I will make the default always type matching. However, one the noise type is more flexible like this, a user who wants real noise with a complex-valued independent variable will be able to easily define a new NoiseProcess which does that. I like the option of "simple defaults, and unlimited tweakability".

dW = (dw1 + i*dw2)/sqrt(2)

Thanks, I'll use this.

It would be great to be able to define a vector drift and factorized diffusion matrix for complex variables using @ode_def_bare, and pass this to SDEProblem, with a some choice as to whether complex or real noises are assumed in the product B*dW, or whether a more direct noise construction is needed.

I think it would be easiest to just allow the direct construction. It's not difficult to write the function, but it would heavily clutter the API to have a lot of "complex-variable specific" parts around. Again, I think it should be supported for users to tweak these parts, but the "simple path" needs to define where it ends, otherwise it will end up not so simple!

[AshtonSBradley]

Yes, I think the simplest situation is: trigger the noise based on integrated type, and allow direct construction for those cases not covered by the default. This means that

1. if I pass a complex variable, I get a default complex noise constructed according to
   dW(t)=[dw1(t)+i*dw2(t)] /sqrt(2)

with with dwj standard Gaussian Wiener processes having
mean(dwj)=0 and
mean(dwj dwk)=dt*\delta_{jk}. For the complex noises as constructed above

mean(dW^*(t)dW(t))=dt

and mean(dW(t)^2) = mean([dW(t)^*]^2)=0.

2. If I pass a real variable (or array of them) for integration, I would have the current behaviour where real noises dwj are returned as elements in the vector dW(t), one for real each variable.

In this context the case of full generality should then, for each complex variable passed, allow access to two independent real noises, one for each real valued degree of freedom. Presumably you could allow an arbitrary number to be called, and this would cover the common case where less noises are needed than variables.

The main reason for this: complex variables allow much simpler "pseudo-code" input. At the moment I would need to reduce my problem to a much more complicated set of equations for each of the real degrees of freedom. For a few variables, no problem, but for multimode fields, this quickly becomes very unwieldy. What you are proposing could be used to construct stochastic field theory problems quite directly, I believe. The real advantage of doing this in Julia: complete freedom of setup of the integrated type. One should not be restricted to certain discretizations of the fields, and can build arbitrary setup and analysis code around the integrator. Super cool!

[ChrisRackauckas]

In this context the case of full generality should then, for each complex variable passed, allow access to two independent real noises, one for each real valued degree of freedom. Presumably you could allow an arbitrary number to be called, and this would cover the common case where less noises are needed than variables.

Expand on this a little bit more?

You mean like allow multiple dW processes? That can be done by making the problem take in a tuple, and have g a tuple which matches its size. If the functions inside of those tuples are evaluated in a recursive manner, the recursion will expand and won't have type-stability or performance issues. It will add a little bit to the code, but shouldn't be bad.

[AshtonSBradley]

So if I understand: I could pass a tuple of coefficients that forms a linear transformation acting on two real noises? Then for each complex variable z[i] = x[i] + i*y[i], I would have (summing over repeated index j)

dz[i] = A(z,t)[i]dt + B(z,t)[i,j]dW(t)[j]

and default would be dW[j] is complex Wiener as above (two real noises for each complex variable), but in general I could pass a tuple of coefficients (a,b) allowing arbitrary linear combinations a*dw1[j] + b*dw2[j] ? I think this could be a nice way to go, making sure redundant calls to randn are suppressed for a or b equal to zero.

The general underlying problem is of the form
dX = A(X,t)dt + B*dW with all of X being real variables, and dW being a vector or real noises. It is just very painful to translate many quantum problems into real variables, and leads to code that is hard to modify and bug test. The complex variable statement is highly convenient, provided sufficient control is possible with the form of B. I think by allowing a tuple, this would expand the ease of use a lot, effectively giving easy access to complex or real noise with the same syntax. It should be equivalent to giving access to the full real valued dW, of length 2N for N complex input variables.

[ChrisRackauckas]

Yes. The proposal is to allow you to specify (g1,g2,g3,...) and then (NoiseProcess1,NoiseProcess2,...) to make

du = f(t,u)dt + g1(t,u)*dW1 + g2(t,u)*dW2 + ...

which is something that was recently lost with the inplace noise process update (though admittedly, the previous way was archaic and limited). Again, each default noise process would default to type-matching, so if u is a vector of Complex128, then by default that's what the dWs would make. Is that the same thing?

I realized that we will need a toggle in SDEProblem to set the noise type. I think it will just be noise_type with a default of eltype(u0). So to change to real-valued noise when u0 is an array of complex, it would be noise_type = Float64 as an additional kwarg. This cannot be done directly from the NoiseProcess because it needs u0 and it needs to be set in order to preallocate the arrays for the RNG.

[AshtonSBradley]

(sorry for some repetition here) I am not so sure that this level of generality is necessary. The general multimode stochastic problem can be reduced to the form

dX = A(X,t)dt + B(X,t)*dW(t)

where X is a vector of variables to be integrated, A is a drift vector of same dimension, and B is a factorisation of the diffusion matrix, i.e. D=B'*B, and dW is a vector of real noises with same dimension as X. We can talk about arrays of variables or arrays of complex variables as convenient inputs, but we can always think of them as a vector to make a clear connection with the noises. So, at least for Markov problems, we should never need more than N noises for a vector X of length N. So, it makes more sense to me, to allow problems to be constructed as

du[i]=A(u,t)[i]dt + B(u,t)[i,j]dW(t)[j]

which seems a bit more specific that what you have written above, that would generate a column vector of noises for each gj(t,u), is that right?

So I think if you replace the product g1(t,u)*dW1+g2(t,u)*dW2+... with the (typically very sparse) matrix product B*dW it would still span the same space of physical multimode diffusion, but with much fewer calls to the noise generator. Triggering noise to be complex based in input type would then be a convenient default for common semi-classical problems, and then access to all real noises would give full control over construction for fully quantum phase space problems.

Maybe I am missing something: is g1 allowed to be a matrix and dW1 a column vector in your proposal?

If that was the case, there could be a nice reason to have the more general constructor you propose by allowing more noise vectors if it means more scope for colored noise, fractional Brownian Motion, etc.

To answer your question:

Again, each default noise process would default to type-matching, so if u is a vector of Complex128, then by default that's what the dWs would make. Is that the same thing?

it would be the same, provided gj(t,u) can be a matrix operation. It would have full flexibility provided one also has the option of explicit construction, i.e. the user is allowed to make some controlled call to randn to create a vector of noises, then carry out arbitrary operations to create some new noise (in general colored).

[ChrisRackauckas]

I see what you're saying. Instead of having dW1 , dW2, ..., you can just have dW a matrix where columns would be those separate noise processes. Indeed, one advantage of this is that then the diffusion function g(t,u) could also be a matrix, and it would be easier to do matrix multiplications like this than it would be to manually decompose that multiplication into a bunch of separate functions. This would likely be much easier to implement inside of the solvers too.

The slightly difficult thing is that we allow more than just equations on vectors. The reason is that using matrices or tensors for u can be beneficial for some representations of SPDEs. However, it seems safe to define:

du[x1,x2,...,xN]=A(u,t)[x1,x2,...,xN]dt + B(u,t)[x1,x2,...,xN,x(N+1)]dW(t)[x(N+1)]

Let's clarify:

Maybe I am missing something: is g1 allowed to be a matrix and dW1 a column vector in your proposal?

In the current setup the g function is only allowed to return a size which matches u, and that matches the initial condition. So if u is a vector (for a vector system of equations), g had to return a vector. It seems we should allow this to return a matrix, with m columns where m is the Ito-dimension. Then, if u is an n-dimensional tensor, g should return and n+1 dimensional tensor where the first n dimensions match u in size and the last dimension is size m, the Ito-dimension.

A necessary input into the SDEProblem would then be the Ito-dimension. It would default to one. The type matching default would just make white noise of that size.

Example:

u = Vector{Float64}(10)
SDEProblem(f,g,u0,tspan)

Makes du[i]=f(u,t)[i]dt + g(u,t)[i]dW(t), but

SDEProblem(f,g,u0,tspan,ito_dims = 5)

Makes du[i]=f(t,u)[i]dt + g(t,u)[i,j]dW(t)[j] where the j=1:5, and dW is a matrix of white noise. Here it would assume that g(t,u,du) is a function that writes inplace into du and size(du) = (10,5), while f(t,u,du) writes into du with size(du) = (10,). Then we also have

u = Matrix{Float64}(10,8)
SDEProblem(f,g,u0,tspan,ito_dims = 5)

Makes du[i,j]=f(t,u)[i,j]dt + g(t,u)[i,j,k]dW(t)[k] where the k=1:5. Here it would assume that g(t,u,du) is a function that writes inplace into du and size(du) = (10,8,5), while f(t,u,du) writes into du with size(du) = (10,8), giving a (10,8) solution at each step.

u = Matrix{Complex128}(10,8)
SDEProblem(f,g,u0,tspan,ito_dims = 5)

would be the same sizing as before, but now dW=[eta(t)+i*eta2(t)] /sqrt(2) for two indepedent BMs eta and eta2, i.e. dW is a complex-valued BM. If that's not what you want, you can do:

u = Matrix{Complex128}(10,8)
SDEProblem(f,g,u0,tspan,ito_dims = 5,noise=my_noise_process)

where we change the noise process by defining my_noise_process. The core of a NoiseProcess will be a function noise_func(t,u,W,rand_vec) which writes inplace into rand_vec (so rand_vec is the size of the dW's, so in this example (10,8,5). rand_vec would be typed to match u, but that shouldn't give problems. For example, to have this equation only have a real-valued BM, we would do:

noise_func(t,u,W,rand_vec) = randn!(rand_vec,Float64)

which will fill it with only real-valued numbers, which would convert to eta + 0*i. It's fine to do this conversion since the promotions in arithmetic would have to do it anyways, so this is likely faster at the cost of some memory. But this will also allow you to generate spatially colored noise, use a separate SDE to generate the noise for another component, generate temporally-colored noise all at once or iteratively using previous values of the process W, etc.

I think this is what we're looking for?

[ChrisRackauckas]

Some caveats to address:

Are element-wise operations on a Vector{Float64}(N) the same speed as Matrix{Float64}(N,1)? (For operations using the linear index)

What do we do about AbstractArray types which cannot extend up a dimension? Good example is https://github.com/JuliaDiffEq/MultiScaleArrays.jl . It seems natural in this case for the return type of g to be a Vector of the multi-scale array, since that will trivially generalize up. However, we wouldn't want to always do that, because instead of a Vector{Vector{Float64}} we will want a matrix, right? Maybe we need finer control over this type, or a specialization on Vector.

If these are worked out, then the above has a good implementation.

[AshtonSBradley]

I think this is likely a good way to go. Most problems would be handled by ito_dims=1, but some very general noise constructions would also be available. I think my only nagging (maybe unfounded) worry regards in-place versus matrix operations.

What if I have

u = Matrix{Complex128}(10,8)

but now instead of

du[i,j]=f(t,u)[i,j]dt + g(t,u)[i,j,k]dW(t)[k]

I want g to be a function of the noises:

du[i,j]=f(t,u)[i,j]dt + g(t,u,dW)[i,j]

For a common problem I look at it, the construction is:

1. sample noises in momentum space (actually in a basis of simple harmonic oscillator modes, in momentum space representation)

2. perform some local (but could be more general) operations in momentum space

3. inverse fourier transform to create a position space coloured noise (non delta-correlated)

4. do some operations in position space (local, but could be more general)

5. project the result into a set of complex coefficients representing the amplitudes of each oscillator mode comprising the stochastic wave function. This is du[i,j].

So getting g to write in place is trivial, but I want to be able to have noises as general inputs to the function g, rather than be forced to write it as a matrix operation (although the latter does cover many cases).

This may sound a bit perverse, but it is the general form of the optimised algorithm. There are other stochastic field theory problems that can be cast in this form. So from a solver point of view, f is just a function of u,t, while g=g(t,u,dW) needs to be able to perform arbitrary operations on the noise.

The problem comes down to having the freedom to a work in a convenient spatial basis, rather than being forced onto a regular spatial grid. For open quantum systems problems this becomes essential...

[ChrisRackauckas]

du[i,j]=f(t,u)[i,j]dt + g(t,u,dW)[i,j]

No, this is a random differential equation (RDE), not an SDE. There is no way to guarantee that the user will treat it as an SDE.

while g=g(t,u,dW) needs to be able to perform arbitrary operations on the noise.

No, it needs to not be able to do that. There are many methods for which more control over the application of the noise process to the function g is required. Milstein methods, stochastic RK, etc. Those cannot be written down if dW cannot be separated. For example, I don't think you can do any of these:

http://epubs.siam.org/doi/pdf/10.1137/09076636X

Anything involving iterated stochastic integrals wouldn't be possible. Only Euler methods will work with this (or of course, RDE methods).

For a common problem I look at it, the construction is:

Are you sure this construction isn't solving a related RDE? Do you have a paper showing the convergence of this? I am curious to see this all written out to see if I can support it properly.

Note that I am interested in the RDE case, and will get something going when Kloden releases his next book:

https://www.researchgate.net/publication/314242937_Random_Ordinary_Differential_Equations_and_their_Numerical_Solution

[AshtonSBradley]

Ok, it may be true that formally it is an RDE once you go to an arbitrary basis, but the underlying object is definitely an SDE. Papers are

http://journals.aps.org/pre/abstract/10.1103/PhysRevE.89.013302

(full numerical method, fairly lengthy)

and for shorter overview

http://journals.aps.org/pra/abstract/10.1103/PhysRevA.86.053634

Even just Euler would be fine though, if that is the only algorithm that is formally valid. Is it possible to let the user do it, and only allow Euler as an integrator? In fact, this is what we have always used as we know it is robust. Semi-implicit Euler for SDE in Stratonovich form.

We have used RK methods quite a bit for weak additive noise in this context, but the full Stochastic GPE needs a lower order method due to multiplicative noise.

I guess the point is this: it is only obviously an SDE in a particular representation, due to the fact that is coming from a projected field theory.

[ChrisRackauckas]

I see. Thanks for the reference. Yeah, even though it's an SDE in some representation, the way it's written is not an SDE unless you can make that dW not inside the projection operator. Because of that, these integration methods don't apply in this representation.

That said, it does follow the general form form for an RDE, so those methods definitely apply. And because the Euler methods still apply to RDEs, IIRC Euler-Maruyama -> Ito and Euler-Heun -> Stratonovich in this form. But there are methods which will do well for these. And adaptivity should carry over using the RODE-SDE pair approach, though a proof of convergence would technically be needed. Reference:

http://publikationen.ub.uni-frankfurt.de/frontdoor/index/index/docId/40146

I think Ito/Stratonovich-driven RODEs should get their own type, since they likely will be important enough. The can plug into the StochsticDiffEq.jl stuff though to get the adaptivity, interpolations, and event handling with a few modifications. But I am not going to commit to tackling this right now: this'll have to wait a bit.

[AshtonSBradley]

Ok, very interesting! I look forward to seeing how this develops. Putting aside the RDE issue (which would be awesome to solve), I think your proposal will expand the scope in very worthwhile ways.

It would be great if Euler was enabled for RODEs in the meantime, and you allowed that level of construction with noise. That would mean we could start writing code to your interface, and then hope to use more interesting algorithms in future as they get implemented. Is there scope for doing that?

[ChrisRackauckas]

Yeah. That's not difficult. I think it would be best to build solvers for Wiener process driven RODEs into StochasticDiffEq.jl. The roadmap would be:

• Setup the RODE types in DiffEqBase.jl
• Open up the StochasticDiffEq.jl dispatch so that way the RODE will go there
• Make a new algorithm for REM (Random Euler-Maruyama). Should do REH (Random Euler-Heun, for Strat problems. Should make EH as well while at it).
• Make a separate dispatch in StochasticDiffEq.jl for the cache construction
• Make a separate dispatch for perform_step on this new cache

That should be all that's needed? Event handling may be broken with this though, but most of the stuff should work.

I need to finish up a paper right now though.

[AshtonSBradley]

+1 for this, and for your paper getting finished

[ChrisRackauckas]

I made a bit of progress here. The complex numbers stuff is implemented, and so is the RODEs, but not the non-diagonal noise things. I have a lot of plans there, and will likely do it today. Note that none of this is "currently available": it's part of a 7 repository API change so please just wait for the release to try it.

Side Note: RODEs and The Projective Method

But I want to point out that what you have isn't easily seen as an RODE either. Those fall into the form:

u' = f(t,u,W)

that is, they use the solution of some random process. Also, in your project formulas, you need the dt too, for example

dψ|M = ≡ P((-i/hbar)*V(r,t)*ψ*dt + i*ψ*dW*M(r,t))

there, using ψ=u, it seems like you need

u' = f(t,u,dt,dW)

Again, I looked that up and am coming up blank. But, there is a convoluted way of writing that as an RODE... and I can use Lipschitz arguments to show certain forms of convergence. So the Euler type method

u_{n+1} = u_n + f(t,u,dt,dW)

where the directly uses the increments in-place actually makes sense from an algorithm perspective (it should converge), but this is a bad idea. No only is there no way to ensure that the "user does this correctly", but also you can't put any algorithms on this. For example, when I say it converges, the simple way of putting in the values for dt and dW into your equation converges in the the Ito sense. You'd have to, in your f, implement a Heun-like update to get Stratanovich convergence. So this doesn't sound good..

True Answer

In the algorithm you want to do, you just need to project after doing the update. Because we know that this converges (under appropriate assumptions), I think it's fine to explicitly do that projection to the SDE. I think the best way to implement this algorithm as an SDE would be in two steps. The first step is to define f and g as though there is no projection. Then use a DiscreteCallback. DiscreteCallbacks can be set to happen after every update using

condition = (t,u,integrator) -> true

Then do your affect:

function affect!(integrator)
  # do something to change integrator.u
  # i.e. do your projection here,
end

Then make this via cb = DiscreteCallback(condition,affect!).

Then you would use

sol = solve(prob,EM(),callback=cb)

That will make it use the Euler-Maruyama method, and project after each update (before each save, see the save_positions arguments on the callback for controlling the saving behavior). Note this converges in the Ito sense. For Stratanovich, you need

sol = solve(prob,EulerHeun(),callback=cb)

and I'll implement the EulerHeun() method shortly.

[AshtonSBradley]

First thoughts: this looks great! There are a bunch of advantages to casting the whole algorithm+formulation in this way. Projected functional calculus tries to achieve this same goal: set up a calculus that is basically as clean and useable as functional calculus (i.e. work as if there are no projectors), but with projectors at appropriate points to achieve a consistent representation of the field theory on a finite basis set.

Nagging worry: this is a shift of perspective so that $\psi(x)$ is the fundamental information in position space, instead of a set of coefficients $c_n$ giving the contribution of each of set of orthonormal eigenfunctions. Down side is I couldn't easily implement terms of different numerical order in a quadrature sense. The algorithm evaluates the projected nonlinear term in the Gross-Pitaevskii equation by transforming $c_n$ to a particular spatial quadrature grid (Gauss-Hermite typically) that will allow efficient computation of the projection of the nonlinear term $|\psi(x)|^2\psi(x)$ onto each basis function.

Let's say I want to implement a three-body loss term, where my nonlinear term in the RHS that I want to project is now of the form $|\psi|^4\psi$. I need to transform to position space on a different Gauss-Hermite grid, with different weights. So my question is, given your above definition, am I still allowed to pass $c_n$ as the fundamental information being propagated? It looks like the answer would be negative ...

Should I go back to position again to get back my field $u=\psi(x)$, that is the quantity being numerically integrated? Which grid should I go onto? The most common one? The grid is not fundamental, only the coefficients.

Put another way, it looks like you are saying "just project at the end of the step" which is cool, but:
Once I project, I am in a different basis representation. I am not really integrating $\psi$ anymore, but rather $c_n$. Does this matter? Perhaps I don't understand the true nature of the callbacks you describe.

Some kind of answer
Thinking through the implications, I think this would mean I would work on a specific quadrature grid that is of sufficient order to integrate all terms in my RHS. That can be done easily, but there will be a memory penalty. If I consider the GPE nonlinearity in 1D, then for N modes, I need a quadrature rule of order 2N to evaluated the RHS exactly, including projection. So it would mean carrying around the position space field on this larger grid. However, since the arrays have to be preallocated for in place operations to get performant code, the only issue would seem to be accepting some increased memory usage for the output fields from the solver (unless I can specify a special kind of output?). For the above example, in 3D, this would amount to saving fields that are 2 x larger in each dimension, so a factor of 8 extra memory compared with coefficient form.

[ChrisRackauckas]

I'm not sure I understand everything you're saying, but:

1. You can save before or after a callback. So you can either choose to save the c_n or psi.
2. You can use a closure or the parameters interface to pass the c_n along. They can also change over time like this.

Perhaps I don't understand the true nature of the callbacks you describe.

Callbacks are for injecting arbitrary code to act on the integrator type. The integrator type holds the full differential equation, so this gives you a way to do all of this special stuff.

[AshtonSBradley]

OK! I think this could be a totally kick ass solution.

One further thought: do the callbacks give enough freedom that the noise triggering can be set by them? If I take the simplest case of additive noise (simple-growth Stochastic Projected Gross-Pitaevskiie equation), the noise is additive in the coefficient space of the representation $\psi(x)=\sum_n c_n \phi_n(x)$, i.e. in setting up the RHS of the SDE, one does something like this:

dc_n = c_n +dW_n

Now if I need to pass the spatial field $u=\psi(x)$, that seems to mean that the noise would called based on that spatial field?

For additive noise this could be just fine (one would need to work out the coefficients of the diffusion matrix for the given quadrature grid). For the multiplicative noise terms, I would need to dig a bit deeper to see if it would allow noise construction.

Put another way - do I have enough access to the individual noises to create spatially coloured noise? So I want to be able to make the following (a true SDE in projected functional calculus):

d\psi(x) = \psi(x)B(x)dW(x)

on a grid x of my choosing (would be a nonuniform Gauss quadrature grid), and have the noise have some interesting spatial correlations:

<dW(x)dW(x')>=M(x-x')dt

For my specific case there is an easy way to construct this noise in momentum space, if I have access to Gaussian random noises triggered by the coefficients c_n. However, this is a fairly generic feature these kinds of problems: I need access to the noises (either real, or complex) per mode triggered by the coefficients, i.e. of the same dimension as c_n.

In practice this means that writing the SDE diffusion term as B*dW for some noise vector and diffusion matrix factorisation B is not so obvious, because it is formulated in a (projected) functional sense, i.e. what I mean by the vector dW(x), is now defined in the functional space in the sense that the vector to be solved is

u --> [d\psi(x); d\psi^*(x)] and the noises would be

dW --> [dW_1(x); dW_2(x)]

[ChrisRackauckas]

Put another way - do I have enough access to the individual noises to create spatially coloured noise?

You could... but you shouldn't do it this way. It will screw with the saving and all sorts of other things. Just create a NoiseProcess.

I implemented everything discussed here except the expanded API for a NoiseProcess. Docs will go up soon, but tags might take awhile. But what's there for NoiseProcess should already be able to handle spatially colored noise.

After this huge tagging set goes through, I'll do

1. Expanded NoiseProcess API (mostly for temporally-colored noise)
2. Help you get something together

Note that for Stratonovich, we only have an Euler-Heun Strong Order 0.5 method. I tried to put a derivative-free Milstein method in the same release, but had some issues:

SciML/StochasticDiffEq.jl#19

I didn't want that blocking the release, so I ignored it for now. But that should be pretty accessible to anyone who wants to help out since it's at the debugging stage.

Let me get some docs up.

[ChrisRackauckas]

Nevermind, I went ahead and did the NoiseProcess stuff. That'll be in this release as well. Docs up soon.

[ChrisRackauckas]

While the tags are waiting to go through, let me give you a brief overview of what this is all looking like. There are some changes to my original plan, but generally the same idea.

Non-Diagonal Noise

This is specified by a keyword argument noise_rate_prototype in the SDEProblem type which specifies the type to build for the noise rate. For example, say we want 4 noise processes for 2 dependent variables, then we can define something like:

f = (t,u,du) -> du.=1.01u
g = function (t,u,du)
  du[1,1] = 0.3u[1]
  du[1,2] = 0.6u[1]
  du[1,3] = 0.9u[1]
  du[1,4] = 0.12u[2]
  du[2,1] = 1.2u[1]
  du[2,2] = 0.2u[2]
  du[2,3] = 0.3u[2]
  du[2,4] = 1.8u[2]
end
prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=zeros(2,4))

and then du*W is the resulting noise. It knows how many noise processes by using the number of columns from noise_rate_prototype. Note that it uses the type you give it. So you can specify du there to be a sparse matrix, for example:

g = function (t,u,du)
  du[1,1] = 0.3u[1]
  du[1,2] = 0.6u[1]
  du[1,3] = 0.9u[1]
  du[1,4] = 0.12u[2]
  du[2,1] = 1.2u[1]
  du[2,2] = 0.2u[2]
  du[2,3] = 0.3u[2]
  du[2,4] = 1.8u[2]
end
prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=sprand(2,4,1.0))

a trivial sparse matrix.

Right now only the EM and the EulerHeun methods are compatible with non-diagonal noise.

Noise Processes

In the SDEProblem type, you can specify the noise. By default, noise=WHITE_NOISE where

WHITE_NOISE = NoiseProcess{:White,true,typeof(white_noise_func_wrapper!)}(white_noise_func_wrapper!)

and

@inline function white_noise_func_wrapper!(rand_vec,integrator)
  wiener_randn!(rand_vec)
end

which is essentially just randn! except it handles complex numbers properly. What this means is that when it wants to generate random numbers, here it is calling wiener_randn!(rand_vec). You can replace this with whatever function you want for generating the random numbers. You can use this to implement the spatial noise. We should probably have a helper function to make making spatial noise earlier.

There is one helper function that exists.

construct_correlated_noisefunc(Γ::AbstractArray)

This constructs a NoiseProcess from a constant Covariance matrix. We probably need more things like this to make it easier, but maybe this helps?

Callbacks

So okay, you have non-diagonal noise being solved, and it's generated by a special NoiseProcess which makes spatially-correlated noise. Now you want to do the projection. This was already implemented before, but let me give you a pointer. Let me conjure up a silly example where we project to the integers after every step. Let's use a linear SDE

using StochasticDiffEq, DiffEqBase

f = (t,u,du) -> du.=1.01u
g = (t,u,du) -> du.=1.01u
u0 = ones(4)
tspan = (0.0,1.0)
prob = SDEProblem(f,g,u0,tspan)

First we have the condition. The callback is called after every step that this is true. Since we want to project after each step, we will just make it:

function condition(t,u,integrator)
  true
end

Now we have our effect. What we want to do is round every value, so:

function affect!(integrator)
  integrator.u .= round.(integrator.u)
end

Now we make the callback:

cb = DiscreteCallback(condition,affect!,save_positions=(false,true))

Notice that we set it to save only after the projection. Now we can look at the cool nonsense that comes out:

sol = solve(prob,EM(),callback=cb,dt=1/10)
using Plots; plot(sol)

Conclusion

Those tools together should be able to solve your problem. If not, I'll need some details and we can try again. The docs should have a lot more details, but feel free to open issues on DiffEqDocs.jl and mention wherever you think the docs are lacking.

[ChrisRackauckas]

Docs are updated but won't build until the tags go through.

[AshtonSBradley]

I notice the docs seem to have updated. The callback stuff works (! super cool !) and I can generate a similar plot.
However, when I run your first code block:

f = (t,u,du) -> du.=1.01u
g = function (t,u,du)
  du[1,1] = 0.3u[1]
  du[1,2] = 0.6u[1]
  du[1,3] = 0.9u[1]
  du[1,4] = 0.12u[2]
  du[2,1] = 1.2u[1]
  du[2,2] = 0.2u[2]
  du[2,3] = 0.3u[2]
  du[2,4] = 1.8u[2]
end
prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=zeros(2,4))

I get

julia> prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=zeros(2,4))
ERROR: MethodError: no method matching DiffEqBase.SDEProblem{uType,tType,isinplace,isinplaceNoise,NoiseClass,F,F2,F3}(::##5#6, ::##7#8, ::Array{Float64,1}, ::Tuple{Float64,Float64}; noise_rate_prototype=[0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0])
Closest candidates are:
  DiffEqBase.SDEProblem{uType,tType,isinplace,isinplaceNoise,NoiseClass,F,F2,F3}(::Any, ::Any, ::Any, ::Any; iip, noise) at /Users/abradley/.julia/v0.5/DiffEqBase/src/problems/sde_problems.jl:20 got unsupported keyword argument "noise_rate_prototype"
  DiffEqBase.SDEProblem{uType,tType,isinplace,isinplaceNoise,NoiseClass,F,F2,F3}{T}(::Any) at sysimg.jl:53 got unsupported keyword argument "noise_rate_prototype"
 in (::Core.#kw#Type)(::Array{Any,1}, ::Type{DiffEqBase.SDEProblem}, ::Function, ::Function, ::Array{Float64,1}, ::Tuple{Float64,Float64}) at ./<missing>:0

Should I be on master here (currently just on latest)?

Another question: I see in the compatibility chart there is no complex number support for ODE. Do you plan to expand complex number support further? It would be great to be able to set up complex variable ODEs and SDEs using ParameterizedFunctions. Can I do this for SDEs yet?

[ChrisRackauckas]

When I run your first code block

It hasn't been released yet. Need to wait on a few things. I'll post when it's released.

Do you plan to expand complex number support further?

Yes, but currently there is no way to calculate Jacobians with complex numbers. So non-stiff methods work, stiff methods won't. Issues for this:

#110
johnmyleswhite/FiniteDiff.jl#14
JuliaDiff/ForwardDiff.jl#157

It would be an easy contribution to those libraries though to set this up, but I haven't had the time to do this.

It would be great to be able to set up ODEs and SDEs using ParameterizedFunctions. Can I do this for SDEs yet?

Yes, just define f and g using ParameterizedFunctions. I don't have a way that supports non-diagonal noise there though. That needs an idea for the syntax overhaul:

SciML/ParameterizedFunctions.jl#17

[ChrisRackauckas]

All that was proposed here is now implemented. Please see the updated docs. To keep things compartmentalized, open a new issue for whatever comes up. Thanks for the discussion and the ideas!

[ChrisRackauckas]

You should be able to give it all a try now. The new docs have built. Let me know if you have any questions!

[AshtonSBradley]

great! thanks so much for this, and for great discussions. I am working on getting some minimal examples running within the new API. More soon!