sciml / scimlbase.jl Goto Github PK

View Code? Open in Web Editor NEW

121.0 13.0 91.0 14.09 MB

The Base interface of the SciML ecosystem

Home Page: https://docs.sciml.ai/SciMLBase/stable

License: MIT License

Julia 100.00%

julia ode sciml dae dde differentialequations ordinary-differential-equations scientific-machine-learning sde

scimlbase.jl's Introduction

SciMLBase

SciMLBase.jl is the core interface definition of the SciML ecosystem. It is a low dependency library made to be depended on by the downstream libraries to supply the common interface and allow for interexchange of mathematical problems.

v2.0 Breaking Changes

The breaking changes in v2.0 are:

IntegralProblem has moved to an interface with IntegralFunction and BatchedIntegralFunction which requires specifying prototypes for the values to be modified instead of nout and batch. #497
ODEProblem was made temporarily into a mutable struct to allow for EnzymeRules support. Using the mutation throws a warning that this is only experimental and should not be relied on. #501
BVProblem now has a new interface for TwoPointBVProblem which splits the bc terms for the two sides, forcing a true two-point BVProblem to allow for further specializations and to allow for wrapping Fortran solvers in the interface. #477
SDEProblem constructor was changed to remove an anti-pattern which required passing the diffusion function g twice, i.e. SDEProblem(SDEFunction(f,g),g, ...). Now this is simply SDEProblem(SDEFunction(f,g),...). #489

scimlbase.jl's People

Contributors

Stargazers

Watchers

Forkers

stjordanis utkarsh530 wi11dey kirillzubov giordano fonsp standardgalactic isaacsas briochemc sebastianm-c andreasnoack chriselrod aml5600 barroff garibarba lamorton sharanry lungd avik-pal nathompson cwiese dannys4 baggepinnen vpuri3 keno nathanaelbosch zoemcc vboussange ggkountouras alyst xtalax brainsmaker frankschae aayushsabharwal daniel-vainsencher yichengdwu rmsrosa oscardssmith wsphillips ranocha hustf thomvet valentinkaisermayer harisorgn arnostrouwen anandijain manuelbb-upb rafaelarutjunjan xlxs4 ven-k gaurav-arya charleskawczynski varlad uncomfyhalomacro jonlym erikqqy torkele ignace-computing yash2798 schmidtjonathan kristofferc thummeto kalelr staticfloat csimal datseris dhairyalgandhi contradict prbzrg lxvm ilianpihlajamaa wsmoses danielvandh lilithhafner pepijndevos qfl3x hadriennu yonatanwesen arnebouillon rzyu45 asinghvi17 h-sax paraspuneetsingh hippyhippohops visr gdalle cedareda termi-official jonathanfischer97 m-bossart abhro

scimlbase.jl's Issues

Segmentation fault with `get_du`

I'm getting segmentation fault when get_du is called on a DiscreteProblem integrator.

using DifferentialEquations

f(u,p,t) = [1]
prob=DiscreteProblem(f, [1], (1., 2.))
i = init(prob, FunctionMap())
get_du(i)

plot(NoiseProcess) direct solve

I noticed that direct solve of NoiseProcess doesn't give a NoiseSolution (doesn't exist), but a NoiseProcess, so the plot recipe doesn't work. It should be really easy to fix as plot(sol.t, sol.u) works. I'm not sure where to add this though.

It looks like the other solutions dispatch on SciMLBase.AbstractTimeseriesSolution, but NoiseProcess does not subtype this

I realize this is probably the wrong place to open, but I found the recipes in here.

μ = 1.0
σ = 2.0
W = GeometricBrownianMotionProcess(μ,σ,0.0,1.0,1.0) # empty plot
prob = NoiseProblem(W,(0.0,1.0))
sol2 = solve(prob;dt=0.1)
plot(sol) # fails
plot(sol.t, sol.u) # works

make solution indexing work with top level system namespacing

If I have a ModelingToolkit ODESystem or such, say sys, with a variable x, then

@variables t, x(t)
... # solve an ODESystem with x
sol[x]

works nicely, however,

sol[sys.x]

doesn't appear to work.

Is it possible to make this work too? This would be convenient for users who don't have the variable x already loaded, say because they are using a generated system.

A similar comment would apply to plotting.

Need to turn off threads in `RecursiveFactorization` if `EnsembleThreads()`

Can potentially stall.

potential logic error in solution interface

It looks to me like this is a mistake:

  if i == nothing
      if issymbollike(i) && Symbol(i) == getindepsym(A)

If i is nothing then it makes no sense to test whether it is symbol-like or is equal to the independent variable -- I think it's probably meant to be sym here instead of i.

remake SecondOrderODEProblem with separate u0 and du0

It was not immediate to me how to remake a SecondOrderODEProblem.

It seemed natural to me to try remake(prob, u0 = [2.0], du0 = [3.0]), but finally @ChrisRackauckas told me I should use remake(prob, u0 = ArrayPartition([2.0], [3.0])).

I was not familiar with ArrayPartition, but it still seems odd to me that u0 = ArrayPartition([2.0], [3.0]) sets both u0 and du0. Anyways, it could be interesting to allow remake(prob, u0 = [2.0], du0 = [3.0]), as well, and possibly rename the keyword for ArrayPartition([2.0], [3.0]), although that would probably be a breaking change.

Chris suggested me to open this issue.

Symbolically index into an EnsembleSolution

With that PR I can index into an ODESolution using a Num
The following code raises an BoundsError: attempt to access Tuple{Int64} at index [0]

using ModelingToolkit
using ModelingToolkit: get_defaults
using OrdinaryDiffEq
using DiffEqSensitivity
using Zygote

Zygote.@adjoint function Base.getindex(VA::ODESolution, sym, j::Int)
    function ODESolution_getindex_pullback(Δ)
        i = SciMLBase.issymbollike(sym) ? SciMLBase.sym_to_index(sym, VA) : sym
        if i === nothing
            getter = SciMLBase.getobserved(VA)
            grz = Zygote.pullback(getter, sym, VA.u[j], VA.prob.p, VA.t[j])[2](Δ)
            du = [k == j ? grz[2] : zero(VA.u[1]) for k in 1:length(VA.u)]
            dp = grz[3] # pullback for p
			dprob = remake(VA.prob, p = dp)
            T = eltype(eltype(VA.u))
            N = length(VA.prob.p)
            Δ′ = ODESolution{T, N, typeof(du), Nothing, Nothing, Nothing, Nothing,
                             typeof(dprob), Nothing, Nothing, Nothing}(du, nothing,
                              nothing, nothing, nothing, dprob, nothing, nothing,
                              VA.dense, 0, nothing, VA.retcode)
			(Δ′, nothing, nothing)
        else
            Δ′ = [m == j ? [i == k ? Δ : zero(VA.u[1][1]) for k in 1:length(VA.u[1])] : zero(VA.u[1]) for m in 1:length(VA.u)]
            (Δ′, nothing, nothing)
        end
    end
    VA[sym, j], ODESolution_getindex_pullback
end

Zygote.@adjoint function Base.getindex(VA::ODESolution, sym)
    function ODESolution_getindex_pullback(Δ)
        i = SciMLBase.issymbollike(sym) ? SciMLBase.sym_to_index(sym, VA) : sym
        if i === nothing
            throw("Zygote AD of purely-symbolic slicing for observed quantities is not yet supported. Work around this by using `A[sym,i]` to access each element sequentially in the function being differentiated.")
        else
            Δ′ = [ [i == k ? Δ[j] : zero(x[1]) for k in 1:length(x)] for (x, j) in zip(VA.u, 1:length(VA))]
            (Δ′, nothing)
        end
    end
    VA[sym], ODESolution_getindex_pullback
end



@variables t
D = Differential(t)

function Neuron(;name)
    @variables v(t)=0.0 I(t)
	@parameters s=0.4
    eqs = D(v) ~ s+I
    ODESystem([eqs],t,[v,I],[s];name)
end

function Network(;name)
    @named n1 = Neuron()
    eqs = n1.I ~ n1.v
    ODESystem([eqs],t,[],[];name,systems=[n1])
end


@named net = Network()
sys = structural_simplify(net)
prob = ODEProblem(sys, ModelingToolkit.get_defaults(sys), (0.0, 1.0), sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)),)

p = [0.13]

function test_loss_ensemble(p, prob)
	prob_func(prob,i,repeat) = remake(prob, p=p)
	output_func(sol,i) = sol, false
	ensemble_prob = EnsembleProblem(prob; prob_func, output_func)
    sol = solve(ensemble_prob, Tsit5(), EnsembleThreads(), trajectories=10)
	@nonamespace n1 = net.n1
	I = [[sol[i][n1.I, t] for t in 1:size(sol[i],2)] for i in 1:size(sol,3)]
	loss=sum(sum.(I))
	return loss
end

test_loss_ensemble(p,prob)
@time gs = Zygote.gradient(p) do p
	test_loss_ensemble(p,prob)
end

BoundsError: attempt to access Tuple{Int64} at index [0]
in eval at base/boot.jl:360 
in top-level scope at base/timing.jl:210
in gradient at Zygote/TaBlo/src/compiler/interface.jl:76
in  at Zygote/TaBlo/src/compiler/interface.jl:41
in  at Zygote/TaBlo/src/compiler/interface2.jl
in Pullback at julia/LTC.jl/src/symbolic_index_issue0.jl:84 
in  at Zygote/TaBlo/src/compiler/interface2.jl
in Pullback at julia/LTC.jl/src/symbolic_index_issue0.jl:75 
in  at Zygote/TaBlo/src/compiler/interface2.jl
in Pullback at DiffEqBase/oe7VF/src/solve.jl:93 
in #1746#back at ZygoteRules/OjfTt/src/adjoint.jl:59 
in  at Zygote/TaBlo/src/lib/lib.jl:203
in  at Zygote/TaBlo/src/compiler/interface2.jl
in Pullback at DiffEqBase/oe7VF/src/solve.jl:96 
in  at ZygoteRules/OjfTt/src/adjoint.jl:59
in #209 at Zygote/TaBlo/src/lib/lib.jl:203 
in  at Zygote/TaBlo/src/compiler/interface2.jl
in Pullback at SciMLBase/oTP8b/src/ensemble/basic_ensemble_solve.jl:103 
in  at Zygote/TaBlo/src/compiler/interface2.jl
in Pullback at SciMLBase/oTP8b/src/ensemble/basic_ensemble_solve.jl:110 
in  at ZygoteRules/OjfTt/src/adjoint.jl:59
in  at DiffEqBase/oe7VF/src/chainrules.jl:62
in collect at base/array.jl:678 
in iterate at base/generator.jl:47 
in  at base/none
in getindex at base/tuple.jl:29

Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.

I also described the problem here
I have troubles with fitting simple data
The code:

using DiffEqFlux
using OrdinaryDiffEq, Flux, Optim, Plots
using Flux, OrdinaryDiffEq
using Zygote 
using DiffEqSensitivity # for ZygteVJP?


abstract type NeuralDELayer <: Function end
basic_tgrad(u,p,t) = zero(u)

struct LTC{M,P,RE,T,TA,AB,A,K} <: NeuralDELayer
    model::M
    p::P # weights
    p_len::Int # for assignment 
    re::RE
    tspan::T
    τ::TA # weights
    τ_len::Int #
    A::AB # weights
    A_len::Int
    args::A
    kwargs::K

    function LTC(model,tspan, τ, A, args...;p = nothing,kwargs...)
        _p,re = Flux.destructure(model) # is it like [p;τ;A] already? 
        if p === nothing
            p = _p
        end
        new{typeof(model),typeof(p),typeof(re),
            typeof(tspan), typeof(τ), typeof(A),
            typeof(args),typeof(kwargs)}(
            model,p, length(p), re,tspan,τ,length(τ),A,length(A),args,kwargs)
    end
end

function (n::LTC)(x)
    function dudt_(u, p, t)
       p_ = @view p[1:n.p_len]
       τ_ = @view p[n.p_len+1:n.p_len+n.τ_len]
       τ_ = Flux.softplus.(τ_) # to ensure τ>=0
       A_ = @view p[n.p_len+n.τ_len+1:end]
       h = -(1 ./τ_+ n.re(p_)(u)) .* u +  n.re(p_)(u) .* A_
    end
    ff = ODEFunction{false}(dudt_,tgrad=basic_tgrad) 
    prob = ODEProblem{false}(ff,x,getfield(n,:tspan), [n.p; n.τ; n.A]) # inital conditions and tspan, etc
    sense = InterpolatingAdjoint(autojacvec=ZygoteVJP()) 
    solve(prob,n.args...;sense=sense,n.kwargs...)
end

Flux.trainable(m::LTC) = (m.p, m.τ, m.A)


# Example 
u0 = Float32[2.; 0.]
datasize = 30
tspan = (0.0f0, 3.0f0)

function trueODEfunc(du,u,p,t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
t = range(tspan[1],tspan[2],length=datasize)
prob = ODEProblem(trueODEfunc,u0,tspan)
ode_data = Array(solve(prob,Tsit5(),saveat=t))
τ = rand(2) # ~size(out) 
A = zeros(2) # ~size(out)

dudt = Chain(x -> x.^3,
             Dense(2,10,tanh),
             Dense(10,2))
ps = Flux.params(n_ode)

pred = n_ode(u0) # Get the prediction using the correct initial condition
scatter(t,ode_data[1,:],label="data")
scatter!(t,pred[1,:],label="prediction")

function predict_n_ode()
  n_ode(u0)
end
loss_n_ode() = sum(abs2,ode_data .- predict_n_ode())


data = Iterators.repeated((), 1000)
opt = ADAM(0.01)
cb = function () #callback function to observe training
    nothing
end
Flux.train!(loss_n_ode, ps, data, opt, cb = cb)

The warning:

┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /home/solar/.julia/packages/SciMLBase/HbD6U/src/integrator_interface.jl:345
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /home/solar/.julia/packages/SciMLBase/HbD6U/src/integrator_interface.jl:345

Interpolating intergrators returns a `Vector{Vector{T}}` not a `DiffEqArray`

Copying an example from test/downstream/integrator_indexing.jl

using ModelingToolkit, OrdinaryDiffEq, RecursiveArrayTools, Test

@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

eqs = [D(x) ~ σ * (y - x),
    D(y) ~ x * (ρ - z) - y,
    D(z) ~ x * y - β * z]

@named lorenz1 = ODESystem(eqs)
@named lorenz2 = ODESystem(eqs)

@parameters γ
@variables a(t) α(t)
connections = [0 ~ lorenz1.x + lorenz2.y + a * γ,
    α ~ 2lorenz1.x + a * γ]
@named sys = ODESystem(connections, t, [a, α], [γ], systems = [lorenz1, lorenz2])
sys_simplified = structural_simplify(sys)

u0 = [lorenz1.x => 1.0,
    lorenz1.y => 0.0,
    lorenz1.z => 0.0,
    lorenz2.x => 0.0,
    lorenz2.y => 1.0,
    lorenz2.z => 0.0,
    a => 2.0]

p = [lorenz1.σ => 10.0,
    lorenz1.ρ => 28.0,
    lorenz1.β => 8 / 3,
    lorenz2.σ => 10.0,
    lorenz2.ρ => 28.0,
    lorenz2.β => 8 / 3,
    γ => 2.0]

tspan = (0.0, 100.0)
prob = ODEProblem(sys_simplified, u0, tspan, p)
integrator = init(prob, Rodas4())
step!(integrator, 100.0, true)

interpolated_integrator = integrator(0.0:1.0:10.0)
@test interpolated_integrator[α] isa Vector   # won't pass because interpolated_integrator is a nested vector

Ambiguity error for `DynamicalDDEProblem`

Running

using Test
using SciMLBase

problems = detect_ambiguities(SciMLBase)

gives

8-element Vector{Tuple{Method, Method}}:
 (DynamicalDDEProblem(f::DynamicalDDEFunction, h, tspan) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:77, DynamicalDDEProblem(f1::Function, f2::Function, args...; kwargs...) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:80)
 (DynamicalDDEProblem(f::DynamicalDDEFunction, h, tspan, p; kwargs...) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:77, DynamicalDDEProblem(f1::Function, f2::Function, args...; kwargs...) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:80)
 ((var"#s255"::Core.var"#Type##kw")(::Any, ::Type{DynamicalDDEProblem}, f::DynamicalDDEFunction, h, tspan, p) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:77, (::Core.var"#Type##kw")(::Any, ::Type{DynamicalDDEProblem}, f1::Function, f2::Function, args...) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:80)
 (DynamicalDDEProblem(f::DynamicalDDEFunction, v0, u0, h, tspan, p; dependent_lags, kwargs...) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:71, DynamicalDDEProblem(f1::Function, f2::Function, args...; kwargs...) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:80)
 (DynamicalDDEProblem(f::DynamicalDDEFunction, v0, u0, h, tspan) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:71, DynamicalDDEProblem(f1::Function, f2::Function, args...; kwargs...) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:80)
 ((var"#s255"::Core.var"#Type##kw")(::Any, ::Type{DynamicalDDEProblem}, f::DynamicalDDEFunction, v0, u0, h, tspan) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:71, (::Core.var"#Type##kw")(::Any, ::Type{DynamicalDDEProblem}, f1::Function, f2::Function, args...) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:80)
 ((var"#s255"::Core.var"#Type##kw")(::Any, ::Type{DynamicalDDEProblem}, f::DynamicalDDEFunction, v0, u0, h, tspan, p) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:71, (::Core.var"#Type##kw")(::Any, ::Type{DynamicalDDEProblem}, f1::Function, f2::Function, args...) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:80)
 ((var"#s255"::Core.var"#Type##kw")(::Any, ::Type{DynamicalDDEProblem}, f::DynamicalDDEFunction, h, tspan) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:77, (::Core.var"#Type##kw")(::Any, ::Type{DynamicalDDEProblem}, f1::Function, f2::Function, args...) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:80)

The above ambiguities refer only to two function definitions:
https://github.com/SciML/SciMLBase.jl/blob/v1.9.0/src/problems/dde_problems.jl#L77-L82

As a MWE:

using SciMLBase

f1(u,t,h) = u
f2(v,h) = v
f = DynamicalDDEFunction(f1, f2)
h(p,t) = (1,1)

DynamicalDDEProblem(f, h, (0,1))

which gives

ERROR: MethodError: DynamicalDDEProblem(::DynamicalDDEFunction{false, DDEFunction{false, typeof(f1), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, DDEFunction{false, typeof(f2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, ::typeof(h), ::Tuple{Int64, Int64}) is ambiguous. Candidates:
  DynamicalDDEProblem(f1::Function, f2::Function, args...; kwargs...) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:80
  DynamicalDDEProblem(f::DynamicalDDEFunction, h, tspan) in SciMLBase at C:\Users\sebastian\.julia\dev\SciMLBase\src\problems\dde_problems.jl:77
Possible fix, define
  DynamicalDDEProblem(::DynamicalDDEFunction, ::Function, ::Any)
Stacktrace:
 [1] top-level scope

The problem appears due to the fact that DynamicalDDEFunction <: Function and the commit c555ead removed the restriction on the h argument and thus in both

function DynamicalDDEProblem(f::DynamicalDDEFunction,h,tspan,p=NullParameters();kwargs...)

and

function DynamicalDDEProblem(f1::Function,f2::Function,args...;kwargs...)

the first argument fits the method, but for the second argument we have Any in the first case (which was more specific) and Function in the second (which is more specific than Any).

One option would be to revert the suggestion by @devmotion in c555ead
The other would be to replace

function DynamicalDDEProblem(f1::Function,f2::Function,args...;kwargs...)

with two separate methods, but this would result in code duplication as mentioned in #14 (comment)

I can open a PR to fix this if you would like, just let me know what approach do you think is best.

Allow symbols for variable choices in interpolations

We now allow sol[x] at time points, but it would be to allow sol(t,save_idxs=x) and sol(t,save_idxs=[x,y]) as well. In addition sol(t, save_idxs = AllObserved()) would be very helpful in some cases too.

kwargs in Sundials throwing Warning

In the latest release from DiffEqBase a warning about future allowed kwargs don't include

advance_to_tstop   
stop_at_next_tstop

Which are used in Sundials.jl

Downstream tests fail

Structural simplify fails in downstream tests with ERROR: LoadError: MethodError: no method matching operation(::Sym{Real, Nothing}).

CC: @ChrisRackauckas

julia> include("test\\downstream\\observables.jl")
[ Info: Precompiling ModelingToolkit [961ee093-0014-501f-94e3-6117800e7a78]
[ Info: Precompiling OrdinaryDiffEq [1dea7af3-3e70-54e6-95c3-0bf5283fa5ed]
ERROR: LoadError: MethodError: no method matching operation(::Sym{Real, Nothing})
Closest candidates are:
  operation(::Term) at C:\Users\Sharan Yalburgi\.julia\packages\SymbolicUtils\JmtMa\src\types.jl:346
  operation(::SymbolicUtils.Add) at C:\Users\Sharan Yalburgi\.julia\packages\SymbolicUtils\JmtMa\src\types.jl:625
  operation(::SymbolicUtils.Mul) at C:\Users\Sharan Yalburgi\.julia\packages\SymbolicUtils\JmtMa\src\types.jl:768
  ...
Stacktrace:
 [1] (::ModelingToolkit.StructuralTransformations.var"#14#20")(x::Sym{Real, Nothing})
   @ ModelingToolkit.StructuralTransformations C:\Users\Sharan Yalburgi\.julia\packages\ModelingToolkit\XlBco\src\structural_transformation\pantelides.jl:64
 [2] filter(f::ModelingToolkit.StructuralTransformations.var"#14#20", a::Vector{Any})
   @ Base .\array.jl:2507
 [3] pantelides_reassemble(sys::ODESystem, eqassoc::Vector{Int64}, assign::Vector{Int64})
   @ ModelingToolkit.StructuralTransformations C:\Users\Sharan Yalburgi\.julia\packages\ModelingToolkit\XlBco\src\structural_transformation\pantelides.jl:64
 [4] dae_index_lowering(sys::ODESystem; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit.StructuralTransformations C:\Users\Sharan Yalburgi\.julia\packages\ModelingToolkit\XlBco\src\structural_transformation\pantelides.jl:153
 [5] dae_index_lowering(sys::ODESystem)
   @ ModelingToolkit.StructuralTransformations C:\Users\Sharan Yalburgi\.julia\packages\ModelingToolkit\XlBco\src\structural_transformation\pantelides.jl:150
 [6] structural_simplify(sys::ODESystem)
   @ ModelingToolkit C:\Users\Sharan Yalburgi\.julia\packages\ModelingToolkit\XlBco\src\systems\abstractsystem.jl:645
 [7] top-level scope
   @ C:\Users\Sharan Yalburgi\.julia\dev\test\SciMLBase.jl\test\downstream\observables.jl:18
 [8] include(fname::String)
   @ Base.MainInclude .\client.jl:444
 [9] top-level scope
   @ REPL[8]:1
in expression starting at C:\Users\Sharan Yalburgi\.julia\dev\test\SciMLBase.jl\test\downstream\observables.jl:18

Different AD backends for objective and constraints

While the objective function is typically scalar, constraint functions may have a very large output arity. There may thus be situations where it's beneficial to differentiate the objective in reverse, but calculate the constraint jacobian in forward mode.

Possibly extend the observed interface to allow for chunked calculations?

The observed interface allows for observed(sym,u,p,t) at a single state point and a single symbol at a time. If this becomes a computational issue, we can extend the interface in two directions:

allobserved(u,p,t) for calculating all of the observed (and state) variables, returning the big array of all things.
multipoint_observed(sym,u,p,t) which allows for checking the symbol once and applying it at many t states simultaneously. This is like what's done in the interpolation because that can specialize the computation a bunch. But is it necessary for speed here?

If we do one and two, we will need multipoint_alloberved(u,p,t) as well. Also, (1) needs a user interface, like allobserved(sol) or allobserved(sol,1:5) (for the first 5 time points).

@YingboMa do you think either of these are necessary?

`calculate_solution_errors!` keeps growing/pushing to `u_analytic`

From the example in the documentation

using StochasticDiffEq, SciMLBase
α=1
β=1
u₀=1/2
f(u,p,t) = α*u
g(u,p,t) = β*u
dt = 1//2^(4)
tspan = (0.0,1.0)
f_analytic(u₀,p,t,W) = u₀*exp((α-(β^2)/2)*t+β*W)
ff = SDEFunction(f,g,analytic=f_analytic)
prob = SDEProblem(ff,g,u₀,(0.0,1.0))

The first call to sol = solve(prob, ...) fills sol.u_analytic to the same length as sol.u, but subsequent calls to calculate_solution_errors!(prob) keep doubling the length of sol.u_analytic.

julia> sol = solve(prob,EM(),dt=dt);

julia> length(sol.u), length(sol.u_analytic)
(17, 17)

julia> SciMLBase.calculate_solution_errors!(sol);
ERROR: DimensionMismatch("dimensions must match: a has dims (Base.OneTo(17),), b has dims (Base.OneTo(34),), mismatch at 1")
Stacktrace:
 [1] promote_shape
   @ ./indices.jl:178 [inlined]
 [2] promote_shape
   @ ./indices.jl:169 [inlined]
 [3] -(A::Vector{Float64}, B::Vector{Float64})
   @ Base ./arraymath.jl:38
 [4] calculate_solution_errors!(sol::RODESolution{Float64, 1, Vector{Float64}, Vector{Float64}, Dict{Symbol, Float64}, Vector{Float64}, NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Nothing}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Nothing}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, Nothing, SDEFunction{false, typeof(f), typeof(g), LinearAlgebra.UniformScaling{Bool}, typeof(f_analytic), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, typeof(g), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Nothing}, EM{true}, StochasticDiffEq.LinearInterpolationData{Vector{Float64}, Vector{Float64}}, DiffEqBase.DEStats}; fill_uanalytic::Bool, timeseries_errors::Bool, dense_errors::Bool)
   @ SciMLBase ~/Documents/git_repositories/julia/SciMLBase/src/solutions/rode_solutions.jl:132
 [5] calculate_solution_errors!(sol::RODESolution{Float64, 1, Vector{Float64}, Vector{Float64}, Dict{Symbol, Float64}, Vector{Float64}, NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Nothing}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Nothing}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, Nothing, SDEFunction{false, typeof(f), typeof(g), LinearAlgebra.UniformScaling{Bool}, typeof(f_analytic), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, typeof(g), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Nothing}, EM{true}, StochasticDiffEq.LinearInterpolationData{Vector{Float64}, Vector{Float64}}, DiffEqBase.DEStats})
   @ SciMLBase ~/Documents/git_repositories/julia/SciMLBase/src/solutions/rode_solutions.jl:110
 [6] top-level scope
   @ REPL[21]:1

julia> length(sol.u), length(sol.u_analytic)
(17, 34)

julia> SciMLBase.calculate_solution_errors!(sol);
ERROR: DimensionMismatch("dimensions must match: a has dims (Base.OneTo(17),), b has dims (Base.OneTo(51),), mismatch at 1")
Stacktrace:
 [1] promote_shape
   @ ./indices.jl:178 [inlined]
 [2] promote_shape
   @ ./indices.jl:169 [inlined]
 [3] -(A::Vector{Float64}, B::Vector{Float64})
   @ Base ./arraymath.jl:38
 [4] calculate_solution_errors!(sol::RODESolution{Float64, 1, Vector{Float64}, Vector{Float64}, Dict{Symbol, Float64}, Vector{Float64}, NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Nothing}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Nothing}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, Nothing, SDEFunction{false, typeof(f), typeof(g), LinearAlgebra.UniformScaling{Bool}, typeof(f_analytic), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, typeof(g), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Nothing}, EM{true}, StochasticDiffEq.LinearInterpolationData{Vector{Float64}, Vector{Float64}}, DiffEqBase.DEStats}; fill_uanalytic::Bool, timeseries_errors::Bool, dense_errors::Bool)
   @ SciMLBase ~/Documents/git_repositories/julia/SciMLBase/src/solutions/rode_solutions.jl:132
 [5] calculate_solution_errors!(sol::RODESolution{Float64, 1, Vector{Float64}, Vector{Float64}, Dict{Symbol, Float64}, Vector{Float64}, NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Nothing}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Float64, Nothing}, false}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, Nothing, SDEFunction{false, typeof(f), typeof(g), LinearAlgebra.UniformScaling{Bool}, typeof(f_analytic), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, typeof(g), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Nothing}, EM{true}, StochasticDiffEq.LinearInterpolationData{Vector{Float64}, Vector{Float64}}, DiffEqBase.DEStats})
   @ SciMLBase ~/Documents/git_repositories/julia/SciMLBase/src/solutions/rode_solutions.jl:110
 [6] top-level scope
   @ REPL[23]:1

julia> length(sol.u), length(sol.u_analytic)
(17, 51)

Since calculate_solution_errors!() is not exported and is always called after build_solution, I am not sure this is considered a bug, but still we might want to empty!(sol.u_analytic) at https://github.com/SciML/SciMLBase.jl/blob/master/src/solutions/rode_solutions.jl#L116 (and similarly for the other types of problem since this is not exclusive to SDEProblems), for the sake of generality.

Register

@JuliaRegistrator register()

Failed to show value: type OptimizationSolution has no field t

In Pluto, When I want to see the result of GalacticOptim optimization, I get this error:

Failed to show value:

type OptimizationSolution has no field t

getproperty@optimization_solutions.jl:39[inlined]
rows(::SciMLBase.OptimizationSolution{Float32, 1, Vector{Float32}, SciMLBase.OptimizationProblem{true, SciMLBase.OptimizationFunction{true, GalacticOptim.AutoZygote, ICNF.var"#f#12"{Float32, ICNF.RNODE{Float32}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float32}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Flux.Optimise.AMSGrad, Float32, Nothing})@tabletraits.jl:26
table_data(::SciMLBase.OptimizationSolution{Float32, 1, Vector{Float32}, SciMLBase.OptimizationProblem{true, SciMLBase.OptimizationFunction{true, GalacticOptim.AutoZygote, ICNF.var"#f#12"{Float32, ICNF.RNODE{Float32}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float32}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Flux.Optimise.AMSGrad, Float32, Nothing}, ::IOContext{IOBuffer})@PlutoRunner.jl:1414
show_richest(::IOContext{IOBuffer}, ::Any)@PlutoRunner.jl:1062
var"#sprint_withreturned#60"(::IOContext{Base.DevNull}, ::Int64, ::typeof(Main.PlutoRunner.sprint_withreturned), ::Function, ::SciMLBase.OptimizationSolution{Float32, 1, Vector{Float32}, SciMLBase.OptimizationProblem{true, SciMLBase.OptimizationFunction{true, GalacticOptim.AutoZygote, ICNF.var"#f#12"{Float32, ICNF.RNODE{Float32}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float32}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Flux.Optimise.AMSGrad, Float32, Nothing})@PlutoRunner.jl:1012
format_output_default(::Any, ::Any)@PlutoRunner.jl:920
var"#format_output#48"(::IOContext{Base.DevNull}, ::typeof(Main.PlutoRunner.format_output), ::Any)@PlutoRunner.jl:937
formatted_result_of(::Base.UUID, ::Bool, ::Vector{String}, ::Nothing, ::Module)@PlutoRunner.jl:841
top-level scope@none:1

bizarre heisenbug performance problem when solving `EnsembleProblem`

I'm getting some extremely bizarre performance characteristics when solving an EnsembleProblem with EnsembleThreads. In my particular case, I had already implemented a parallel solve with parallel threads using ThreadsX that solved in about 25 s for the examples I was trying. Expecting EnsembleProblem to be about the same or slightly faster, I was surprised to find it not only take a really long time but actually run out of memory (on 16 GB of RAM) before it ever got there.

I've done a ton of experimenting trying to isolate this but it has been extremely difficult. I have tried reproducing just about all aspects of my problem to create a MWE, but my attempted MWE was never slow, even at the point where I was solving the same equations (you can see this attempt here, again this was NOT slow). So, I went the other way, and started trying to simplify the context of my original problem.

I believe I now have simplified it about as much as possible, see my repo here. Running my solve method is incredibly slow and inconsistent, anywhere from 10 to 20 s for a mere 100 trajectories. At this point, I am only solving a simple 4-dimensional harmonic oscillator and it is almost as slow as solving for the geodesics on a wormhole manifold.

The only thing I can conclude at this point is that whatever is happening is only happening when I run this inside my package, but not if I just include a file with it.

Something else rather bizarre is that if I attempt to profile this using @profile the problem is partially alleviated (though I think still a bit slow):

julia> @time solve(C);
 21.132324 seconds (4.65 k allocations: 756.469 KiB)

julia> @profile @time solve(C);
  0.017434 seconds (4.66 k allocations: 756.516 KiB)

Because of this I have gotten exactly nowhere trying to track this down via profiling.

This is occurring on both Julia 1.6.2 and 1.7-beta2.

Add deprecation to `.minimizer` of OptimizationProblem

Scheduled for early 2022.

save_idxs in combination with PlotRecipe prints wrong symbols

When saving k indices, the first k symbols are used for the plot, instead of the symbols at the correct index.

using OrdinaryDiffEq, Plots
function func!(du,u,p,t)
  du[1] = u[1]
  du[2] = u[2]
  du[3] = - u[3]
end
odefunc! = ODEFunction(func!, syms = [:a,:b,:c])
odeprob  = ODEProblem(odefunc!,[1.0;1.0; 1.0],(0.0,1.0),nothing)
odesol   = solve(odeprob, Tsit5(), save_idxs = [1,3])
plot(odesol)

Here the idxs [1,3] are saved, but the plot shows the symbols [:a, :b].

Errors in EnsembleAnalysis

Taking my usual tired SIR example; this time as a jump system solved by SSA, and run for 1000 replicates in an ensemble:

using ModelingToolkit
using DiffEqJump
using DifferentialEquations

@parameters t β c γ
@variables S(t) I(t) R(t)

rxs = [Reaction(β*c/(S+I+R), [S,I], [I], [1,1], [2])
       Reaction(γ, [I], [R])]
rs  = ReactionSystem(rxs, t, [S,I,R], [β,c,γ])

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
ts = 0:δt:tmax
u0 = [S => 990, I => 10, R => 0]
p = [β =>0.05, c => 10.0, γ => 0.25]

jumpsys = convert(JumpSystem, rs)
dprob = DiscreteProblem(jumpsys, u0, tspan, p)
jprob = JumpProblem(jumpsys, dprob, Direct())
ensemble_jprob = EnsembleProblem(jprob)
ensemble_jsol = solve(ensemble_jprob,SSAStepper(),trajectories=1000)

The following works fine.

timeseries_point_meanvar(ensemble_jsol,ts)

However, the following gives an error:

timeseries_point_meancov(ensemble_jsol,ts)

ERROR: MethodError: no method matching timeseries_point_meancov(::EnsembleSolution{Int64,3,Array{ODESolution{Int64,2,Array{Array{Int64,1},1},Nothing,Nothing,Array{Float64,1},Nothing,DiscreteProblem{Array{Int64,1},Tuple{Float64,Float64},true,Array{Float64,1},DiscreteFunction{true,DiscreteFunction{true,SciMLBase.var"#140#141",Nothing,Nothing},Nothing,Array{Symbol,1}},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}},SSAStepper,SciMLBase.ConstantInterpolation{Array{Float64,1},Array{Array{Int64,1},1}},DiffEqBase.DEStats},1}}, ::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}})
Closest candidates are:
  timeseries_point_meancov(::Any, ::Any, ::Any) at /home/simon/.julia/packages/SciMLBase/Afx1r/src/ensemble/ensemble_analysis.jl:143
Stacktrace:
 [1] top-level scope at none:1

Incorrect x-axis limits for solution plot w/ function of states

Description
When plotting an ODE solution using a function f such as in vars = (f,2,3) the x-axis is set on the range of the state variable corresponding to the first function argument, i.e. u[vars[2]]. However, it should rather be set based on the range of the actually plotted values, which depend on the return value of f.

Actual Behavior

using Plots, OrdinaryDiffEq

f(x,p,t) = p*x
prob = ODEProblem(f,[1.0,2.0,3.0],(0.0,1.0),-.2)
sol = solve(prob, Tsit5())

f4(y,z) = (y,z) # f4 keeps order
f5(z,y) = (y,z) # f5 inverts order
f6(t,x,y,z) = (y,z) # f6 has unused arguments

# Working as expected with x-axis range of [1.5, 2.2]:
plot(sol, vars = (2,3))                # CORRECT: x-axis range [1.5, 2.2] equal to range of u[2]
plot(sol, vars = (f4,2,3))             # CORRECT: x-axis range [1.5, 2.2] equal to range of u[2]

# Not working as expected:
plot(sol, vars = (f5,3,2))             # WRONG: x-axis range [2.3, 3.3] equal to range of u[3]
plot(sol, vars = (f6,0,1,2,3))         # WRONG: x-axis range [0, 1]     equal to range of time and label "t"
plot(sol, vars = [(f4,2,3), (f5,3,2)]) # WRONG: x-axis range [1.5, 3.3] equal to overall range of u[2] and u[3]

Expected Behavior
Following statements should yield a plot with x-axis range [1.5, 2.2] and without axis label "t"

plot(sol, vars = (f5,3,2))             
plot(sol, vars = (f6,0,1,2,3))         
plot(sol, vars = [(f4,2,3), (f5,3,2)])

Related issues
OPEN: SciML/DiffEqBase.jl#591
CLOSED: SciML/DiffEqBase.jl#35

Possible solution
I pinpointed the behavior to the lines:

SciMLBase.jl/src/solutions/solution_interface.jl

Lines 204 to 212 in e44240f

 else 

 mins = minimum(sol[int_vars[1][2],:]) 

 maxs = maximum(sol[int_vars[1][2],:]) 

 for iv in int_vars 

 mins = min(mins,minimum(sol[iv[2],:])) 

 maxs = max(maxs,maximum(sol[iv[2],:])) 

 end 

 xlims --> ((1-sign(mins)*axis_safety)*mins,(1+sign(maxs)*axis_safety)*maxs) 

 end

Instead one could use plot_vecs, which acutally gets printed. Something along the lines below. However, I'm unsure whether this breaks other use cases.

mins = minimum(plot_vecs[1])
maxs = maximum(plot_vecs[1])

idxs vs. vars for accessing solution values

It would be nice if the same keyword could be used in plotting and in evaluating solutions for specific components, i.e. sol(t; vars=[...]) and plot(sol; vars=[...]) or vice-versa with idxs.

Compat chart link warning is missing

The link provided by the warning in src/utils.jl:100, namely http://docs.juliadiffeq.org/dev/basics/compatibility_chart.html, doesn't exist.

Add median to `EnsembleSummary`

Right now, EnsembleSummary contains the mean, variance, and quantiles. The corresponding plot recipe plots (by default) the mean, and the quantiles as ribbons. This can lead to unexpected plots for cases with strong outliers, where the line lies outside of the ribbon.

Proposal: Add the median to the ensemble summary. Then plot either the median (line) and the quantiles (ribbon) for ci_type=:quantile or the mean and variance as currently done when ci_type=:SEM.

Create an interface for PDE solutions

The idea is to create a PDESolution(prob,discretizer,res) that takes your problem, your discretizer, and your result. This interface can then have common actions like plot, get a certain variable at a certain interval, and so on.

It would also be great if the interface could make base manipulation to a solution similar to the TimeseriesSolution:

sol.y(0.1)
sol.approx(dx)
plot(sol)
sol.res

This would provide SciML with a standard way of analyzing PDEs – regardless of the method used to solve it.

First batch not counted in simulation time when solving EnsembleProblems.

Hi,

I found a minor issue when solving EnsembleProblems.

Due to the way the @elapsed are placed in this function :

 __solve(prob::AbstractEnsembleProblem,  alg::Union{DEAlgorithm,Nothing}, ensemblealg::BasicEnsembleAlgorithm; kwargs...)

if there are mutiple bathes or prob.reduction !== DEFAULT_REDUCTION then the first batch does not count in the simulation time (see below).

SciMLBase.jl/src/ensemble/basic_ensemble_solve.jl

Lines 55 to 82 in 9788d32

 if num_batches == 1 && prob.reduction === DEFAULT_REDUCTION 

 elapsed_time = @elapsed u = solve_batch(prob,alg,ensemblealg,1:trajectories,pmap_batch_size;kwargs...) 

 _u = tighten_container_eltype(u) 

 return EnsembleSolution(_u,elapsed_time,true) 

 end 

 converged::Bool = false 

 i = 1 

 II = (batch_size*(i-1)+1):batch_size*i 

 batch_data = solve_batch(prob,alg,ensemblealg,II,pmap_batch_size;kwargs...) 

 u = prob.u_init === nothing ? similar(batch_data, 0) : prob.u_init 

 u,converged = prob.reduction(u,batch_data,II) 

 elapsed_time = @elapsed for i in 2:num_batches 

 converged && break 

 if i == num_batches 

 II = (batch_size*(i-1)+1):trajectories 

 else 

 II = (batch_size*(i-1)+1):batch_size*i 

 end 

 batch_data = solve_batch(prob,alg,ensemblealg,II,pmap_batch_size;kwargs...) 

 u,converged = prob.reduction(u,batch_data,II) 

 end 

 _u = tighten_container_eltype(u) 

 return EnsembleSolution(_u,elapsed_time,converged)

This is especially bad when reduction is needed but there is only one batch (e.g. when trying to avoid a race condition with EnsembleThreads() ) or if the number of batches is low.

I believe an easy fix would be moving the second @elapsed before the first batch using a begin ... end block :

https://github.com/brainsMAKER/SciMLBase.jl/blob/7d1494e118b427c99f9eb7b611698f4d216709c5/src/ensemble/basic_ensemble_solve.jl#L61-L81

Depend on Symbolics

In order to replace issymbollike & solve method ambiguities for AD rules, SciMLBase could depend on Symbolics. Cf. this remark. by @sharanry in #85.

Currently, we can't do this b/c Symbolics.jl depends on SciMLBase here.

issymbollike from SciMLBase is extended in Symbolics: here

The original definition of issymbollike is in RecursiveArrayTools here. SciMLBase depends on RecursiveArrayTools here.

Symbolics could define an AbstractSymbol type, RecursiveArrayTools could use this for dispatching instead of branching on issymbollike, but in order to make this work we need to make sure Symbolics doesn't depend on SciMLBase.

Edit: It appears that the only mention of SciMLBase within Symbolics is where issymbollike is extended. I'm gonna do a smoke test & make sure.

kwargs field makes `IntegralProblem` type unstable

function copy_(prob)
    prob_ = IntegralProblem(prob.f, prob.lb, prob.ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...) 
end
@code_warntype copy_(prob)

SciML/Integrals.jl#65

CC: @ChrisRackauckas @YingboMa

Updating SciML to ChainRules

The list of things to handle is:

The solve.jl one is a bit nasty because it calls all of the adjoints in DiffEqSensitivity.jl, so it's somewhat breaking, so we can... just do it really quickly.

@frankschae @YingboMa

Stackoverflow with unexpected argument type to solution indexing

When I accidentally tried to index an ODESolution using a character instead of a symbol, I got a stackoverflow.

using ModelingToolkit, OrdinaryDiffEq
@parameters α β δ γ
@variables t x(t) y(t) dx(t)
D = Differential(t)
eqs = [
  dx ~ α*x - β*x*y, 
  D(x) ~ dx,
  D(y) ~ -δ*y + γ*x*y
]
@named lv = ODESystem(eqs)
syss = structural_simplify(lv) 
parms = [α => 1.5, β => 1.0, δ => 3.0, γ => 1.0]
x0 = [x => 1.0, y => 1.0]
tsteps = 0.0:0.1:10.0
prob = ODEProblem(syss, x0, extrema(tsteps), parms, jac = true)
soltrue = solve(prob,  Tsit5(), saveat = tsteps);
soltrue['a'] #Stackoverflow

The culprit is solution_interface.jl#60, which causes infinite recursion if you hit it with something that is neither an integer nor a symbol.

`EnsembleSolution` subset using `Num`s

sim[x1], where x1<: Num gives an ERROR: ArgumentError: invalid index: x1(t) of type Num currently (SciMLBase v1.49.1).

Using [s[x1] for s in sim] gives the desired output, but it will not be compatible with downstream analysis tools, like EnsembleSummary. I tried using

function ensemble_subset(sim, states)
    u = [DiffEqArray(s[states], s.t) for s in sim]
    EnsembleSolution(u, sim.elapsedTime, sim.converged)
end

and this works for Nums, but breaks plotting for Vector{Num} because of the plot recipe assumes that we have the shape of an ODE solution which needs transposing for plotting
https://github.com/SciML/RecursiveArrayTools.jl/blob/47bdb6eb3ca9af8def40e6806dc03a0f902a4209/src/vector_of_array.jl#L301-L308
due to the N+1 in
https://github.com/SciML/RecursiveArrayTools.jl/blob/47bdb6eb3ca9af8def40e6806dc03a0f902a4209/src/vector_of_array.jl#L97

This way the DiffEqArray{T,1} dispatch that was added in SciML/RecursiveArrayTools.jl#43 is not used and we get the generic fallback.

I've tried an other workaround by remaking the ODESolution

remake_sol(sol, u) = SciMLBase.build_solution(sol.prob, sol.alg, sol.t, u)

but this breaks because a lot of functions will assume that all states are present and I'm not sure it's a good idea.

This issue seems related to #87

Rewrite default plotting system

As noted in https://github.com/SciML/SciMLBase.jl/pull/9/files#issuecomment-770289578, the default plotting system is somewhat of a mess because it's just faking functionality that we now have, like:

Indexing via symbols
Observed variables
Independent variable symbol handling

It should get a refactor to greatly decrease the complexity by using these aspects of the system.

Figures

It figures. Made using LucidChart.

`set_u!` with symbols

Currently modifying the state of an integrator requires passing a Vector to set_u!, which is nontrivial when dealing with larger systems. Is there a way to do this symbolically, that I can't find? E.g. set_u!(integrator, x, 1.23)? Would this be feasible to implement?

Finish the OptimizationProblem docstring

Left incomplete at the end of #157 . @Vaibhavdixit02 could I enlist you to help this get finished?

store pointer to original `DiffEqArrayOperator` (if not overwritten) in `FactorizedDiffEqArrayOperator`

so we don't have to convert to AbstractMatrix every time we call *, mul!

SciMLBase.jl/src/operators/basic_operators.jl

Lines 121 to 127 in bfc028d

 struct FactorizedDiffEqArrayOperator{T<:Number, FType <: Union{ 

 Factorization{T}, Diagonal{T}, Bidiagonal{T}, 

 Adjoint{T,<:Factorization{T}}, 

 } 

 } <: AbstractDiffEqLinearOperator{T} 

 F::FType 

 end

SciMLBase.jl/src/operators/common_defaults.jl

Lines 13 to 18 in bfc028d

 for op in (:*, :/, :\) 

 @eval Base.$op(L::AbstractDiffEqLinearOperator, x::AbstractArray) = $op(convert(AbstractMatrix,L), x) 

 @eval Base.$op(x::AbstractArray, L::AbstractDiffEqLinearOperator) = $op(x, convert(AbstractMatrix,L)) 

 @eval Base.$op(L::DiffEqArrayOperator, x::Number) = $op(convert(AbstractMatrix,L), x) 

 @eval Base.$op(x::Number, L::DiffEqArrayOperator) = $op(x, convert(AbstractMatrix,L)) 

 end

LabelledArrays support for Statistics Summary

Hi !

Currently all SciMLBase.EnsembleAnalysis functions output LabelledArrays if the solution contains LabelledArrays with two exceptions: median and quantile (both timestep and timepoint versions). The MWE is:

using SciMLBase.EnsembleAnalysis, OrdinaryDiffEq, LabelledArrays

prob = ODEProblem((u,p,t)->1.01u,LVector(u1=0.5),(0.0,1.0))

function prob_func(prob,i,repeat)
  remake(prob,u0=rand()*prob.u0)
end

ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
sim = solve(ensemble_prob,Tsit5(),EnsembleDistributed(),trajectories=10)

The result of timestep_mean is a LabelledArray but timestep_median timestep_quantile produce just Vector{Float64}

typeof(timestep_mean(sim,2))
# LArray{Float64, 1, Vector{Float64}, (:u1,)}

typeof(timestep_median(sim,2))
# Vector{Float64} (alias for Array{Float64, 1})

The problem is both in vecarr_to_vectors functions, which converts LVector to Vector and loop, which computes median, quantile.

SciMLBase.jl/src/ensemble/ensemble_analysis.jl

Line 11 in 67d805e

return vecarr_to_vectors(VectorOfArray(arr))

https://github.com/SciML/RecursiveArrayTools.jl/blob/deeb0f69ee97b283a423a15ecc1c34721ec695cc/src/vector_of_array.jl#L210

I have added methods to fix it in my code but it seems a larger question of RecursiveArrayTools compatibility with LabelledArrays

Many functions on `AbstractODESolution` return actual `ODESolution`s

https://github.com/SciML/SciMLBase.jl/blob/master/src/solutions/ode_solutions.jl defines a range of functions, which operate on AbstractODESolutions, and yet return ODESolutions. This is particularly confusing for cases such as solution_new_retcode or solution_new_tslocation. For example, if I implement my own MySolution <: AbstractODESolution then it silently gets transformed into a ODESolution. I would argue that defining these functions on ODESolution would therefore be more appropriate, so that an error is raised.

The relevant functions here are:

build_solution
solution_new_retcode
solution_new_tslocation
solution_slice
sensitivity_solution

Documentation of the fields in the problem types

#151 added a bunch of documentation for the problem types, but took some shortcuts when doing the fields. The only one that did it properly was the ODEProblem:

SciMLBase.jl/src/problems/ode_problems.jl

Lines 69 to 82 in 80d0756

 struct ODEProblem{uType,tType,isinplace,P,F,K,PT} <: 

 AbstractODEProblem{uType,tType,isinplace} 

 """The ODE is `du = f(u,p,t)` for out-of-place and f(du,u,p,t) for in-place.""" 

 f::F 

 """The initial condition is `u(tspan[1]) = u0`.""" 

 u0::uType 

 """The solution `u(t)` will be computed for `tspan[1] ≤ t ≤ tspan[2]`.""" 

 tspan::tType 

 """Constant parameters to be supplied as the second argument of `f`.""" 

 p::P 

 """A callback to be applied to every solver which uses the problem.""" 

 kwargs::K 

 problem_type::PT 

 """An internal argument for storing traits about the solving process."""

All of the others didn't do it via docstrings on the fields. We should get that fixed, it's a nice first issue since it's quite cookie-cutter fix!

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Defining an `ODEProblem` with kwarg `cb` doesn't apply callback or throw error

I hope this is the right place to open this issue.

There is a slight mismatch in the interfaces of ODEProblem and DiffEqFlux's sciml_train. The former allows for callbacks through the kwarg callback while the latter uses cb.
I fell into the "trap" of using cb on an ODEProblem and spent quite a while wondering why my callbacks weren't triggered.

Maybe cb could be made into an alias for callback?

DAEProblem without differential initial conditions

It would be nice to have a DAEProblem constructor that did not require the user to specify du0, as this is often redundant.
Or are there good reasons to require this?

Overload constructorof on SciMLProblems

It would be somewhat like #54 but note that the type parameters of SciMLProblems don't have the same order as their fields. One needs to be careful with the overloading.

Downstream interface changes

NonlinearSolve should do a few wrappers to NonlinearProblem/NonlinearSolve
SteadyStateDiffEq should use NonlinearProblem/NonlinearSolve in SSRootFind

ERROR: LoadError: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 652 and 139")

I am trying to run an UDE of a Rosenzweig-MacArthur model with noisy data. I get the following Warning and Error when I try to run it.
┌ Warning: Interrupted. Larger maxiters is needed. └ @ SciMLBase ~/.julia/packages/SciMLBase/UIp7W/src/integrator_interface.jl:331 ERROR: LoadError: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 652 and 139")

Here's the data for further checking
RM_data.csv

And the relevant code

X = CSV.read("Data/RM_data.csv",DataFrame)
println(first(X))
ts = X.timestamp
Xs = X.value1
Ys = X.value2
Zs = X.value3

##Generate subsamples of data
TS_length = length(Xs)
println(TS_length)
N_samples = 5
subsampleLength = floor(Int64,TS_length/N_samples)

#Subsample the data
#We assume we have only 5 measurements, evenly distributed
TS = zeros(N_samples,subsampleLength)
XS = zeros(N_samples,subsampleLength)
YS = zeros(N_samples,subsampleLength)
ZS = zeros(N_samples,subsampleLength)

for i in 1:N_samples
    println(i)
    XS[i, :] = Xs[((i-1)*subsampleLength+1):(i*subsampleLength)]
    YS[i, :] = Ys[((i-1)*subsampleLength+1):(i*subsampleLength)]
    ZS[i, :] = Zs[((i-1)*subsampleLength+1):(i*subsampleLength)]
    TS[i, :] = ts[((i-1)*subsampleLength+1):(i*subsampleLength)]
end

#Define nonlinear network
# Gaussian RBF as activation
rbf(x) = exp.(-(x.^2))

# Define the network 3->5->5->3
U = FastChain(
    FastDense(3,5,rbf), FastDense(5,5, rbf), FastDense(5,3)
)

#Define the initial parameters
p = vcat(rand(Float32,8), initial_params(U))
#p = rand(Float32,8)

#Define hybrid model
function ude_dynamics!(du,u, p, t, p_true)
    r = p[1]
    K = p[2]
    a2 = p[3]
    c2 = p[4]
    d2 = p[5]
    a3 = p[6]
    d3 = p[7]
    c3 = p[8]

    û = U(u, p[9:end]) # Network prediction
    #Define the known dynamics
    du[1] = r*u[1]*(1-u[1]/K) - a2*u[1]*u[2]/(1+u[1]) + û[1]
    du[2] = c2*a2*u[1]*u[2]/(1+u[1]) - d2*u[2] - a3*u[2]*u[3]/(1+u[2]) + û[2]
    du[3] = c3*a3*u[2]*u[3]/(1+u[2]) - d3*u[3] + û[3]
end

# Closure with the known parameter
nn_dynamics!(du,u,p,t) = ude_dynamics!(du,u,p,t,nothing)
# Define the problem
prob_nn = ODEProblem(nn_dynamics!,[XS[1,1],YS[1,1],ZS[1,1]], (ts[1],ts[end]), p)

## Function to train the network
# Define a predictor
function predict(θ, X = X[:,1], T = ts)
    Array(solve(prob_nn, Vern7(), u0 = X, p=θ,
                tspan = (T[1], T[end]), saveat = T,
                abstol=1e-6, reltol=1e-6,
                sensealg = ForwardDiffSensitivity()
                ))
end

# Multiple shooting like loss
function loss(θ)
    # Start with a regularization on the network
    l = convert(eltype(θ), 1e-3)*sum(abs2, θ) ./ length(θ)
    for i in 1:size(XS,1)
        XX= predict(θ, [XS[i,1], YS[i,1], ZS[i,1]], TS[i, :])
        # Full prediction in x
        l += sum(abs2, XS[i,:] .- XX[1,:])
        l += sum(abs2, YS[i, :] .- XX[2, :])
        l += sum(abs2, ZS[i, :] .- XX[3, :])
    end
    return l
end

# Container to track the losses
losses = Float32[]

# Callback to show the loss during training
callback(θ,l) = begin
    push!(losses, l)
    if length(losses)%50==0
        println("Current loss after $(length(losses)) iterations: $(losses[end])")
    end
    false
end

## Training

# First train with ADAM for better convergence -> move the parameters into a
# favourable starting positing for BFGS
res1  = DiffEqFlux.sciml_train(loss, p, ADAM(0.1f0), cb=callback, maxiters = 10000)
# Train with BFGS
res2 = DiffEqFlux.sciml_train(loss, res1.minimizer, BFGS(initial_stepnorm=0.01f0), cb=callback, maxiters = 10000, g_tol = 1e-10)
println("Final training loss after $(length(losses)) iterations: $(losses[end])")

# Plot the losses
pl_losses = plot(1:10000, losses[1:10000], yaxis = :log10, xaxis = :log10, xlabel = "Iterations", ylabel = "Loss", label = "ADAM", color = :blue)
plot!(10001:length(losses), losses[10001:end], yaxis = :log10, xaxis = :log10, xlabel = "Iterations", ylabel = "Loss", label = "BFGS", color = :red)
#savefig(pl_losses, joinpath(pwd(), "plots", "$(svname)_losses.pdf"))

# Rename the best candidate
p_trained = res2.minimizer

I don't get any print of the callback, which makes me assume that it fails on the first iteration (especially since increasing the number of iterations doesn't change the result). I tried running the equivalent NODE (commenting any reference to û) and got prints of the callback and the process finishes, which tells me the issue is being caused by the û term.

Remove getindepsym_defaultt

The only reason it still exists is for backwards compatbility, i.e. plot(sol,vars=(:t,:x) where independent variables is not set. After #9 it should be done by matching the indepvar, in which case :t doesn't need to be a hardcoded default and should get deleted. This should be scheduled for about a year after DSLs like MTK and ParameterizedFunctions.jl and everything else stabilizes and starts appropriately defining the independent variables.

indexing NonlinearSolution with vector of Num errors

using ModelingToolkit, OrdinaryDiffEq, SteadyStateDiffEq

@parameters t
@variables x(t) y(t)
D = Differential(t)

eqs = [D(x) ~ -x,
    D(y) ~ -y
]

@named sys = ODESystem(eqs)
prob = SteadyStateProblem(sys, [1.0, 1.0])
sol = solve(prob, DynamicSS(Tsit5()))
@test_throws StackOverflowError sol[[x, y]]

versions

  [961ee093] ModelingToolkit v8.14.1
  [1dea7af3] OrdinaryDiffEq v6.17.0
  [9672c7b4] SteadyStateDiffEq v1.8.0

Finish the QuadratureProblem docstring

Left unfinished at the end of #157 . @utkarsh530 could I enlist you to help get this one finished?

	else
	mins = minimum(sol[int_vars[1][2],:])
	maxs = maximum(sol[int_vars[1][2],:])
	for iv in int_vars
	mins = min(mins,minimum(sol[iv[2],:]))
	maxs = max(maxs,maximum(sol[iv[2],:]))
	end
	xlims --> ((1-sign(mins)axis_safety)mins,(1+sign(maxs)axis_safety)maxs)
	end

	if num_batches == 1 && prob.reduction === DEFAULT_REDUCTION
	elapsed_time = @elapsed u = solve_batch(prob,alg,ensemblealg,1:trajectories,pmap_batch_size;kwargs...)
	_u = tighten_container_eltype(u)
	return EnsembleSolution(_u,elapsed_time,true)
	end

	converged::Bool = false
	i = 1
	II = (batch_size(i-1)+1):batch_sizei

	batch_data = solve_batch(prob,alg,ensemblealg,II,pmap_batch_size;kwargs...)

	u = prob.u_init === nothing ? similar(batch_data, 0) : prob.u_init
	u,converged = prob.reduction(u,batch_data,II)
	elapsed_time = @elapsed for i in 2:num_batches
	converged && break
	if i == num_batches
	II = (batch_size*(i-1)+1):trajectories
	else
	II = (batch_size(i-1)+1):batch_sizei
	end
	batch_data = solve_batch(prob,alg,ensemblealg,II,pmap_batch_size;kwargs...)
	u,converged = prob.reduction(u,batch_data,II)
	end

	_u = tighten_container_eltype(u)

	return EnsembleSolution(_u,elapsed_time,converged)

	struct FactorizedDiffEqArrayOperator{T<:Number, FType <: Union{
	Factorization{T}, Diagonal{T}, Bidiagonal{T},
	Adjoint{T,<:Factorization{T}},
	}
	} <: AbstractDiffEqLinearOperator{T}
	F::FType
	end

	for op in (:*, :/, :\)
	@eval Base.$op(L::AbstractDiffEqLinearOperator, x::AbstractArray) = $op(convert(AbstractMatrix,L), x)
	@eval Base.$op(x::AbstractArray, L::AbstractDiffEqLinearOperator) = $op(x, convert(AbstractMatrix,L))
	@eval Base.$op(L::DiffEqArrayOperator, x::Number) = $op(convert(AbstractMatrix,L), x)
	@eval Base.$op(x::Number, L::DiffEqArrayOperator) = $op(x, convert(AbstractMatrix,L))
	end

	struct ODEProblem{uType,tType,isinplace,P,F,K,PT} <:
	AbstractODEProblem{uType,tType,isinplace}
	"""The ODE is `du = f(u,p,t)` for out-of-place and f(du,u,p,t) for in-place."""
	f::F
	"""The initial condition is `u(tspan[1]) = u0`."""
	u0::uType
	"""The solution `u(t)` will be computed for `tspan[1] ≤ t ≤ tspan[2]`."""
	tspan::tType
	"""Constant parameters to be supplied as the second argument of `f`."""
	p::P
	"""A callback to be applied to every solver which uses the problem."""
	kwargs::K
	problem_type::PT
	"""An internal argument for storing traits about the solving process."""