Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Setup DiffEqGPU integration #122

Merged
merged 4 commits into from
Oct 20, 2023
Merged

[WIP] Setup DiffEqGPU integration #122

merged 4 commits into from
Oct 20, 2023

Conversation

ChrisRackauckas
Copy link
Member

MWE:

from diffeqpy import de

def f(u,p,t):
    x, y, z = u
    sigma, rho, beta = p
    return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]

u0 = [1.0,0.0,0.0]
tspan = (0., 100.)
p = [10.0,28.0,8/3]
prob = de.ODEProblem(f, u0, tspan, p)
fast_prob = de.jit(prob)
sol = de.solve(fast_prob,saveat=0.01)

import random
def prob_func(prob,i,rep):
  de.remake(prob,u0=[random.uniform(0, 1)*u0[i] for i in range(0,3)],
            p=[random.uniform(0, 1)*p[i] for i in range(0,3)])

ensembleprob = de.EnsembleProblem(fast_prob, prob_func = prob_func, safetycopy=False)
sol = de.solve(ensembleprob,de.Tsit5(),de.EnsembleSerial(),trajectories=10000,saveat=0.01)

from diffeqpy import cuda

sol = de.solve(ensembleprob,cuda.GPUTsit5(),cuda.EnsembleGPUKernel(cuda.CUDABackend()),trajectories=10000,saveat=0.01)

Currently fails at the first ensemble solve:

sol = de.solve(ensembleprob,de.Tsit5(),de.EnsembleSerial(),trajectories=10000,saveat=0.01)
>>> sol = de.solve(ensembleprob,de.Tsit5(),de.EnsembleSerial(),trajectories=10000,saveat=0.01)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Users\accou\.julia\packages\PythonCall\qTEA1\src\jlwrap\any.jl", line 208, in __call__
    return self._jl_callmethod($(pyjl_methodnum(pyjlany_call)), args, kwargs)
       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
juliacall.JuliaError: MethodError: no method matching init(::Py, ::OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; saveat::Float64)

Closest candidates are:
  init(!Matched::SciMLBase.OptimizationProblem, ::Any, !Matched::Any...; kwargs...)
   @ SciMLBase C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\solve.jl:146
  init(!Matched::SciMLBase.PDEProblem, ::SciMLBase.AbstractDEAlgorithm, !Matched::Any...; kwargs...)
   @ DiffEqBase C:\Users\accou\.julia\packages\DiffEqBase\xSmHR\src\solve.jl:1116
  init(!Matched::SciMLBase.AbstractJumpProblem, ::Any...; kwargs...)
   @ DiffEqBase C:\Users\accou\.julia\packages\DiffEqBase\xSmHR\src\solve.jl:499
  ...

Stacktrace:
  [1] solve(::Py, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ CommonSolve C:\Users\accou\.julia\packages\CommonSolve\JfpfI\src\CommonSolve.jl:23
  [2] batch_func(i::Int64, prob::SciMLBase.EnsembleProblem{SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Py, typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ SciMLBase C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\ensemble\basic_ensemble_solve.jl:100
  [3] (::SciMLBase.var"#604#605"{Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}}, SciMLBase.EnsembleProblem{SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Py, typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}})(i::Int64)
    @ SciMLBase C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\ensemble\basic_ensemble_solve.jl:154
  [4] responsible_map
    @ C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\ensemble\basic_ensemble_solve.jl:147 [inlined]
  [5] #solve_batch#603
    @ C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\ensemble\basic_ensemble_solve.jl:153 [inlined]
  [6] solve_batch
    @ C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\ensemble\basic_ensemble_solve.jl:152 [inlined]
  [7] macro expansion
    @ .\timing.jl:393 [inlined]
  [8] __solve(prob::SciMLBase.EnsembleProblem{SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Py, typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, ensemblealg::SciMLBase.EnsembleSerial; trajectories::Int64, batch_size::Int64, pmap_batch_size::Int64, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ SciMLBase C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\ensemble\basic_ensemble_solve.jl:64
  [9] solve(::SciMLBase.EnsembleProblem{SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Py, typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, ::OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:saveat, :trajectories), Tuple{Float64, Int64}}})
    @ DiffEqBase C:\Users\accou\.julia\packages\DiffEqBase\xSmHR\src\solve.jl:1053
 [10] pyjlany_call(self::typeof(CommonSolve.solve), args_::Py, kwargs_::Py)
    @ PythonCall C:\Users\accou\.julia\packages\PythonCall\qTEA1\src\jlwrap\any.jl:34
 [11] _pyjl_callmethod(f::Any, self_::Ptr{PythonCall.C.PyObject}, args_::Ptr{PythonCall.C.PyObject}, nargs::Int64)
    @ PythonCall C:\Users\accou\.julia\packages\PythonCall\qTEA1\src\jlwrap\base.jl:69
 [12] _pyjl_callmethod(o::Ptr{PythonCall.C.PyObject}, args::Ptr{PythonCall.C.PyObject})
    @ PythonCall.C C:\Users\accou\.julia\packages\PythonCall\qTEA1\src\cpython\jlwrap.jl:47

MWE:

```py
from diffeqpy import de

def f(u,p,t):
    x, y, z = u
    sigma, rho, beta = p
    return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]

u0 = [1.0,0.0,0.0]
tspan = (0., 100.)
p = [10.0,28.0,8/3]
prob = de.ODEProblem(f, u0, tspan, p)
fast_prob = de.jit(prob)
sol = de.solve(fast_prob,saveat=0.01)

import random
def prob_func(prob,i,rep):
  de.remake(prob,u0=[random.uniform(0, 1)*u0[i] for i in range(0,3)],
            p=[random.uniform(0, 1)*p[i] for i in range(0,3)])

ensembleprob = de.EnsembleProblem(fast_prob, prob_func = prob_func, safetycopy=False)
sol = de.solve(ensembleprob,de.Tsit5(),de.EnsembleSerial(),trajectories=10000,saveat=0.01)

from diffeqpy import cuda

sol = de.solve(ensembleprob,cuda.GPUTsit5(),cuda.EnsembleGPUKernel(cuda.CUDABackend()),trajectories=10000,saveat=0.01)
```

Currently fails at the first ensemble solve:

```py
sol = de.solve(ensembleprob,de.Tsit5(),de.EnsembleSerial(),trajectories=10000,saveat=0.01)
```

```
>>> sol = de.solve(ensembleprob,de.Tsit5(),de.EnsembleSerial(),trajectories=10000,saveat=0.01)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Users\accou\.julia\packages\PythonCall\qTEA1\src\jlwrap\any.jl", line 208, in __call__
    return self._jl_callmethod($(pyjl_methodnum(pyjlany_call)), args, kwargs)
       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
juliacall.JuliaError: MethodError: no method matching init(::Py, ::OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; saveat::Float64)

Closest candidates are:
  init(!Matched::SciMLBase.OptimizationProblem, ::Any, !Matched::Any...; kwargs...)
   @ SciMLBase C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\solve.jl:146
  init(!Matched::SciMLBase.PDEProblem, ::SciMLBase.AbstractDEAlgorithm, !Matched::Any...; kwargs...)
   @ DiffEqBase C:\Users\accou\.julia\packages\DiffEqBase\xSmHR\src\solve.jl:1116
  init(!Matched::SciMLBase.AbstractJumpProblem, ::Any...; kwargs...)
   @ DiffEqBase C:\Users\accou\.julia\packages\DiffEqBase\xSmHR\src\solve.jl:499
  ...

Stacktrace:
  [1] solve(::Py, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ CommonSolve C:\Users\accou\.julia\packages\CommonSolve\JfpfI\src\CommonSolve.jl:23
  [2] batch_func(i::Int64, prob::SciMLBase.EnsembleProblem{SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Py, typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ SciMLBase C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\ensemble\basic_ensemble_solve.jl:100
  [3] (::SciMLBase.var"#604#605"{Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}}, SciMLBase.EnsembleProblem{SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Py, typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}})(i::Int64)
    @ SciMLBase C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\ensemble\basic_ensemble_solve.jl:154
  [4] responsible_map
    @ C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\ensemble\basic_ensemble_solve.jl:147 [inlined]
  [5] #solve_batch#603
    @ C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\ensemble\basic_ensemble_solve.jl:153 [inlined]
  [6] solve_batch
    @ C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\ensemble\basic_ensemble_solve.jl:152 [inlined]
  [7] macro expansion
    @ .\timing.jl:393 [inlined]
  [8] __solve(prob::SciMLBase.EnsembleProblem{SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Py, typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, ensemblealg::SciMLBase.EnsembleSerial; trajectories::Int64, batch_size::Int64, pmap_batch_size::Int64, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ SciMLBase C:\Users\accou\.julia\packages\SciMLBase\McEqc\src\ensemble\basic_ensemble_solve.jl:64
  [9] solve(::SciMLBase.EnsembleProblem{SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Py, typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, ::OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:saveat, :trajectories), Tuple{Float64, Int64}}})
    @ DiffEqBase C:\Users\accou\.julia\packages\DiffEqBase\xSmHR\src\solve.jl:1053
 [10] pyjlany_call(self::typeof(CommonSolve.solve), args_::Py, kwargs_::Py)
    @ PythonCall C:\Users\accou\.julia\packages\PythonCall\qTEA1\src\jlwrap\any.jl:34
 [11] _pyjl_callmethod(f::Any, self_::Ptr{PythonCall.C.PyObject}, args_::Ptr{PythonCall.C.PyObject}, nargs::Int64)
    @ PythonCall C:\Users\accou\.julia\packages\PythonCall\qTEA1\src\jlwrap\base.jl:69
 [12] _pyjl_callmethod(o::Ptr{PythonCall.C.PyObject}, args::Ptr{PythonCall.C.PyObject})
    @ PythonCall.C C:\Users\accou\.julia\packages\PythonCall\qTEA1\src\cpython\jlwrap.jl:47
```
@ChrisRackauckas
Copy link
Member Author

@LilithHafner do you know why that conversion of Py -> EnsembleProblem isn't done?

@LilithHafner
Copy link
Member

I think the issue is that the keyword argument prob_func in

ensembleprob = de.EnsembleProblem(fast_prob, prob_func = prob_func, safetycopy=False)`

is not converted.

If I change that line to

ensembleprob = de.EnsembleProblem(fast_prob, safetycopy=False)`

Then the non-cuda solve works (or maybe it hangs, idk. It works when I lower trajectories to 1000).

@ChrisRackauckas
Copy link
Member Author

from diffeqpy import de

def f(u,p,t):
    x, y, z = u
    sigma, rho, beta = p
    return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]

u0 = [1.0,0.0,0.0]
tspan = (0., 100.)
p = [10.0,28.0,8/3]
prob = de.ODEProblem(f, u0, tspan, p)
fast_prob = de.jit(prob)
sol = de.solve(fast_prob,saveat=0.01)

import random
def prob_func(prob,i,rep):
  de.remake(prob,u0=[random.uniform(0, 1)*u0[i] for i in range(0,3)],
            p=[random.uniform(0, 1)*p[i] for i in range(0,3)])

ensembleprob = de.EnsembleProblem(fast_prob, safetycopy=False)
sol = de.solve(ensembleprob,cuda.GPUTsit5(),cuda.EnsembleGPUKernel(CUDABackend()),trajectories=10000,saveat=0.01)
>>> sol = de.solve(ensembleprob,cuda.GPUTsit5(),cuda.EnsembleGPUKernel(CUDABackend()),trajectories=10000,saveat=0.01)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Users\accou\.julia\packages\PythonCall\qTEA1\src\jlwrap\any.jl", line 208, in __call__
    return self._jl_callmethod($(pyjl_methodnum(pyjlany_call)), args, kwargs)
       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
juliacall.JuliaError: CuArray only supports element types that are allocated inline.
Vector{Float64} is a mutable type

Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] check_eltype(T::Type)
    @ CUDA C:\Users\accou\.julia\packages\CUDA\35NC6\src\array.jl:70
  [3] CUDA.CuArray{Vector{Float64}, 2, CUDA.Mem.DeviceBuffer}(#unused#::UndefInitializer, dims::Tuple{Int64, Int64})
    @ CUDA C:\Users\accou\.julia\packages\CUDA\35NC6\src\array.jl:85
  [4] (CUDA.CuArray{Vector{Float64}, 2})(#unused#::UndefInitializer, dims::Tuple{Int64, Int64})
    @ CUDA C:\Users\accou\.julia\packages\CUDA\35NC6\src\array.jl:176
  [5] (CUDA.CuArray{Vector{Float64}})(#unused#::UndefInitializer, dims::Tuple{Int64, Int64})
    @ CUDA C:\Users\accou\.julia\packages\CUDA\35NC6\src\array.jl:187
  [6] allocate(#unused#::CUDA.CUDAKernels.CUDABackend, #unused#::Type{Vector{Float64}}, dims::Tuple{Int64, Int64})
    @ CUDA.CUDAKernels C:\Users\accou\.julia\packages\CUDA\35NC6\src\CUDAKernels.jl:17
  [7] vectorized_asolve(probs::CUDA.CuArray{DiffEqGPU.ImmutableODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, 1, CUDA.Mem.DeviceBuffer}, prob::SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::DiffEqGPU.GPUTsit5; dt::Float32, saveat::Float64, save_everystep::Bool, abstol::Float32, reltol::Float32, debug::Bool, callback::SciMLBase.CallbackSet{Tuple{}, Tuple{}}, tstops::Nothing, kwargs::Base.Pairs{Symbol, DiffEqGPU.var"#114#120", Tuple{Symbol}, NamedTuple{(:unstable_check,), Tuple{DiffEqGPU.var"#114#120"}}})
    @ DiffEqGPU C:\Users\accou\.julia\dev\DiffEqGPU\src\ensemblegpukernel\lowerlevel_solve.jl:178
  [8] batch_solve_up_kernel(ensembleprob::SciMLBase.EnsembleProblem{SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, typeof(SciMLBase.DEFAULT_PROB_FUNC), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, probs::Vector{DiffEqGPU.ImmutableODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, alg::DiffEqGPU.GPUTsit5, ensemblealg::DiffEqGPU.EnsembleGPUKernel{CUDA.CUDAKernels.CUDABackend}, I::UnitRange{Int64}, adaptive::Bool; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:saveat, :unstable_check), Tuple{Float64, DiffEqGPU.var"#114#120"}}})
    @ DiffEqGPU C:\Users\accou\.julia\dev\DiffEqGPU\src\solve.jl:267
  [9] batch_solve(ensembleprob::SciMLBase.EnsembleProblem{SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, typeof(SciMLBase.DEFAULT_PROB_FUNC), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::DiffEqGPU.GPUTsit5, ensemblealg::DiffEqGPU.EnsembleGPUKernel{CUDA.CUDAKernels.CUDABackend}, I::UnitRange{Int64}, adaptive::Bool; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:unstable_check, :saveat), Tuple{DiffEqGPU.var"#114#120", Float64}}})
    @ DiffEqGPU C:\Users\accou\.julia\dev\DiffEqGPU\src\solve.jl:168
 [10] macro expansion
    @ .\timing.jl:393 [inlined]
 [11] __solve(ensembleprob::SciMLBase.EnsembleProblem{SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#545"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x946926fe, 0xab7b8dc2, 0x707dbba2, 0x8b060826, 0xe63fdf5b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xaf9b47a5, 0x8d86fb18, 0x65c64d40, 0xf8480c5a, 0x8614cffa), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, typeof(SciMLBase.DEFAULT_PROB_FUNC), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::DiffEqGPU.GPUTsit5, ensemblealg::DiffEqGPU.EnsembleGPUKernel{CUDA.CUDAKernels.CUDABackend}; trajectories::Int64, batch_size::Int64, unstable_check::Function, adaptive::Bool, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ DiffEqGPU C:\Users\accou\.julia\dev\DiffEqGPU\src\solve.jl:55
 [12] __solve
    @ C:\Users\accou\.julia\dev\DiffEqGPU\src\solve.jl:1 [inlined]
 [13] #solve#44
    @ C:\Users\accou\.julia\packages\DiffEqBase\xSmHR\src\solve.jl:1053 [inlined]
 [14] pyjlany_call(self::typeof(CommonSolve.solve), args_::Py, kwargs_::Py)
    @ PythonCall C:\Users\accou\.julia\packages\PythonCall\qTEA1\src\jlwrap\any.jl:34
 [15] _pyjl_callmethod(f::Any, self_::Ptr{PythonCall.C.PyObject}, args_::Ptr{PythonCall.C.PyObject}, nargs::Int64)
    @ PythonCall C:\Users\accou\.julia\packages\PythonCall\qTEA1\src\jlwrap\base.jl:69
 [16] _pyjl_callmethod(o::Ptr{PythonCall.C.PyObject}, args::Ptr{PythonCall.C.PyObject})
    @ PythonCall.C C:\Users\accou\.julia\packages\PythonCall\qTEA1\src\cpython\jlwrap.jl:47

@utkarsh530 can DiffEqGPU not auto-convert to SA?

@ChrisRackauckas
Copy link
Member Author

from diffeqpy import de

def f(u,p,t):
    x, y, z = u
    sigma, rho, beta = p
    return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]

u0 = [1.0,0.0,0.0]
tspan = (0., 100.)
p = [10.0,28.0,8/3]
prob = de.ODEProblem(f, u0, tspan, p)
fast_prob = de.jit(prob)
sol = de.solve(fast_prob,saveat=0.01)

import random
def prob_func(prob,i,rep):
  de.remake(prob,u0=[random.uniform(0, 1)*u0[i] for i in range(0,3)],
            p=[random.uniform(0, 1)*p[i] for i in range(0,3)])

ensembleprob = de.EnsembleProblem(fast_prob, safetycopy=False)
sol = de.solve(ensembleprob,de.Tsit5(),de.EnsembleSerial(),trajectories=10000,saveat=0.01)

from diffeqpy import cuda
sol = de.solve(ensembleprob,cuda.GPUTsit5(),cuda.EnsembleGPUKernel(CUDABackend()),trajectories=10000,saveat=0.01)

This works with SciML/DiffEqGPU.jl#308

@ChrisRackauckas ChrisRackauckas merged commit a7afab6 into master Oct 20, 2023
2 checks passed
@ChrisRackauckas ChrisRackauckas deleted the cuda branch October 20, 2023 04:04
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants