feat: initial implementation of SCCNonlinearProblem codegen

[AayushSabharwal]

I've checked this on the trivial case in SCCNonlinearSolve's tests. However, there is one major caveat: There can be multiple orderings of the SCCs. As a result, solving the same system with NonlinearProblem and SCCNonlinearProblem leads to different ordering of the solution vector, because one uses the order of unknowns in the system and the other uses the order of SCCs. I'm figuring out how to use the same sorted order of SCCs in the simplified system.

[AayushSabharwal]

Example:

julia> @variables u(t)[1:8] [irreducible = true]
julia> @mtkbuild sys = NonlinearSystem([
                 0 ~ u[5]^2 + u[4]
                 0 ~ u[3]^2 + u[5]
                 0 ~ cos(u[2]) - u[1]
                 0 ~ 2u[4] + u[3] + 1.0
                 0 ~ u[1] + u[2] + u[3] + u[4] + u[5]    + 2.0u[6] + 2.5u[7] + 1.5u[8]
                 0 ~ u[1] + u[2] + u[3] + 2.0u[4] + u[5] + 4.0u[6] - 1.5u[7] + 1.5u[8]
                 0 ~ u[1] + 2.0u[2] + 3.0u[3] + 5.0u[4] + 6.0u[5] + u[6] - u[7] - u[8]
                 0 ~ sin(u[1] + u[2]) + u[2]
             ])
julia> prob =  SCCNonlinearProblem{true}(sys, unknowns(sys) .=> zeros(8));
julia> sol = solve(prob, NewtonRaphson())
julia> sol
retcode: Success
u: 8-element Vector{Float64}:
 -0.443600343333633
  0.9032122743166665
 -0.41964337760708054
 -0.17610056436947866
 -0.6477988712610427
 -4.179650698308442
  1.0580471331304233
  2.204144548445586
julia> nprob = NonlinearProblem(sys, unknowns(sys) .=> zeros(8))
julia> nsol = solve(nprob, NewtonRaphson())
retcode: Success
u: 8-element Vector{Float64}:
 -0.41964337760708054
 -0.17610056436947866
 -0.6477988712610426
 -0.44360034333363063
  0.903212274316663
 -4.179650698308441
  1.0580471331304235
  2.204144548445586

Note how both solutions have the same answer, but in different permutations.

[ChrisRackauckas]

But SII should take care of that?

[ChrisRackauckas]

1. Is there a way to toggle the initialization problem to be SCC or not?

[AayushSabharwal]

This is rebased on top of #3235. That should be merged before this one. Also requires SciML/SciMLBase.jl#878 and SciML/SciMLBase.jl#878 to use SCCNonlinearProblem for initialization. However, doing so by default is breaking unless NonlinearSolve also pulls in SCCNonlinearSolve or MTK does.

[AayushSabharwal]

[email protected] was tagged a few hours ago, I don't know why CI can't find it

[ChrisRackauckas]

We should separately validate use_scc && neqs == nunknown and error saying that SCCs can only be used if fully determined, falling back to non-SCC codegen.

[AayushSabharwal]

Added the warnings in 9943532

[ChrisRackauckas]

what's this case?

[AayushSabharwal]

If the user provides all the initial conditions to ODEProblem, we don't construct an initialization system and thus the InitializationSystemMetadata above this is not available. The metadata is used to merge the u0 and p maps provided to remake with the ones used to construct the existing initialization problem since otherwise the initialization system for the remade problem is underdetermined. It's also used to populate guesses to prevent missing initial conditions for the nonlinear problem (if a guess was provided to ODEProblem, it wouldn't be available anywhere else without this). Since we don't have this metadata, we can't perform the merge, so the code here is a way of trying to figure out the extra initial conditions/guesses necessary.

1. If the user doesn't provide u0 to remake, they want to retain the current u0 state. Hence, we fill u0map with the values of differential variables allowing initialization to solve for the algebraic ones.
2. If the user did provide u0, we can't fill the differential variables because they might have provided the value of an algebraic variable to remake and expect us to solve for the differential ones. This gets tricky, because there's no guarantee the user provided enough initial conditions to fully determine the initialization and we don't have the information to figure out which of the previous u0 they want to retain. We can solve this by allowing the initialization to be underdetermined and setting guesses of all the unknowns. If initialization can't solve for any of them, they retain the guess value.
3. Similarly for p, if the user doesn't provide it to remake, they want to retain values so we provide values for all parameters.
4. If the user provided p to remake, we treat solvable parameters like unknowns and provide guesses (since we don't know which ones the user wants to solve for) and set values for the rest.

[AayushSabharwal]

Looks like this doesn't fully work in case 2. The extra unknowns aren't involved in the equations of the initialization system, so structural_simplify removes them and we throw an IncompleteInitializationError since the simplified system doesn't have all the unknowns of the ODE.

[AayushSabharwal]

Turns out this was the wrong solution to a different problem.

[AayushSabharwal]

SciMLBase/Downstream will fail because it uses its own Project.toml and doesn't dev this branch