Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

chore: set indexing style #385

Closed
wants to merge 1 commit into from

Conversation

DhairyaLGandhi
Copy link
Member

@DhairyaLGandhi DhairyaLGandhi commented Jun 5, 2024

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

ReverseDiff restricts array types to IndexLinear() indexing style as written in their docs: https://juliadiff.org/ReverseDiff.jl/dev/limits/

Not following this causes failures with using ReverseDiff:

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

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

@mtkbuild sys = ODESystem(eqs, t)

ModelingToolkit.observed(sys)

u0 = [D(x) => 2.0,
    x => 1.0,
    y => 0.0,
    z => 0.0,]

p ==> 28.0,
    ρ => 10.0,
    β => 8 / 3,]

tspan = (0.0, 100.0)
prob = ODEProblem(sys, u0, tspan, p, jac = true)
sol = solve(prob, Tsit5())

mtkp = parameter_values(sol)
gt = rand(873)

gradient(mtkparams) do p
    sol = solve(prob, Rosenbrock23(), p = p, sensealg = SciMLSensitivity.GaussAdjoint(autojacvec = ReverseDiffVJP()))
    mean(abs.(sol[sys.x] .- gt))
end

cc @AayushSabharwal

Add any other context about the problem here.

@ChrisRackauckas
Copy link
Member

But its index style is naturally Cartesian?

@DhairyaLGandhi
Copy link
Member Author

DhairyaLGandhi commented Jun 5, 2024

Ah, that is right, maybe we can restrict it to the case of Vector of Vector? The reason is that when we call canonicalize we get an ArrayPartition, and that runs into this ReverseDiff limitation as I diagnosed this today with @AayushSabharwal. It runs into this assert https://github.com/JuliaDiff/ReverseDiff.jl/blob/bff801b75ea23ead2ece6f594ce1deee79898ba4/src/tracked.jl#L77

@ChrisRackauckas
Copy link
Member

Where does ReverseDiff come into play?

@DhairyaLGandhi
Copy link
Member Author

DhairyaLGandhi commented Jun 5, 2024

This is when we pass sensealg = GaussAdjoint(autojacvec = ReverseDiffVJP())

gradient(mtkparams) do p
    sol = solve(prob, Rosenbrock23(), p = p, sensealg = SciMLSensitivity.GaussAdjoint(autojacvec = ReverseDiffVJP()))
    mean(abs.(sol[sys.x] .- gt))
end

@ChrisRackauckas
Copy link
Member

With Array or ComponentArray? Is there a case that regresses/breaks or is it that case is not handled for the newer stuff?

@DhairyaLGandhi
Copy link
Member Author

It's not handled when we integrate SciMLStructures and SII with Sensitivity.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants