Skip to content

Capture observed dependence on unknowns in simplified nonlinear systems #2863

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

Merged
merged 12 commits into from
Jul 16, 2024
Next Next commit
Substitute in observed before calculating jacobian
  • Loading branch information
hersle committed Jul 16, 2024
commit f736a7857a883f4e46eee2ed0b6d79ee2210a63c
5 changes: 4 additions & 1 deletion src/systems/nonlinear/nonlinearsystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -193,8 +193,11 @@ function calculate_jacobian(sys::NonlinearSystem; sparse = false, simplify = fal
return cache[1]
end

rhs = [eq.rhs for eq in equations(sys)]
# observed equations may depend on unknowns, so substitute them in first
obs = map(eq -> eq.lhs => eq.rhs, observed(sys))
rhs = map(eq -> substitute(eq.rhs, obs), equations(sys))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This won't catch a chain of observed, it would need a fixed point sub I think

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Like in the most recent commit? Can you think of a simple example on to test?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

z ~ y, y ~ x, x^2 - 2z` should hit it, since if you just plug in z you get y.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not quite after structural simplification, it seems.
But I added 0 ~ x^2 + 2*z + y, z ~ y, y ~ x, which seems to do the trick.

vals = [dv for dv in unknowns(sys)]

if sparse
jac = sparsejacobian(rhs, vals, simplify = simplify)
else
Expand Down