Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Activity
  • 08:46
    robertfeldt opened #143
  • 08:15
    robertfeldt closed #142
  • 05:45
    amadeusine starred SciML/RecursiveArrayTools.jl
  • 00:32
    ChrisRackauckas commented #345
  • Apr 01 21:48
  • Apr 01 20:47
    tkf commented #241
  • Apr 01 19:25
    ViralBShah closed #238
  • Apr 01 19:25
    ViralBShah commented #238
  • Apr 01 19:24
    ViralBShah closed #227
  • Apr 01 19:24
    ViralBShah commented #227
  • Apr 01 17:55
    aarmey opened #233
  • Apr 01 17:46
    FHell commented #583
  • Apr 01 17:45
    FHell commented #583
  • Apr 01 16:38
    ChrisRackauckas commented #241
  • Apr 01 16:33

    ChrisRackauckas on master

    Fix reduction func in ensembles… Merge pull request #345 from jl… (compare)

  • Apr 01 16:33
    ChrisRackauckas closed #345
  • Apr 01 16:32
    ChrisRackauckas commented #583
  • Apr 01 16:12
    omalleyc edited #584
  • Apr 01 16:11
    omalleyc opened #584
  • Apr 01 15:53

    ChrisRackauckas on master

    Small fixes to the DAE example … Merge pull request #346 from as… (compare)

BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> oh yes, something about the KLU binding seems off in the latest Julia versions
[slack] <rveltz> :thumbsup:
[slack] <rveltz> I am always worried my question is dumb
BridgingBot
@GitterIRCbot

[slack] <Thomas Propson> Hi all, I am new to both Julia and the DiffEq package, but I have experience with differential equations and autodiff. I am attempting to do an adjoint sensitivity analysis on a cost function on the final state of my system. Following the documentation [0] it looks like I want to use the Zygote package and the concrete_solve interface. It looks like I could do something similar with the adjoint_sensitivities interface, but I can’t express the derivative of my cost function easily and want to backprop rather than forwardprop or finite difference. I am not sure what I am doing wrong. This is not the problem that I actually want to solve but this is the minimum example I put together to show what I am trying to do.
```# foo.jl
using DiffEqSensitivity
using DifferentialEquations
using Zygote

module Foo

function rhs(u, p, t)
return u * p
end

function cost(u, p, prob)
sol = concrete_solve(prob, Tsit5(), u, p,
sensealg=QuadratureAdjoint())
final_u = sol.u[length(sol.u)]
return final_u * 2
end

function main()
u0 = 1.0
p = 1.0
end_time = 1.0
tspan = (0.0, end_time)
prob = ODEProblem(rhs, u0, tspan, p=p)
res = Zygote.gradient((u, p) -> cost(u, p, prob), u0, p)
println(res)
end

endThe full stack trace is:julia> include("foo.jl")
WARNING: replacing module Foo.
Main.Foo

julia> Main.Foo.main()
ERROR: BoundsError: attempt to access 0-element Array{Any,1} at index []
Stacktrace:
[1] throw_boundserror(::Array{Any,1}, ::Tuple{}) at ./abstractarray.jl:537
[2] checkbounds at ./abstractarray.jl:502 [inlined]
[3] _getindex at ./abstractarray.jl:1002 [inlined]
[4] getindex(::Array{Any,1}) at ./abstractarray.jl:980
[5] (::Zygote.var"#back#183"{:u,Zygote.Context,RecursiveArrayTools.DiffEqArray{Float64,1,Array{Float64,1},Array{Float64,1}},Array{Float64,1}})(::Array{Float64,1}) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/lib/lib.jl:199
[6] #355#back at /home/tcpropson/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:49 [inlined]
[7] cost at /home/tcpropson/repos/qoc/_test/foo.jl:89 [inlined]
[8] (::typeof(∂(cost)))(::Float64) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface2.jl:0
[9] #1 at /home/tcpropson/repos/qoc/_test/foo.jl:99 [inlined]
[10] (::typeof(∂(λ)))(::Float64) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface2.jl:0
[11] (::Zygote.var"#36#37"{typeof(∂(λ))})(::Float64) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface.jl:36
[12] gradient(::Function, ::Float64, ::Vararg{Float64,N} where N) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface.jl:45
[13] main() at /home/tcpropson/repos/qoc/_test/foo.jl:99
[14] top-level scope at REPL[94]:1```
[0] https://docs.sciml.ai/dev/analysis/sensitivity/#concrete_solve-Examples-1

[slack] <chrisrackauckas> See if `
Array(concrete_solve(prob, Tsit5(), u, p, sensealg=QuadratureAdjoint()))
fixes it
[slack] <Thomas Propson> No, unfortunately. Thank you for your quick response. The new error is:
julia> Main.Foo.main() ERROR: MethodError: no method matching vec(::Float64) Closest candidates are: vec(::LinearAlgebra.Transpose{T,#s662} where #s662<:(AbstractArray{T,1} where T) where T) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/adjtrans.jl:201 vec(::SparseArrays.AbstractSparseArray{Tv,Ti,1} where Ti where Tv) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/SparseArrays/src/sparsevector.jl:912 vec(::StaticArrays.StaticArray) at /home/tcpropson/.julia/packages/StaticArrays/1g9bq/src/abstractarray.jl:163 ... Stacktrace: [1] _broadcast_getindex_evalf at ./broadcast.jl:631 [inlined] [2] _broadcast_getindex at ./broadcast.jl:604 [inlined] [3] getindex at ./broadcast.jl:564 [inlined] [4] copy at ./broadcast.jl:854 [inlined] [5] materialize at ./broadcast.jl:820 [inlined] [6] Array(::RecursiveArrayTools.DiffEqArray{Float64,1,Array{Float64,1},Array{Float64,1}}) at /home/tcpropson/.julia/packages/RecursiveArrayTools/euiv5/src/vector_of_array.jl:13 [7] adjoint at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/lib/array.jl:7 [inlined] [8] _pullback(::Zygote.Context, ::Type{Array}, ::RecursiveArrayTools.DiffEqArray{Float64,1,Array{Float64,1},Array{Float64,1}}) at /home/tcpropson/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:47 [9] cost at /home/tcpropson/repos/qoc/_test/foo.jl:87 [inlined] [10] _pullback(::Zygote.Context, ::typeof(Main.Foo.cost), ::Float64, ::Float64, ::DiffEqBase.ODEProblem{Float64,Tuple{Float64,Float64},false,DiffEqBase.NullParameters,DiffEqBase.ODEFunction{false,typeof(Main.Foo.rhs),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:p,),Tuple{Float64}}},DiffEqBase.StandardODEProblem}) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface2.jl:0 [11] #1 at /home/tcpropson/repos/qoc/_test/foo.jl:99 [inlined] [12] _pullback(::Zygote.Context, ::Main.Foo.var"#1#2"{DiffEqBase.ODEProblem{Float64,Tuple{Float64,Float64},false,DiffEqBase.NullParameters,DiffEqBase.ODEFunction{false,typeof(Main.Foo.rhs),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:p,),Tuple{Float64}}},DiffEqBase.StandardODEProblem}}, ::Float64, ::Float64) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface2.jl:0 [13] _pullback(::Function, ::Float64, ::Float64) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface.jl:29 [14] pullback(::Function, ::Float64, ::Float64) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface.jl:35 [15] gradient(::Function, ::Float64, ::Vararg{Float64,N} where N) at /home/tcpropson/.julia/packages/Zygote/KNUTW/src/compiler/interface.jl:44 [16] main() at /home/tcpropson/repos/qoc/_test/foo.jl:99 [17] top-level scope at REPL[96]:1
[slack] <chrisrackauckas> oh I see what it is
[slack] <chrisrackauckas> the issue is that the adjoint methods don't work with non-mutable values right now, so you need to use an array instead of a scalar.
BridgingBot
@GitterIRCbot
[slack] <Thomas Propson> as in, make my state u an array?
[slack] <chrisrackauckas> yes
BridgingBot
@GitterIRCbot
[slack] <Thomas Propson> I think my problem is that I can’t grab the final state the way I’m currently trying. For instance, this works:
function cost(u0, p, prob) sol = concrete_solve(prob, Tsit5(), u0, p, saveat=0.1, sensealg=QuadratureAdjoint()) return sum(sol) end
But this doesn’t:
function cost(u0, p, prob) sol = concrete_solve(prob, Tsit5(), u0, p, saveat=0.1, sensealg=QuadratureAdjoint()) final_u = sol.u[length(sol.u)] return final_u[1] end
My question is: how do I get the gradients with respect to my cost function evaluated only at the final state.
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> interesting, it looks like that should work
pjentsch0
@pjentsch0
Hello! I'm trying to figure out why my DiffEq function is allocating. If I post the code here could someone take a look?
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> I don't have very much time, but could take a peak.
pjentsch0
@pjentsch0
The functions are just conditionals
function τ(t::Float64, p::parameters)::Float64
return t<p.t_drug ? 0 : p.ε_D p.τ_0
end
function ψ(t::Float64,p::parameters)::Float64
return t < p.t_vac ? 0 : p.ε_V
p.ψ_0
end
results of @btime:
1.960 μs (28 allocations: 1.56 KiB)
I know it's not allocating much but I feel like this model should be much faster.
@ChrisRackauckas your diffeq library is awesome and it makes my research with SDEs muuuuch easier, thank you!
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> that's the results of @btime directly on this function?
pjentsch0
@pjentsch0
yes
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> Maybe you want to use DiffEqBase.@.. du[n_g+1:n_g*2] = p.r * S*p.β_x_I - p.σ.*E
pjentsch0
@pjentsch0
@btime system(u_0,u_0,params,0.0) is what I use to time it
use that macro rather than the broadcasting?
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> that macro will do broadcasting, just possibly fix the allocations
[slack] <chrisrackauckas> for long expressions
[slack] <chrisrackauckas> I don't know if that's long enough, but worth a try
[slack] <chrisrackauckas> or maybe you might need to @views
pjentsch0
@pjentsch0
ok, I'll give it a try. Thanks!
BridgingBot
@GitterIRCbot

[slack] <rveltz> Hi,

Is this a bug?

julia> Dirichlet0BC(ComplexF64) ERROR: MethodError: no method matching RobinBC(::Tuple{Complex{Float64},Complex{Float64},Complex{Float64}}, ::Tuple{Complex{Float64},Complex{Float64},Complex{Float64}}, ::Float64, ::Float64)

BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> possibly
[slack] <chrisrackauckas> I think it may have too tight of a type on it
BridgingBot
@GitterIRCbot
[slack] <rveltz> thank you for your answer. Also, how can you impose multiple BC, can you compose them?
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> there's a multi-BC constructor
HelgavonLichtenstein
@HelgavonLichtenstein
Dear @kanav any luck with the Non bracketing interval passed in bisection method error? I wish I could help, but this is really out of my league.
C O Malley
@omalleyc
Hi, in the release of the latest version it mentions that the adjoint methods can now support singular mass matrices. Are there any examples for this?
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> @omalleyc
[slack] <chrisrackauckas> ``` function rober(du,u,p,t)
y₁,y₂,y₃ = u
k₁,k₂,k₃ = p
du[1] = -k₁y₁+k₃y₂y₃
du[2] = k₁
y₁-k₂y₂^2-k₃y₂*y₃
du[3] = y₁ + y₂ + y₃ - 1
nothing
end
M = [1. 0 0
0 1. 0
0 0 0]
f = ODEFunction(rober,mass_matrix=M)
p = [0.04,3e7,1e4]
prob_singular_mm = ODEProblem(f,[1.0,0.0,0.0],(0.0,100),p)
sol_singular_mm = solve(prob_singular_mm,Rodas5(),reltol=1e-8,abstol=1e-8)
ts = [50, sol_singular_mm.t[end]]
dg_singular(out,u,p,t,i) = (fill!(out, 0); out[end] = -1)
_, res = adjoint_sensitivities(sol_singular_mm,alg,dg_singular,ts,abstol=1e-5,reltol=1e-5,sensealg=QuadratureAdjoint())
reference_sol = ForwardDiff.gradient(p->G(p, prob_singular_mm, ts, sol->sum(last, sol.u)), vec(p))
@test res' ≈ reference_sol rtol = 1e-2

_, res_interp = adjoint_sensitivities(sol_singular_mm,alg,dg_singular,ts,abstol=1e-5,reltol=1e-5,sensealg=InterpolatingAdjoint())
@test res_interp ≈ res rtol = 1e-2
_, res_interp2 = adjoint_sensitivities(sol_singular_mm,alg,dg_singular,ts,abstol=1e-5,reltol=1e-5,sensealg=InterpolatingAdjoint(checkpointing=true),checkpoints=sol_singular_mm.t[1:10:end])
@test res_interp2 ≈ res rtol = 1e-2```
BridgingBot
@GitterIRCbot
[slack] <kanav> @HelgavonLichtenstein sorry I was busy in college projects, I ran your latest code on latest DifferentialEquations and I don’t see any error. Did you try updating DifferentialEquations? Or maybe just DiffEqBase
C O Malley
@omalleyc
@ChrisRackauckas , thanks for the example, but I am getting a type error for the adjoint_sensitivites function, "expected Float64, got ForwardDiff.Dual{Nothing,Float64,1}"
C O Malley
@omalleyc
@ChrisRackauckas , ah it seemed I needed to set autodiff=false in alg, then it works. Should the function ODELocalSensitivityProblem also work or just adjoint_sensitivity?
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> @omalleyc ODELocalSensitivityProblem with ForwardDiffSensitivity will work
[slack] <chrisrackauckas> The other forward mode choice, open an issue on it
[slack] <chrisrackauckas> it's easy to fix but I haven't done it yet.
C O Malley
@omalleyc
@ChrisRackauckas I get another error using ODELocalSensitivityProblem with ForwardDiffSensitivity that appears to be caused by the mass matrix: "ERROR: DimensionMismatch("W: (Base.OneTo(12), Base.OneTo(12)), mass matrix: (Base.OneTo(6), Base.OneTo(6))")". I will raise this as an issue. Thanks for your help
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> thanks yeah, I just hadn't thought of doing those together, but mathematically it's fairly simple so I could get that fixed up in no time, just not today.
C O Malley
@omalleyc
@ChrisRackauckas no problem, thanks again!
BridgingBot
@GitterIRCbot
[slack] <rveltz> thank you
BridgingBot
@GitterIRCbot

[slack] <Thomas Propson> Has anyone tried to use the adjoint sensitivities interface on complex inputs? It appears that the Tracker.Tracked objects don’t have operations defined for complex objects. Maybe one of the other vjp options [0] works, but I don’t see how I am supposed to set this. I’m happy to poke around but I am new to Julia and could use some advice about where to look.

[0] https://docs.sciml.ai/dev/analysis/sensitivity/#Internal-Automatic-Differentiation-Options-(ADKwargs)-1

BridgingBot
@GitterIRCbot

[slack] <Thomas Propson> I see this issue [1], but I’m wondering if it is outdated because libraries like Zygote now support complex derivatives [2]?

[1] SciML/DifferentialEquations.jl#462

[2] https://fluxml.ai/Zygote.jl/latest/complex/