-
Notifications
You must be signed in to change notification settings - Fork 1
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
ForwardDiff Decapode Gradient Fail Example #182
Comments
So Decapodes has two different code generation options, one that uses PreallocationTools and the DiffCaches, and one that doesn't. This is an example of code that does use PreallocationTools, and ForwardDiff doesn't work. |
@ChrisRackauckas I think the next steps for Decapode calibration stuff is to make sure that ForwardDiff works in situations like this, and then making sure that For context, I have a heat equation decapode example that uses PreallocationTools, ForwardDiff works, and |
@ChrisRackauckas I think this is the most pared down I can get it while still getting the error. using Catlab
using Catlab.Graphics
using CombinatorialSpaces
using Decapodes
using ComponentArrays
using ForwardDiff
using Zygote
using SciMLSensitivity
# External Dependencies
using MLStyle
using MultiScaleArrays
using LinearAlgebra
using OrdinaryDiffEq
using JLD2
using SparseArrays
using Statistics
using GeometryBasics: Point2, Point3
Point2D = Point2{Float64};
Point3D = Point3{Float64};
using DiagrammaticEquations
using DiagrammaticEquations.Deca
halfar_eq2 = @decapode begin
h::Form0
Γ::Form1
n::Constant
ḣ == ∂ₜ(h)
ḣ == ∘(⋆, d, ⋆)(Γ * d(h) * avg₀₁(mag(♯(d(h)))^(n - 1)) * avg₀₁(h^(n + 2)))
end
glens_law = @decapode begin
Γ::Form1
(A, ρ, g, n)::Constant
Γ == (2.0 / (n + 2.0)) * A * (ρ * g)^n
end
@info("Decapodes Defined")
ice_dynamics_composition_diagram = @relation () begin
dynamics(Γ, n)
stress(Γ, n)
end
ice_dynamics_cospan = oapply(ice_dynamics_composition_diagram,
[Open(halfar_eq2, [:Γ, :n]),
Open(glens_law, [:Γ, :n])])
ice_dynamics = apex(ice_dynamics_cospan)
ice_dynamics1D = expand_operators(ice_dynamics)
infer_types!(ice_dynamics1D, op1_inf_rules_1D, op2_inf_rules_1D)
resolve_overloads!(ice_dynamics1D, op1_res_rules_1D, op2_res_rules_1D)
s_prime = EmbeddedDeltaSet1D{Bool,Point2D}()
add_vertices!(s_prime, 10, point=Point2D.(range(-2.0, 2.0, length= 10), 0.0))
add_edges!(s_prime, 1:nv(s_prime)-1, 2:nv(s_prime))
orient!(s_prime)
s = EmbeddedDeltaDualComplex1D{Bool,Float64,Point2D}(s_prime)
subdivide_duals!(s, Circumcenter())
function generate(sd, my_symbol; hodge=GeometricHodge())
op = @match my_symbol begin
:♯ => x -> begin
# This is an implementation of the "sharp" operator from the exterior
# calculus, which takes co-vector fields to vector fields.
# This could be up-streamed to the CombinatorialSpaces.jl library. (i.e.
# this operation is not bespoke to this simulation.)
e_vecs = map(edges(sd)) do e
point(sd, sd[e, :∂v0]) - point(sd, sd[e, :∂v1])
end
neighbors = map(vertices(sd)) do v
union(incident(sd, v, :∂v0), incident(sd, v, :∂v1))
end
n_vecs = map(neighbors) do es
[e_vecs[e] for e in es]
end
map(neighbors, n_vecs) do es, nvs
sum([nv * norm(nv) * x[e] for (e, nv) in zip(es, nvs)]) / sum(norm.(nvs))
end
end
:mag => x -> norm.(x)
x => error("Unmatched operator $my_symbol")
end
return (args...) -> op(args...)
end
decapode_code = gensim(ice_dynamics1D, dimension=1, preallocate=true)
file = open("ice_sheet1D_alloc.jl", "w")
write(file, string("decapode_f = ", decapode_code))
close(file)
include("ice_sheet1D_alloc.jl")
fₘ = decapode_f(s, generate)
h₀ = map(x -> exp(-2 * x[1]^2), point(s_prime))
flow_rate, ice_density, u_init_arr = 1e-3, 910.0, h₀
n = 3.0
ρ = ice_density
g = 9.8101
A = fill(flow_rate, ne(s))
tₑ = 10.0
u₀ = ComponentArray(dynamics_h=u_init_arr)
constants_and_parameters = ComponentArray(
n=n,
stress_ρ=ρ,
stress_g=g,
stress_A=A)
data_prob = ODEProblem{true,SciMLBase.FullSpecialize}(fₘ, u₀, (0, tₑ), constants_and_parameters)
decapode_sol = solve(data_prob, Tsit5())
reference_dat = last(decapode_sol).dynamics_h
function loss(u) #only compares last time step
newp = ComponentArray(n=n, stress_ρ=u[1], stress_g=g, stress_A=A)
prob = remake(data_prob, p=newp)
sol = solve(prob, Rodas4(autodiff=false), sensealg=GaussAdjoint())
current_dat = last(sol).dynamics_h
sum(abs2, reference_dat .- current_dat)
end
#loss([800.0])
Zygote.gradient(loss, [800.0]) |
The code generation gives something like this decapode_f = begin
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:758 =#
(mesh, operators, hodge = GeometricHodge())->begin
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:758 =#
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:759 =#
begin
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:235 =#
(var"GenSim-M_d₀", d₀) = default_dec_matrix_generate(mesh, :d₀, hodge)
(var"GenSim-M_avg₀₁", avg₀₁) = default_dec_matrix_generate(mesh, :avg₀₁, hodge)
(var"GenSim-M_⋆₁", ⋆₁) = default_dec_matrix_generate(mesh, :⋆₁, hodge)
(var"GenSim-M_dual_d₀", dual_d₀) = default_dec_matrix_generate(mesh, :dual_d₀, hodge)
(var"GenSim-M_⋆₀⁻¹", ⋆₀⁻¹) = default_dec_matrix_generate(mesh, :⋆₀⁻¹, hodge)
♯ = operators(mesh, :♯)
mag = operators(mesh, :mag)
end
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:760 =#
begin
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:619 =#
var"GenSim-M_GenSim-ConMat_1" = var"GenSim-M_⋆₀⁻¹" * var"GenSim-M_dual_d₀" * var"GenSim-M_⋆₁"
var"GenSim-ConMat_1" = (x->var"GenSim-M_GenSim-ConMat_1" * x)
end
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:761 =#
begin
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:140 =#
var"__dynamics_•1" = Decapodes.FixedSizeDiffCache(Vector{Float64}(undef, nparts(mesh, :E)))
var"__dynamics_•6" = Decapodes.FixedSizeDiffCache(Vector{Float64}(undef, nparts(mesh, :E)))
__dynamics_mult_3 = Decapodes.FixedSizeDiffCache(Vector{Float64}(undef, nparts(mesh, :E)))
__dynamics_ḣ = Decapodes.FixedSizeDiffCache(Vector{Float64}(undef, nparts(mesh, :V)))
end
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:762 =#
f(du, u, p, t) = begin
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:762 =#
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:763 =#
begin
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:297 =#
dynamics_h = u.dynamics_h
n = p.n
stress_A = p.stress_A
stress_ρ = p.stress_ρ
stress_g = p.stress_g
var"1" = 1.0
var"2" = 2.0
var"2.0" = 2.0
end
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:764 =#
begin
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:534 =#
var"dynamics_•1" = Decapodes.get_tmp(var"__dynamics_•1", u)
var"dynamics_•6" = Decapodes.get_tmp(var"__dynamics_•6", u)
dynamics_mult_3 = Decapodes.get_tmp(__dynamics_mult_3, u)
dynamics_ḣ = Decapodes.get_tmp(__dynamics_ḣ, u)
end
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:765 =#
mul!(var"dynamics_•1", var"GenSim-M_d₀", dynamics_h)
mul!(var"dynamics_•6", var"GenSim-M_d₀", dynamics_h)
var"dynamics_•5" = ♯(var"dynamics_•6")
var"dynamics_•4" = mag(var"dynamics_•5")
var"dynamics_•7" = n .- var"1"
var"dynamics_•3" = var"dynamics_•4" .^ var"dynamics_•7"
var"stress_•3" = stress_ρ .* stress_g
var"stress_•2" = var"stress_•3" .^ n
dynamics_sum_1 = (.+)(n, var"2")
stress_sum_1 = (.+)(n, var"2.0")
var"dynamics_•2" = avg₀₁(var"dynamics_•3")
var"dynamics_•9" = dynamics_h .^ dynamics_sum_1
var"stress_•1" = var"2.0" ./ stress_sum_1
stress_mult_1 = var"stress_•1" .* stress_A
Γ = stress_mult_1 .* var"stress_•2"
var"dynamics_•8" = avg₀₁(var"dynamics_•9")
dynamics_mult_3 .= Γ .* var"dynamics_•1"
dynamics_mult_1 = dynamics_mult_3 .* var"dynamics_•2"
dynamics_mult_2 = dynamics_mult_1 .* var"dynamics_•8"
mul!(dynamics_ḣ, var"GenSim-M_GenSim-ConMat_1", dynamics_mult_2)
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:766 =#
begin
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:338 =#
setproperty!(du, :dynamics_h, dynamics_ḣ)
end
#= /home/jadonclugston/.julia/packages/Decapodes/MGJA6/src/simulation.jl:767 =#
return nothing
end
end
end where |
The issue comes from function reinterpret(::Type{T}, a::A) where {T,N,S,A<:AbstractArray{S, N}}
function thrownonint(S::Type, T::Type, dim)
@noinline
throw(ArgumentError("""
cannot reinterpret an `$(S)` array to `$(T)` whose first dimension has size `$(dim)`.
The resulting array would have non-integral first dimension.
"""))
end
function throwaxes1(S::Type, T::Type, ax1)
@noinline
throw(ArgumentError("cannot reinterpret a `$(S)` array to `$(T)` when the first axis is $ax1. Try reshaping first."))
end
isbitstype(T) || throwbits(S, T, T)
isbitstype(S) || throwbits(S, T, S)
(N != 0 || sizeof(T) == sizeof(S)) || throwsize0(S, T, "different")
if N != 0 && sizeof(S) != sizeof(T)
ax1 = axes(a)[1]
dim = length(ax1)
if issingletontype(T)
issingletontype(S) || throwsingleton(S, T)
else
rem(dim*sizeof(S),sizeof(T)) == 0 || thrownonint(S, T, dim)
end
first(ax1) == 1 || throwaxes1(S, T, ax1)
end
readable = array_subpadding(T, S)
writable = array_subpadding(S, T)
new{T, N, S, A, false}(a, readable, writable)
end
this line is being triggered: ~~but as far as I can tell from debugging, so I'm not sure why the error should be triggering~~ actually I was doing something wrong, I think so
|
|
More direct: using Catlab
using Catlab.Graphics
using CombinatorialSpaces
using Decapodes
using ComponentArrays
using ForwardDiff
using Zygote
using SciMLSensitivity
# External Dependencies
using MLStyle
using MultiScaleArrays
using LinearAlgebra
using OrdinaryDiffEq
using JLD2
using SparseArrays
using Statistics
using GeometryBasics: Point2, Point3
Point2D = Point2{Float64};
Point3D = Point3{Float64};
using DiagrammaticEquations
using DiagrammaticEquations.Deca
halfar_eq2 = @decapode begin
h::Form0
Γ::Form1
n::Constant
ḣ == ∂ₜ(h)
ḣ == ∘(⋆, d, ⋆)(Γ * d(h) * avg₀₁(mag(♯(d(h)))^(n - 1)) * avg₀₁(h^(n + 2)))
end
glens_law = @decapode begin
Γ::Form1
(A, ρ, g, n)::Constant
Γ == (2.0 / (n + 2.0)) * A * (ρ * g)^n
end
@info("Decapodes Defined")
ice_dynamics_composition_diagram = @relation () begin
dynamics(Γ, n)
stress(Γ, n)
end
ice_dynamics_cospan = oapply(ice_dynamics_composition_diagram,
[Open(halfar_eq2, [:Γ, :n]),
Open(glens_law, [:Γ, :n])])
ice_dynamics = apex(ice_dynamics_cospan)
ice_dynamics1D = expand_operators(ice_dynamics)
infer_types!(ice_dynamics1D, op1_inf_rules_1D, op2_inf_rules_1D)
resolve_overloads!(ice_dynamics1D, op1_res_rules_1D, op2_res_rules_1D)
s_prime = EmbeddedDeltaSet1D{Bool,Point2D}()
add_vertices!(s_prime, 10, point=Point2D.(range(-2.0, 2.0, length= 10), 0.0))
add_edges!(s_prime, 1:nv(s_prime)-1, 2:nv(s_prime))
orient!(s_prime)
s = EmbeddedDeltaDualComplex1D{Bool,Float64,Point2D}(s_prime)
subdivide_duals!(s, Circumcenter())
function generate(sd, my_symbol; hodge=GeometricHodge())
op = @match my_symbol begin
:♯ => x -> begin
# This is an implementation of the "sharp" operator from the exterior
# calculus, which takes co-vector fields to vector fields.
# This could be up-streamed to the CombinatorialSpaces.jl library. (i.e.
# this operation is not bespoke to this simulation.)
e_vecs = map(edges(sd)) do e
point(sd, sd[e, :∂v0]) - point(sd, sd[e, :∂v1])
end
neighbors = map(vertices(sd)) do v
union(incident(sd, v, :∂v0), incident(sd, v, :∂v1))
end
n_vecs = map(neighbors) do es
[e_vecs[e] for e in es]
end
map(neighbors, n_vecs) do es, nvs
sum([nv * norm(nv) * x[e] for (e, nv) in zip(es, nvs)]) / sum(norm.(nvs))
end
end
:mag => x -> norm.(x)
x => error("Unmatched operator $my_symbol")
end
return (args...) -> op(args...)
end
decapode_code = gensim(ice_dynamics1D, dimension=1, preallocate=false)
#file = open("ice_sheet1D_alloc.jl", "w")
#write(file, string("decapode_f = ", decapode_code))
#close(file)
include("ice_sheet1D_alloc.jl")
fₘ = decapode_f(s, generate)
h₀ = map(x -> exp(-2 * x[1]^2), point(s_prime))
flow_rate, ice_density, u_init_arr = 1e-3, 910.0, h₀
n = 3.0
ρ = ice_density
g = 9.8101
A = fill(flow_rate, ne(s))
tₑ = 10.0
u₀ = ComponentArray(dynamics_h=u_init_arr)
constants_and_parameters = ComponentArray(
n=n,
stress_ρ=ρ,
stress_g=g,
stress_A=A)
data_prob = ODEProblem{true,SciMLBase.FullSpecialize}(fₘ, u₀, (0, tₑ), constants_and_parameters)
function steploss(du,u) #only compares last time step
data_prob.f(du,u,data_prob.p,data_prob.tspan[1])
nothing
end
u = copy(data_prob.u0)
du = zero(u)
d_u = zero(u)
d_du = zero(u)
using Enzyme
Enzyme.autodiff(Reverse, steploss, Duplicated(du, d_du), Duplicated(u,d_u))
using Cthulhu
@descend data_prob.f(du,u,data_prob.p,data_prob.tspan[1]) |
@wsmoses it seems the sharp function here gives Enzyme some issues. :♯ => x -> begin
# This is an implementation of the "sharp" operator from the exterior
# calculus, which takes co-vector fields to vector fields.
# This could be up-streamed to the CombinatorialSpaces.jl library. (i.e.
# this operation is not bespoke to this simulation.)
e_vecs = map(edges(sd)) do e
point(sd, sd[e, :∂v0]) - point(sd, sd[e, :∂v1])
end
neighbors = map(vertices(sd)) do v
union(incident(sd, v, :∂v0), incident(sd, v, :∂v1))
end
n_vecs = map(neighbors) do es
[e_vecs[e] for e in es]
end
map(neighbors, n_vecs) do es, nvs
sum([nv * norm(nv) * x[e] for (e, nv) in zip(es, nvs)]) / sum(norm.(nvs))
end
end @jpfairbanks if I understand correctly, |
For debugging I started manually modifying the generated function: decapode_f = begin
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:758 =#
(mesh, operators, hodge = GeometricHodge())->begin
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:758 =#
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:759 =#
begin
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:235 =#
(var"GenSim-M_d₀", d₀) = default_dec_matrix_generate(mesh, :d₀, hodge)
(var"GenSim-M_avg₀₁", avg₀₁) = default_dec_matrix_generate(mesh, :avg₀₁, hodge)
(var"GenSim-M_⋆₁", ⋆₁) = default_dec_matrix_generate(mesh, :⋆₁, hodge)
(var"GenSim-M_dual_d₀", dual_d₀) = default_dec_matrix_generate(mesh, :dual_d₀, hodge)
(var"GenSim-M_⋆₀⁻¹", ⋆₀⁻¹) = default_dec_matrix_generate(mesh, :⋆₀⁻¹, hodge)
♯ = x -> begin
# This is an implementation of the "sharp" operator from the exterior
# calculus, which takes co-vector fields to vector fields.
# This could be up-streamed to the CombinatorialSpaces.jl library. (i.e.
# this operation is not bespoke to this simulation.)
e_vecs = map(edges(mesh)) do e
point(mesh, mesh[e, :∂v0]) - point(mesh, mesh[e, :∂v1])
end
neighbors = map(vertices(mesh)) do v
union(incident(mesh, v, :∂v0), incident(mesh, v, :∂v1))
end
n_vecs = map(neighbors) do es
[e_vecs[e] for e in es]
end
res = zeros(Point2D, length(n_vecs))
map!(res, neighbors, n_vecs) do es, nvs
sum([nv * norm(nv) * x[e] for (e, nv) in zip(es, nvs)]) / sum(norm.(nvs))
end
res
end
mag = x -> norm.(x)
end
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:760 =#
begin
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:619 =#
var"GenSim-M_GenSim-ConMat_1" = var"GenSim-M_⋆₀⁻¹" * var"GenSim-M_dual_d₀" * var"GenSim-M_⋆₁"
var"GenSim-ConMat_1" = (x->var"GenSim-M_GenSim-ConMat_1" * x)
end
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:761 =#
begin
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:140 =#
end
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:762 =#
let ♯=♯, mag=mag
f(du, u, p, t) = begin
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:762 =#
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:763 =#
begin
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:297 =#
dynamics_h = u.dynamics_h
n = p.n
stress_A = p.stress_A
stress_ρ = p.stress_ρ
stress_g = p.stress_g
var"1" = 1.0
var"2" = 2.0
var"2.0" = 2.0
end
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:764 =#
begin
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:534 =#
end
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:765 =#
var"dynamics_•1" = d₀(dynamics_h)
var"dynamics_•6" = d₀(dynamics_h)
var"dynamics_•5" = ♯(var"dynamics_•6")
var"dynamics_•4" = mag(var"dynamics_•5")
var"dynamics_•7" = n .- var"1"
var"dynamics_•3" = var"dynamics_•4" .^ var"dynamics_•7"
var"stress_•3" = stress_ρ .* stress_g
var"stress_•2" = var"stress_•3" .^ n
dynamics_sum_1 = (.+)(n, var"2")
stress_sum_1 = (.+)(n, var"2.0")
var"dynamics_•2" = avg₀₁(var"dynamics_•3")
var"dynamics_•9" = dynamics_h .^ dynamics_sum_1
var"stress_•1" = var"2.0" ./ stress_sum_1
stress_mult_1 = var"stress_•1" .* stress_A
Γ = stress_mult_1 .* var"stress_•2"
var"dynamics_•8" = avg₀₁(var"dynamics_•9")
dynamics_mult_3 = Γ .* var"dynamics_•1"
dynamics_mult_1 = dynamics_mult_3 .* var"dynamics_•2"
dynamics_mult_2 = dynamics_mult_1 .* var"dynamics_•8"
dynamics_ḣ = var"GenSim-ConMat_1"(dynamics_mult_2)
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:766 =#
begin
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:338 =#
setproperty!(du, :dynamics_h, dynamics_ḣ)
end
#= /Users/chrisrackauckas/.julia/packages/Decapodes/MGJA6/src/simulation.jl:767 =#
return nothing
end
end
end
end This is what is being differentiated and giving Enzyme the type analysis error. res = zeros(Point2D, length(n_vecs))
map!(res, neighbors, n_vecs) do es, nvs
sum([nv * norm(nv) * x[e] for (e, nv) in zip(es, nvs)]) / sum(norm.(nvs))
end
res
In-place map forces the return value, but that still errors the same. This means it's pretty isolated down to the constructs sd[e, :∂v0] and incident(sd, v, :∂v0) being the weird ones, so I'm pretty sure however they are working is something Enzyme disagrees with. |
@wsmoses I think it's a missing BLAS rule? ERROR: MethodError: no method matching augmented_primal(::EnzymeCore.EnzymeRules.RevConfigWidth{…}, ::Const{…}, ::Type{…}, ::Duplicated{…}, ::Duplicated{…}, ::Duplicated{…}, ::Const{…}, ::Const{…})
Closest candidates are:
augmented_primal(::EnzymeCore.EnzymeRules.RevConfig, ::Const{typeof(mul!)}, ::Type{RT}, ::Annotation{<:StridedVecOrMat}, ::Const{<:Union{SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, SubArray{Tv, 2, <:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, Tuple{Base.Slice{Base.OneTo{Int64}}, I}} where I<:AbstractUnitRange} where {Tv, Ti}}, ::Annotation{<:StridedVecOrMat}, ::Annotation{<:Number}, ::Annotation{<:Number}) where RT
@ Enzyme ~/.julia/packages/Enzyme/azJki/src/internal_rules.jl:732
augmented_primal(::EnzymeCore.EnzymeRules.RevConfig, ::Const{Type{BigFloat}}, ::Type{<:Union{BatchDuplicated, BatchDuplicatedNoNeed, Duplicated, DuplicatedNoNeed}}, ::Any...)
@ Enzyme ~/.julia/packages/Enzyme/azJki/src/internal_rules.jl:1404
augmented_primal(::Any, ::Const{typeof(QuadGK.quadgk)}, ::Type{RT}, ::Any, ::Annotation{T}...; kws...) where {RT, T}
@ QuadGKEnzymeExt ~/.julia/packages/QuadGK/BjmU0/ext/QuadGKEnzymeExt.jl:6
...
Stacktrace:
[1] custom_rule_method_error
@ ~/.julia/packages/Enzyme/azJki/src/rules/customrules.jl:452 [inlined]
[2] mul!
@ ~/.julia/juliaup/julia-1.10.6+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:237 [inlined]
[3] *
@ ~/.julia/juliaup/julia-1.10.6+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:57 [inlined]
[4] #15
@ ~/.julia/packages/Decapodes/MGJA6/src/operators.jl:99 [inlined]
[5] f
@ ~/Desktop/askem/ice_sheet1D_alloc.jl:72
[6] ODEFunction
@ ~/.julia/packages/SciMLBase/BWkqx/src/scimlfunctions.jl:2362 [inlined]
[7] ODEFunction
@ ~/.julia/packages/SciMLBase/BWkqx/src/scimlfunctions.jl:0 [inlined]
[8] augmented_julia_ODEFunction_16361_inner_1wrap
@ ~/.julia/packages/SciMLBase/BWkqx/src/scimlfunctions.jl:0
[9] macro expansion
@ ~/.julia/packages/Enzyme/azJki/src/compiler.jl:8369 [inlined]
[10] enzyme_call
@ ~/.julia/packages/Enzyme/azJki/src/compiler.jl:7932 [inlined]
[11] AugmentedForwardThunk
@ ~/.julia/packages/Enzyme/azJki/src/compiler.jl:7769 [inlined]
[12] runtime_generic_augfwd(activity::Type{…}, runtimeActivity::Val{…}, width::Val{…}, ModifiedBetween::Val{…}, RT::Val{…}, f::ODEFunction{…}, df::ODEFunction{…}, primal_1::ComponentVector{…}, shadow_1_1::ComponentVector{…}, primal_2::ComponentVector{…}, shadow_2_1::ComponentVector{…}, primal_3::ComponentVector{…}, shadow_3_1::ComponentVector{…}, primal_4::Float64, shadow_4_1::Base.RefValue{…})
@ Enzyme.Compiler ~/.julia/packages/Enzyme/azJki/src/rules/jitrules.jl:483
[13] steploss
@ ~/Desktop/askem/mwe.jl:126 [inlined]
[14] steploss
@ ~/Desktop/askem/mwe.jl:0 [inlined]
[15] diffejulia_steploss_11669_inner_1wrap
@ ~/Desktop/askem/mwe.jl:0
[16] macro expansion
@ ~/.julia/packages/Enzyme/azJki/src/compiler.jl:8369 [inlined]
[17] enzyme_call
@ ~/.julia/packages/Enzyme/azJki/src/compiler.jl:7932 [inlined]
[18] CombinedAdjointThunk
@ ~/.julia/packages/Enzyme/azJki/src/compiler.jl:7705 [inlined]
[19] autodiff
@ ~/.julia/packages/Enzyme/azJki/src/Enzyme.jl:491 [inlined]
[20] autodiff
@ ~/.julia/packages/Enzyme/azJki/src/Enzyme.jl:537 [inlined]
[21] autodiff(::ReverseMode{…}, ::typeof(steploss), ::Duplicated{…}, ::Duplicated{…})
@ Enzyme ~/.julia/packages/Enzyme/azJki/src/Enzyme.jl:504
[22] top-level scope
@ ~/Desktop/askem/mwe.jl:138
Some type information was truncated. Use `show(err)` to see complete types. Got this by adding: using Enzyme
Enzyme.API.strictAliasing!(false)
Enzyme.autodiff(Reverse, steploss, Duplicated(du, d_du), Duplicated(u,d_u)) https://github.com/JuliaLang/julia/blob/v1.10.6/stdlib/LinearAlgebra/src/matmul.jl#L237 |
Can you show the full types of the MethodError? |
julia> show(err)
1-element ExceptionStack:
LoadError: MethodError: no method matching augmented_primal(::EnzymeCore.EnzymeRules.RevConfigWidth{1, true, true, (false, true, true, false, false, false), false}, ::Const{typeof(mul!)}, ::Type{Duplicated{Vector{Float64}}}, ::Duplicated{Vector{Float64}}, ::Duplicated{SparseMatrixCSC{Int8, Int32}}, ::Duplicated{SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}, ::Const{Bool}, ::Const{Bool})
Closest candidates are:
augmented_primal(::EnzymeCore.EnzymeRules.RevConfig, ::Const{typeof(mul!)}, ::Type{RT}, ::Annotation{<:StridedVecOrMat}, ::Const{<:Union{SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, SubArray{Tv, 2, <:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, Tuple{Base.Slice{Base.OneTo{Int64}}, I}} where I<:AbstractUnitRange} where {Tv, Ti}}, ::Annotation{<:StridedVecOrMat}, ::Annotation{<:Number}, ::Annotation{<:Number}) where RT
@ Enzyme ~/.julia/packages/Enzyme/azJki/src/internal_rules.jl:732
augmented_primal(::EnzymeCore.EnzymeRules.RevConfig, ::Const{Type{BigFloat}}, ::Type{<:Union{BatchDuplicated, BatchDuplicatedNoNeed, Duplicated, DuplicatedNoNeed}}, ::Any...)
@ Enzyme ~/.julia/packages/Enzyme/azJki/src/internal_rules.jl:1404
augmented_primal(::Any, ::Const{typeof(QuadGK.quadgk)}, ::Type{RT}, ::Any, ::Annotation{T}...; kws...) where {RT, T}
@ QuadGKEnzymeExt ~/.julia/packages/QuadGK/BjmU0/ext/QuadGKEnzymeExt.jl:6
...
Stacktrace:
[1] custom_rule_method_error
@ ~/.julia/packages/Enzyme/azJki/src/rules/customrules.jl:452 [inlined]
[2] mul!
@ ~/.julia/juliaup/julia-1.10.6+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:237 [inlined]
[3] *
@ ~/.julia/juliaup/julia-1.10.6+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:57 [inlined]
[4] #15
@ ~/.julia/packages/Decapodes/MGJA6/src/operators.jl:99 [inlined]
[5] f
@ ~/Desktop/askem/ice_sheet1D_alloc.jl:72
[6] ODEFunction
@ ~/.julia/packages/SciMLBase/BWkqx/src/scimlfunctions.jl:2362 [inlined]
[7] ODEFunction
@ ~/.julia/packages/SciMLBase/BWkqx/src/scimlfunctions.jl:0 [inlined]
[8] augmented_julia_ODEFunction_17737_inner_1wrap
@ ~/.julia/packages/SciMLBase/BWkqx/src/scimlfunctions.jl:0
[9] macro expansion
@ ~/.julia/packages/Enzyme/azJki/src/compiler.jl:8369 [inlined]
[10] enzyme_call
@ ~/.julia/packages/Enzyme/azJki/src/compiler.jl:7932 [inlined]
[11] AugmentedForwardThunk
@ ~/.julia/packages/Enzyme/azJki/src/compiler.jl:7769 [inlined]
[12] runtime_generic_augfwd(activity::Type{Val{(true, true, true, true, true)}}, runtimeActivity::Val{false}, width::Val{1}, ModifiedBetween::Val{(true, true, true, true, true)}, RT::Val{@NamedTuple{1, 2, 3}}, f::ODEFunction{true, SciMLBase.FullSpecialize, var"#f#53"{var"#42#52"{SparseMatrixCSC{Float64, Int64}}, Decapodes.var"#37#38"{SparseMatrixCSC{Float64, Int32}}, Decapodes.var"#15#16"{SparseMatrixCSC{Int8, Int32}}, var"#34#44"{EmbeddedDeltaDualComplex1D{Bool, Float64, GeometryBasics.Point{2, Float64}}}, var"#41#51"}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, df::ODEFunction{true, SciMLBase.FullSpecialize, var"#f#53"{var"#42#52"{SparseMatrixCSC{Float64, Int64}}, Decapodes.var"#37#38"{SparseMatrixCSC{Float64, Int32}}, Decapodes.var"#15#16"{SparseMatrixCSC{Int8, Int32}}, var"#34#44"{EmbeddedDeltaDualComplex1D{Bool, Float64, GeometryBasics.Point{2, Float64}}}, var"#41#51"}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, primal_1::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(dynamics_h = 1:10,)}}}, shadow_1_1::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(dynamics_h = 1:10,)}}}, primal_2::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(dynamics_h = 1:10,)}}}, shadow_2_1::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(dynamics_h = 1:10,)}}}, primal_3::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(n = 1, stress_ρ = 2, stress_g = 3, stress_A = 4:12)}}}, shadow_3_1::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(n = 1, stress_ρ = 2, stress_g = 3, stress_A = 4:12)}}}, primal_4::Float64, shadow_4_1::Base.RefValue{Float64})
@ Enzyme.Compiler ~/.julia/packages/Enzyme/azJki/src/rules/jitrules.jl:483
[13] steploss
@ ~/Desktop/askem/mwe.jl:126 [inlined]
[14] steploss
@ ~/Desktop/askem/mwe.jl:0 [inlined]
[15] diffejulia_steploss_13475_inner_1wrap
@ ~/Desktop/askem/mwe.jl:0
[16] macro expansion
@ ~/.julia/packages/Enzyme/azJki/src/compiler.jl:8369 [inlined]
[17] enzyme_call
@ ~/.julia/packages/Enzyme/azJki/src/compiler.jl:7932 [inlined]
[18] CombinedAdjointThunk
@ ~/.julia/packages/Enzyme/azJki/src/compiler.jl:7705 [inlined]
[19] autodiff
@ ~/.julia/packages/Enzyme/azJki/src/Enzyme.jl:491 [inlined]
[20] autodiff
@ ~/.julia/packages/Enzyme/azJki/src/Enzyme.jl:537 [inlined]
[21] autodiff(::ReverseMode{false, false, FFIABI, false, false}, ::typeof(steploss), ::Duplicated{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(dynamics_h = 1:10,)}}}}, ::Duplicated{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(dynamics_h = 1:10,)}}}})
@ Enzyme ~/.julia/packages/Enzyme/azJki/src/Enzyme.jl:504
[22] top-level scope
@ ~/Desktop/askem/mwe.jl:138
[23] eval
@ ./boot.jl:385 [inlined]
[24] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:2076
[25] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::@Kwargs{})
@ Base ./essentials.jl:892
[26] invokelatest(::Any, ::Any, ::Vararg{Any})
@ Base ./essentials.jl:889
[27] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/eval.jl:271
[28] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/eval.jl:181
[29] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/repl.jl:276
[30] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/eval.jl:179
[31] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/repl.jl:38
[32] (::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/eval.jl:150
[33] with_logstate(f::Function, logstate::Any)
@ Base.CoreLogging ./logging.jl:515
[34] with_logger
@ ./logging.jl:627 [inlined]
[35] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/eval.jl:263
[36] #invokelatest#2
@ ./essentials.jl:892 [inlined]
[37] invokelatest(::Any)
@ Base ./essentials.jl:889
[38] (::VSCodeServer.var"#64#65")()
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/eval.jl:34
in expression starting at /Users/chrisrackauckas/Desktop/askem/mwe.jl:138 Looks like it is due to |
@jpfairbanks Can you highlight where the sparse matrix is used here? It looks like it's a missing rule in Enzyme EnzymeAD/Enzyme.jl#2013 |
If you do
FBDF()
the forwarddiff fails. The forwarddiff also fails if the sensealg usesautodiff = true
.The text was updated successfully, but these errors were encountered: