SDEProblem with OrnsteinUhlenbeckProcess! does not work for scalar state

[tkf]

This works:

noise = OrnsteinUhlenbeckProcess!(1.0, 0.0, 1.0, 0.0, [0.0], [0.0])
prob = SDEProblem((t, u) -> -u, (t, u) -> 1.0, [0.0], (0.0, 1.0), noise=noise)
sol = solve(prob, dt=1e-3, adaptive=false)

But this doesn't (notice W0, Z0 and u0 are scalar this time):

noise = OrnsteinUhlenbeckProcess!(1.0, 0.0, 1.0, 0.0, 0.0, 0.0)
prob = SDEProblem((t, u) -> -u, (t, u) -> 1.0, 0.0, (0.0, 1.0), noise=noise)
sol = solve(prob, dt=1e-3, adaptive=false)

Here is the whole stack trace:

ERROR: MethodError: no method matching randn!(::RandomNumbers.Xorshifts.Xoroshiro128Plus, ::Float64)
Closest candidates are:
  randn!(::AbstractRNG, ::SA<:StaticArrays.StaticArray) where SA<:StaticArrays.StaticArray at ~/.julia/v0.6/StaticArrays/src/arraymath.jl:133
  randn!(::AbstractRNG, ::AbstractArray{T,N} where N) where T at random.jl:1283
Stacktrace:
 [1] (::DiffEqNoiseProcess.OrnsteinUhlenbeck!{Float64,Float64,Float64})(::Float64, ::DiffEqNoiseProcess.NoiseProcess{Float64,1,Float64,Float64,Float64,Array{Float64,1},DiffEqNoiseProcess.OrnsteinUhlenbeck!{Float64,Float64,Float64},Void,true,DataStructures.Stack{Tuple{Float64,Float64,Float64}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Float64}},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus}, ::Float64, ::RandomNumbers.Xorshifts.Xoroshiro128Plus) at ~/.julia/v0.6/DiffEqNoiseProcess/src/ornstein_uhlenbeck.jl:45
 [2] setup_next_step! at ~/.julia/v0.6/DiffEqNoiseProcess/src/noise_interfaces/noise_process_interface.jl:154 [inlined]
 [3] #init#59(::Float64, ::Int64, ::Bool, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Void, ::Bool, ::Void, ::Bool, ::Bool, ::Bool, ::Bool, ::Rational{Int64}, ::Float64, ::Float64, ::Rational{Int64}, ::Rational{Int64}, ::Int64, ::Rational{Int64}, ::Bool, ::Rational{Int64}, ::Rational{Int64}, ::Rational{Int64}, ::Float64, ::Float64, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK, ::DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::Bool, ::DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE, ::String, ::Void, ::Void, ::Bool, ::Bool, ::Bool, ::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.SDEProblem{Float64,Float64,false,DiffEqNoiseProcess.NoiseProcess{Float64,1,Float64,Float64,Float64,Array{Float64,1},DiffEqNoiseProcess.OrnsteinUhlenbeck!{Float64,Float64,Float64},Void,true,DataStructures.Stack{Tuple{Float64,Float64,Float64}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Float64}},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},##89#91,##90#92,Void,UniformScaling{Int64},Void}, ::StochasticDiffEq.SRIW1, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ~/.julia/v0.6/StochasticDiffEq/src/solve.jl:377
 [4] (::DiffEqBase.#kw##init)(::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.SDEProblem{Float64,Float64,false,DiffEqNoiseProcess.NoiseProcess{Float64,1,Float64,Float64,Float64,Array{Float64,1},DiffEqNoiseProcess.OrnsteinUhlenbeck!{Float64,Float64,Float64},Void,true,DataStructures.Stack{Tuple{Float64,Float64,Float64}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Float64}},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},##89#91,##90#92,Void,UniformScaling{Int64},Void}, ::StochasticDiffEq.SRIW1, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0
 [5] #solve#58(::Array{Any,1}, ::Function, ::DiffEqBase.SDEProblem{Float64,Float64,false,DiffEqNoiseProcess.NoiseProcess{Float64,1,Float64,Float64,Float64,Array{Float64,1},DiffEqNoiseProcess.OrnsteinUhlenbeck!{Float64,Float64,Float64},Void,true,DataStructures.Stack{Tuple{Float64,Float64,Float64}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Float64}},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},##89#91,##90#92,Void,UniformScaling{Int64},Void}, ::StochasticDiffEq.SRIW1, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ~/.julia/v0.6/StochasticDiffEq/src/solve.jl:6
 [6] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.SDEProblem{Float64,Float64,false,DiffEqNoiseProcess.NoiseProcess{Float64,1,Float64,Float64,Float64,Array{Float64,1},DiffEqNoiseProcess.OrnsteinUhlenbeck!{Float64,Float64,Float64},Void,true,DataStructures.Stack{Tuple{Float64,Float64,Float64}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Float64}},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},##89#91,##90#92,Void,UniformScaling{Int64},Void}, ::StochasticDiffEq.SRIW1, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0 (repeats 2 times)
 [7] #solve#2(::Bool, ::Array{Any,1}, ::Function, ::DiffEqBase.SDEProblem{Float64,Float64,false,DiffEqNoiseProcess.NoiseProcess{Float64,1,Float64,Float64,Float64,Array{Float64,1},DiffEqNoiseProcess.OrnsteinUhlenbeck!{Float64,Float64,Float64},Void,true,DataStructures.Stack{Tuple{Float64,Float64,Float64}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Float64}},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},##89#91,##90#92,Void,UniformScaling{Int64},Void}, ::Void) at ~/.julia/v0.6/DifferentialEquations/src/default_solve.jl:14
 [8] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.SDEProblem{Float64,Float64,false,DiffEqNoiseProcess.NoiseProcess{Float64,1,Float64,Float64,Float64,Array{Float64,1},DiffEqNoiseProcess.OrnsteinUhlenbeck!{Float64,Float64,Float64},Void,true,DataStructures.Stack{Tuple{Float64,Float64,Float64}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Float64}},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},##89#91,##90#92,Void,UniformScaling{Int64},Void}, ::Void) at ./<missing>:0
 [9] #solve#1(::Bool, ::Array{Any,1}, ::Function, ::DiffEqBase.SDEProblem{Float64,Float64,false,DiffEqNoiseProcess.NoiseProcess{Float64,1,Float64,Float64,Float64,Array{Float64,1},DiffEqNoiseProcess.OrnsteinUhlenbeck!{Float64,Float64,Float64},Void,true,DataStructures.Stack{Tuple{Float64,Float64,Float64}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Float64}},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},##89#91,##90#92,Void,UniformScaling{Int64},Void}) at ~/.julia/v0.6/DifferentialEquations/src/default_solve.jl:5
 [10] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.SDEProblem{Float64,Float64,false,DiffEqNoiseProcess.NoiseProcess{Float64,1,Float64,Float64,Float64,Array{Float64,1},DiffEqNoiseProcess.OrnsteinUhlenbeck!{Float64,Float64,Float64},Void,true,DataStructures.Stack{Tuple{Float64,Float64,Float64}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Float64}},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},##89#91,##90#92,Void,UniformScaling{Int64},Void}) at ./<missing>:0
 [11] macro expansion at ./REPL.jl:97 [inlined]
 [12] (::Base.REPL.##1#2{Base.REPL.REPLBackend})() at ./event.jl:73

Tested against current master (d424fcd). DiffEqNoiseProcess.jl is v0.5.0.

I found the issue 29 but not sure if it is related.

[ChrisRackauckas]

The in place version will always require that the state is a mutable type. For scalars, use the out-of-place version:

noise = OrnsteinUhlenbeckProcess(1.0, 0.0, 1.0, 0.0, 0.0, 0.0)

Note that for this to be correct though, I'm pretty sure you're going to need to use Euler-Maruyama:

sol = solve(prob, EM(),dt=1e-3, adaptive=false)

The higher order methods will probably converge with only order 0.5 with colored noise. I'm not sure the
Wiktorsson approximations still hold when it's not white noise. I should probably find a way to make sure the defaults take note of that and throw an error otherwise.

[tkf]

Thanks. I should have tried the out-of-place version. Though I think it'd be better if NoiseProcess constructor raises an error when I try to construct inplace version with immutable types.

And also thanks for the advice on the method. I wasn't paying attention to the method at all at this point since I wanted to make the code run in the first place (so that I can see what the sol.alg is).

[ChrisRackauckas]

Though I think it'd be better if NoiseProcess constructor raises an error when I try to construct inplace version with immutable types.

Julia doesn't have a compile-time trait for mutability. Hopefully it's a 1.x thing, then we can throw a good error for it.

JuliaLang/julia#21869

[tkf]

I know it's not cool but can't you do that at run time in the constructor manually?

[ChrisRackauckas]

Not in a way that's always correct. A good example is here:

SciML/RecursiveArrayTools.jl#19

Many times a mutable array can be wrapped in an immutable wrapper and setindex! is then defined on that. I think a GPUArray is actually a useful example which ends up being mutable if checked incorrectly.

[tkf]

I see. Thanks.