add rebuild function

[isaacsas]

Since remake won't handle updating conserved constants in problems by design, we can make our own wrapper over it we tell Catalyst users to use that will be more user friendly, (say rebuild).

All it needs to do is wrap remake, check if a user is updating a species a conserved constant depends on, and then calculate the new constant value. Then we call remake adding the new constant value into the updated p passed by the user (or making an updated p if none is passed). If a user passes a manual value for the constant we can just forward everything to remake directly (which seems to work).

Then we can update the docs and tests to use rebuild instead of remake.

[ChrisRackauckas]

Since remake won't handle updating conserved constants in problems by design, we can make our own wrapper over it we tell Catalyst users to use that will be more user friendly, (say rebuild).

remake does?

[isaacsas]

I thought the whole discussion we just had in the other issue was that it is deferred to integrator initialization, so passing a fully specified set of independent parameters and initial conditions to ODEProblem does not result in a problem storing the same u0 and p as calling remake with the same input?

[ChrisRackauckas]

Yes it needs integrator initialization because generally the u0 and p requires resolving both callbacks and initialization. This has always been true in general though. For ODEs without callbacks the problem is described fully explicitly and so it's not generally an issue, but with some types of callbacks and with DAEs in general the system is only described implicitly and thus you need an implicit solve in order to resolve the initial state before the integration. Since it requires a solver, and problem constructors do not do solving, it is deferred until the start of the integration.

All it needs to do is wrap remake, check if a user is updating a species a conserved constant depends on, and then calculate the new constant value. Then we call remake adding the new constant value into the updated p passed by the user (or making an updated p if none is passed). If a user passes a manual value for the constant we can just forward everything to remake directly (which seems to work).

The issue though is that it's not really well-defined what the "right" value would be for a p that is a discrete quantity that is updated throughout the simulation. And since it generally changes, you always need to update p right at the start of the simulation for such quantities. So I presume here that you want to compute the p(0)?

Is the plan to effectively then call integ = init(prob,Rodas5P()) and then use the integ.u and integ.p to update a problem and return that?

[isaacsas]

I understand why one needs initialization.

I don't understand why ODEProblem(prob, u0, tspan, p) and remake(prob; u0, p) for the same u0 and p, i.e. when each is passed mappings that should fully-determine the parameters/initial conditions (for example, in the absence of callbacks and DAEs), are not generating problems with the same stored u0s and ps. If a given p value can be determined when creating an ODEProblem, and the same input mapping is passed to remake, why can that same p value not be determined in remake? Why is different initialization run in ODEProblem vs remake? I thought that point of remake is it should be equivalent to setting up a new ODEProblem (again, given the same input mappings)?

For the past year the advice to users has been, if you need to call f(du,u,p,t) yourself for some application, the recommended way to get at a correct p object is to build a problem yourself and then use the p it stores. There has been no other exported interface to get the new parameters structure that I am aware of. In that time, remake was the way to get a correctly updated p one could pass into that function too if you needed to change a parameter (again, I understand this will not be right if one has callbacks with initialization or DAEs but I don't see either of those as being particularly common in Catalyst to date -- there wasn't even an interface to add initialization to a callback defined via MTK previously).

[isaacsas]

Just to be clear; we are talking about the parameters representing conserved constants that arise from linear conservation laws among species, all of which must be defined via ODEs (not DAEs). How would our manually evaluating and setting such a p as part of remake have any issue with subsequent initialization? If it is problematic for remake to do this wouldn't it also be problematic that ODEProblem is currently (correctly) setting the value of this p from the u0 too?

We would now simply be initializing the value of such ps before calling remake to be the same as what ODEProblem currently generates.

[ChrisRackauckas]

How would our manually evaluating and setting such a p as part of remake have any issue with subsequent initialization? If it is problematic for remake to do this wouldn't it also be problematic that ODEProblem is currently (correctly) setting the value of this p from the u0 too?

In general solving for p cannot be done in isolation but requires also solving for u0. For example the conservation laws is a good example of this: the parameter must be equal to the sum of some values from u0. Thus in general, you have to solve a big nonlinear system in order to know what the u0 are and what the p's are, resolving all of the initial relations.

Now with the conservation laws and DAEs you have here, it just so happens that all of the DAE equations are linear and all of the relationships between the u0 and p are linear, so in theory you could build a matrix of this and do a single tearing/linear solve to learn the exact analytical solution to this relationship. But of course you can see how in general if someone specifies p ~ x + 1 (interpreted as x(0)) and x ~ p^2, you need to fully run the DAE initialization process. And the problem building state is supposed to be a non-solving stage, no dependencies on NonlinearSolve etc., so it can't be in ODEProblem, so that has to be in solve.

For the past year the advice to users has been, if you need to call f(du,u,p,t) yourself for some application, the recommended way to get at a correct p object is to build a problem yourself and then use the p it stores.

I see, that's a bit of an unfortunate. How common is that in the codebase and users? For all DAEs it won't generally give the derivative but a mix of derivatives and algebraic residuals so it's really only a useful call on ODEs.

[isaacsas]

For the past year the advice to users has been, if you need to call f(du,u,p,t) yourself for some application, the recommended way to get at a correct p object is to build a problem yourself and then use the p it stores.

I see, that's a bit of an unfortunate. How common is that in the codebase and users? For all DAEs it won't generally give the derivative but a mix of derivatives and algebraic residuals so it's really only a useful call on ODEs.

I have no idea. Prior to the past year, MTK itself told people to use varmap_to_vars. Then that stopped working and the way to get a valid parameter object if you need to use generated ODE functions yourself outside of SciML solvers was via a problem.

I thought one of the goals of MTK was to be composable with many Julia libraries, not just those in SciML? What if one wants to use an optimization package outside of Optimization.jl? Or a non-wrapped dynamical systems analysis library? Does initialization even work for all kinds of SciML integrators, or does it require them to opt into it? While making things work for DAEs is great new functionality, there is a large world of non-DAE models and tooling for which it would also be great to still have a relatively simple way to pass MTK generated functions and parameters. The BifurcationKit extension just directly builds a NonlinearFunction that is passed on -- it doesn't use any integrator initialization to setup guesses or parameters. That is an extension for which one would very much want to use a Catalyst system where conservation laws have been eliminated to ensure a non-singular Jacobian. (And that thankfully still works since the extension uses the older varmap_to_vars approach, which correctly evaluates the conservation law constants.)

But of course you can see how in general if someone specifies p ~ x + 1 (interpreted as x(0)) and x ~ p^2, you need to fully run the DAE initialization process.

Yes, I can see that and understand the broader issue in the context of DAEs and callback initialization. However, we are not discussing the general situation. The parameters we are discussing are introduced internally by Catalyst after system creation and are solely functions of species that are all defined via ODEs. A user can't use them in equations / initial conditions / callbacks / defining other parameters as they don't exist at the time the user is building their system. (We have even discussed attaching metadata to them to make clear they have a special meaning.)

Since integrator initialization with directly constructed ODEProblems works fine even with the current pre-initialization calculation, which correctly sets these parameters' values in the problem, I can't see what harm there would be in lightly wrapping remake to also evaluate them (and then explicitly pass through their calculated values to remake). Given that ODEProblem and varmap_to_vars are able to evaluate them without integrator initialization our wrapper should also be able to. Initialization would still run during subsequent integrator setup for any SciML integrator that supports it, ensuring that the broader collection of parameters and initial conditions is consistent.

[TorkelE]

I have to admit things are going around a bit in my mind with what should be possible and what shouldn't be possible. However, there are a couple of cases where

oprob = ODEProblem(rs, u0, tspan, p)
oprob = remake(oprob; u0 = u0_new, p = p_new)
sol = solve(oprob)

generates the wrong simulation. I think we more or less know which cases these are. So, for a starter, creating a remake which throws errors/warnings here seems like a direct improvement. Then we can try to expand it later on/update remake/whatever is needed. However, there seems to be a relatively straightforward thing to do here to make stuff more robust.

[ChrisRackauckas]

generates the wrong simulation.

Those have all been fixed now, right? Can you share any cases that are not fixed? My understanding is that all of those issues are now closed and your cases are now test cases.

[ChrisRackauckas]

I thought one of the goals of MTK was to be composable with many Julia libraries, not just those in SciML?

It is, how is it not? MTK is still compatible with the whole list. Though in order to support more advanced features, we do need to make sure that more solvers get updated to support the more advanced features. That said OrdinaryDiffEq is really the only set of solvers for which mass matrix DAEs are supported, except a few in ODEInterface.jl which are a bit slow (and will appropriately error in this case, and we should add the initializeprob/initializeprobmap handling when we get the time).

What if one wants to use an optimization package outside of Optimization.jl?

I don't see the issue here?

Does initialization even work for all kinds of SciML integrators, or does it require them to opt into it?

Only Sundials and OrdinaryDiffEq have full support for DAEs, and both support all of MTK. What is the issue here?

While making things work for DAEs is great new functionality, there is a large world of non-DAE models and tooling for which it would also be great to still have a relatively simple way to pass MTK generated functions and parameters.

Agreed. And we still do that? In fact, a lot of the recent changes were required for cases that @TorkelE put together for Catalyst.

The BifurcationKit extension just directly builds a NonlinearFunction that is passed on -- it doesn't use any integrator initialization to setup guesses or parameters. That is an extension for which one would very much want to use a Catalyst system where conservation laws have been eliminated to ensure a non-singular Jacobian.

We should improve that. @AayushSabharwal can you make sure the problem building here correctly uses the initializeprobmap? The homotopy continuation piece as well.

BTW Catalyst likely doesn't need the HomotopyContinuation extension with the MTK level one now?

(And that thankfully still works since the extension uses the older varmap_to_vars approach, which correctly evaluates the conservation law constants.)

It doesn't work for implicitly defined parameters, like I discussed above, which is why we're trying to move away from it. We need to put correctness first.

Yes, I can see that and understand the broader issue in the context of DAEs and callback initialization.

And Catalyst.jl where you define parameters in terms of other parameters, as coming from @TorkelE's examples that drove this in the first place. I don't understand the argument: this is required even for just ODEs, even for Catalyst to be correct! How can you resolve a nonlinear system of equations without calling a solver for a nonlinear system of equations? You can end up in that case even in Catalyst, as shown from the issues that drove the development in the first place!

Given that ODEProblem and varmap_to_vars are able to evaluate them without integrator initialization our wrapper should also be able to.

The explicit cases, but it still gets the implicit cases wrong so we shouldn't rely on that.

If you want to do this right, then what you'd do is: https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/lib/OrdinaryDiffEqCore/src/initialize_dae.jl#L152-L191. In fact, it may make sense to move OverrideInit and CheckInit to DiffEqBase so that the code can more easily be moved. But basically what you want to do is solve the prob.f.initializeprob and then integrator.u = prob.f.initializeprobmap(nlsol) integrator.p = prob.f.initializeprobpmap(prob, nlsol) are the u0 and p that satisfy the initialization. Note that since it can in theory be numerical, you do need to check the residual of the numerical solve and figure out what to do if the set of equations is not solvable (error?)

[AayushSabharwal]

I have no idea. Prior to the past year, MTK itself told people to use varmap_to_vars. Then that stopped working and the way to get a valid parameter object if you need to use generated ODE functions yourself outside of SciML solvers was via a problem.

varmap_to_vars still works just the way it did with v8 (with some bug fixes encountered along the way). I intentionally haven't touched it even with the new problem construction rewrite because existing users rely on its old semantics. If you want a valid parameter object with a split = true system, you don't need to construct a problem. MTKParameters is an exported and documented constructor.

As for remake and ODEProblem not being identical, they aren't. Since forever, if you only pass u0 to remake, it will not update p (and vice-versa). The example you gave in the MTK issue for remake not updating the parameter was found to be a bug, and is something I have on the list of things I need to fix. If you have dependent defaults and pass the relevant keyword to remake, it will respect those defaults and re-create the initialization system so that regardless of remake behavior/bugs, the solve will respect your initial conditions or retcode saying they can't be satisfied.

We should improve that. @AayushSabharwal can you make sure the problem building here correctly uses the initializeprobmap?

Sure.

The homotopy continuation piece as well.

The homotopy continuation interface is purely a NonlinearProblem alternative, it doesn't solve initialization?

[ChrisRackauckas]

But nonlinearproblem would have to resolve the parameters? So would homotopy continuation then?

[isaacsas]

Given that ODEProblem and varmap_to_vars are able to evaluate them without integrator initialization our wrapper should also be able to.

The explicit cases, but it still gets the implicit cases wrong so we shouldn't rely on that.

So perhaps what is needed is detection if everything can be evaluated explicitly, and if so that can be done without having to have the full initialization process? It seems like that is what ODEProblem currently does but not remake? (Or perhaps the remake issue is a bug from what @AayushSabharwal said -- I am not sure what is supposed to always work outside of integrator initialization now after all the options / cases that have been discussed across these issues.)

[isaacsas]

Here is a simple example showing that initialization is not working except in a newly created ODEProblem:

using Catalyst, SciMLBase, OrdinaryDiffEqTsit5, Test
rn = @reaction_network begin
       (k1,k2), A <--> B
       end
osys = convert(ODESystem, rn; remove_conserved = true)
osys = complete(osys)

oprob = ODEProblem(osys, [osys.A => 1.0, osys.B => 1.0], (0.0, 1.0), [osys.k1 => 2.0, osys.k2 => 3.0])
sol = solve(oprob, Tsit5())
# these pass
@test sol(1.0; idxs = osys.A + osys.B) ≈ 2
@test sol(1.0; idxs = osys.Γ[1]) ≈ 2

oprob2 = remake(oprob; u0 = [osys.B => 5.0]) 
sol2 = solve(oprob2, Tsit5())
# these fail
@test sol2(1.0; idxs = osys.A + osys.B) ≈ 6.0
@test sol2(1.0; idxs = osys.Γ[1]) ≈ 6.0

oprob3 = ODEProblem(osys, [osys.A => 1.0, osys.B => 5.0], (0.0, 1.0), [osys.k1 => 2.0, osys.k2 => 3.0])
sol3 = solve(oprob3, Tsit5())
# these pass
@test sol3(1.0; idxs = osys.A + osys.B) ≈ 6.0
@test sol3(1.0; idxs = osys.Γ[1]) ≈ 6.0

oprob4 = remake(oprob; u0 = [osys.B => 5.0], use_defaults = true) 
sol4 = solve(oprob4, Tsit5())
# these fail
@test sol4(1.0; idxs = osys.A + osys.B) ≈ 6.0
@test sol4(1.0; idxs = osys.Γ[1]) ≈ 6.0

[isaacsas]

oprob5 = remake(oprob; u0 = [osys.B => 5.0], p = [osys.Γ[1] => 6.0]) 
sol5 = solve(oprob5, Tsit5())
# these pass
@test sol5(1.0; idxs = osys.A + osys.B) ≈ 6.0
@test sol5(1.0; idxs = osys.Γ[1]) ≈ 6.0

So if we just wrapped remake with a function that calculates osys.Γ[1] and forwards it to remake we could have this all work.

edit: I guess it hasn't been a year since this particular issue was discovered, only six months, so I take that back. My apologies. I was conflating it with various other things that broke last fall and have since been fixed or changed.

[isaacsas]

And here is what NonlinearSolve does

nsys = convert(NonlinearSystem, rn; remove_conserved = true)
nsys = complete(nsys)

nprob = NonlinearProblem(nsys, [nsys.A => 1.0, nsys.B => 1.0], [nsys.k1 => 2.0, nsys.k2 => 3.0])
sol = solve(nprob)
# these pass
@test sol[nsys.A + nsys.B] ≈ 2
@test sol.ps[nsys.Γ[1]] ≈ 2

nprob2 = remake(nprob; u0 = [nsys.B => 5.0]) 
sol2 = solve(nprob2)
# these fail
@test sol2[nsys.A + nsys.B] ≈ 6.0
@test sol2.ps[nsys.Γ[1]] ≈ 6.0

nprob3 = NonlinearProblem(nsys, [nsys.A => 1.0, nsys.B => 5.0], [nsys.k1 => 2.0, nsys.k2 => 3.0])
sol3 = solve(nprob3)
# these pass
@test sol3[nsys.A + nsys.B] ≈ 6.0
@test sol3.ps[nsys.Γ[1]] ≈ 6.0

nprob4 = remake(nprob; u0 = [nsys.B => 5.0], use_defaults = true) 
sol4 = solve(nprob4)
# these fail
@test sol4[nsys.A + nsys.B] ≈ 6.0
@test sol4.ps[nsys.Γ[1]] ≈ 6.0

nprob5 = remake(nprob; u0 = [nsys.B => 5.0], p = [nsys.Γ[1] => 6.0]) 
sol5 = solve(nprob5)
# these pass
@test sol5[nsys.A + nsys.B] ≈ 6.0
@test sol5.ps[nsys.Γ[1]] ≈ 6.0

[AayushSabharwal]

But nonlinearproblem would have to resolve the parameters? So would homotopy continuation then?

We don't have initializeprob infrastructure anywhere except for ODEs/DAEs though. I can add that, but then every solver library needs to hook into it? generate_initializesystem also just assumes a time-dependent system.

Here is a simple example showing that initialization is not working except in a newly created ODEProblem:

There's two parts to the problem here:

• Remake won't touch parameters because you're not passing p. It never used to, this is just how it works. If you want it to solve for Gamma after you do remake without passing p, you need to give Gamma a guess and it'll solve for it when the problem is initialized.
• The fact that Gamma is still not updated even if you pass p = Dict() is a bug (assuming Gamma[1] => A + B is a default)

Also, could you check the retcode of the failing ODE solves? I have a feeling it'll be InitialFailure.

[isaacsas]

Thanks for the suggestions. Trying to use guesses here are the updated results -- none seem to work (I assume NonlinearProblem will have analogous issues). I didn't try passing guesses to the system during construction as that isn't a use-case that we can handle seemlessly for users (i.e. Catalyst won't necessarily know the appropriate type information to pass a guess to the system). I tried passing guesses to ODEProblem, remake, and solve but they have no-impact (the first) or give a kwarg error (the latter two). If passing guesses to remake worked we could just make a simple rebuild for now in Catalyst that wraps remake with an added guess on gamma to get this out quickly.

using Catalyst, SciMLBase, OrdinaryDiffEqTsit5, Test, NonlinearSolve
rn = @reaction_network begin
       (k1,k2), A <--> B
       end
osys = convert(ODESystem, rn; remove_conserved = true)
osys = complete(osys)

oprob = ODEProblem(osys, [osys.A => 1.0, osys.B => 1.0], (0.0, 1.0), [osys.k1 => 2.0, osys.k2 => 3.0])
sol = solve(oprob, Tsit5())
# these pass
@test sol(1.0; idxs = osys.A + osys.B) ≈ 2
@test sol(1.0; idxs = osys.Γ[1]) ≈ 2

oprob2 = remake(oprob; u0 = [osys.B => 5.0]) 
sol2 = solve(oprob2, Tsit5())
# these fail
@test sol2(1.0; idxs = osys.A + osys.B) ≈ 6.0
@test sol2(1.0; idxs = osys.Γ[1]) ≈ 6.0

oprob3 = ODEProblem(osys, [osys.A => 1.0, osys.B => 5.0], (0.0, 1.0), [osys.k1 => 2.0, osys.k2 => 3.0])
sol3 = solve(oprob3, Tsit5())
# these pass
@test sol3(1.0; idxs = osys.A + osys.B) ≈ 6.0
@test sol3(1.0; idxs = osys.Γ[1]) ≈ 6.0

oprob4 = remake(oprob; u0 = [osys.B => 5.0], use_defaults = true) 
sol4 = solve(oprob4, Tsit5())
# these fail
@test sol4(1.0; idxs = osys.A + osys.B) ≈ 6.0
@test sol4(1.0; idxs = osys.Γ[1]) ≈ 6.0

oprob5 = remake(oprob; u0 = [osys.B => 5.0], p = [osys.Γ[1] => 6.0]) 
sol5 = solve(oprob5, Tsit5())
# these pass
@test sol5(1.0; idxs = osys.A + osys.B) ≈ 6.0
@test sol5(1.0; idxs = osys.Γ[1]) ≈ 6.0


# guesses not defined for remake
oprob6 = remake(oprob; u0 = [osys.B => 5.0], guesses = [osys.Γ[1] => 1.0]) 
# errors about invalid kwarg
sol6 = solve(oprob6, Tsit5())

oprob7 = ODEProblem(osys, [osys.A => 1.0, osys.B => 5.0], (0.0, 1.0), 
    [osys.k1 => 2.0, osys.k2 => 3.0], guesses = [osys.Γ[1] => 1.0])
sol7 = solve(oprob7, Tsit5())
# these pass
@test sol7(1.0; idxs = osys.A + osys.B) ≈ 6.0
@test sol7(1.0; idxs = osys.Γ[1]) ≈ 6.0


oprob8 = remake(oprob7; u0 = [osys.B => 7.0])
sol8 = solve(oprob8, Tsit5())
# errors
@test sol8(1.0; idxs = osys.A + osys.B) ≈ 8.0
@test sol8(1.0; idxs = osys.Γ[1]) ≈ 8.0

oprob9 = remake(oprob7; u0 = [osys.B => 7.0])
# errors about invalid kwarg
sol9 = solve(oprob9, Tsit5(); guesses = [osys.Γ[1] => 1.0])

oprob10 = remake(oprob7; u0 = [osys.B => 7.0], use_defaults = true)
sol10 = solve(oprob10, Tsit5())
# errors
@test sol10(1.0; idxs = osys.A + osys.B) ≈ 8.0
@test sol10(1.0; idxs = osys.Γ[1]) ≈ 8.0

All the cases where solve runs, i.e. when not passing it an invalid kwarg, report success.

[isaacsas]

I'm not sure why remake's historically not updating p in this case is a reason to keep that behavior? If anything, @TorkelE and my repeated questions about this shows that some users might find this convention confusing. As far as I understand, it also means that for this simple explicitly defined parameter case, the integrator has to opt into initialization to ultimately get the right value, instead of the p it uses already being correct. That then limits the scope of where conservation law elimination can be used.

[ChrisRackauckas]

As far as I understand, it also means that for this simple explicitly defined parameter case, the integrator has to opt into initialization to ultimately get the right value, instead of the p it uses already being correct.

Proving explicit can be really difficult in general though. We don't have a good cycle detection for that. And then you also cannot make good guarantees either, since if someone does make a cycle then you would have to run the solve to resolve the parameters (and for a nonlinear problem, it is part of the nonlinear solve). Currently the cycle detection uses structural simplify and should make sure that we reuse some of the internals there for pre-simplification if possible?

If passing guesses to remake worked we could just make a simple rebuild for now in Catalyst that wraps remake with an added guess on gamma to get this out quickly.

See what I mentioned above about that. NonlinearProblem, OptimizationProblem, BifurcationProblem, HomotopyProblem etc. all need to get some upgrades in order to ensure correctness with respect to these kinds of constraints. For NonlinearProblem, it would either need to be a two-stage nonlinear solve with a pre-solve for the parameters and then a solve for the state, or you could lift the parameters into the state and let it fully simplify together. The pre-solve is likely more efficient though for the same reason SCCs are. For Optimization, you could just extend the constraint system and lift parameters into the state. For HomotopyProblem and BifurcationProblem, you can use the same initializationproblem pieces and run that prior to BifurcationProblem/HomotopyProblem.

[AayushSabharwal]

Thanks for the suggestions. Trying to use guesses here are the updated results -- none seem to work

I don't see where you've provided guesses (except at the very end, to ODEProblem) though? They're part of the system, like defaults. remake isn't at the level of abstraction of guesses, which is why you can't pass them to it.

oprob7 = ODEProblem(osys, [osys.A => 1.0, osys.B => 5.0], (0.0, 1.0), 
   [osys.k1 => 2.0, osys.k2 => 3.0], guesses = [osys.Γ[1] => 1.0])

oprob8 = remake(oprob7; u0 = [osys.B => 7.0])
sol8 = solve(oprob8, Tsit5())

@ChrisRackauckas this is an interesting case. Gamma becomes a solvable parameter when constructing ODEProblem, because a guess is provided for it. When remake is called, that guess is not available because it is not part of the system and hence Gamma isn't solvable anymore and the initprob is infeasible (A ~ 1, B ~ 7, A + B ~ Gamma but Gamma is a parameter). I guess we need to propagate this information so that Gamma is solvable in the remade initprob as well?

I'm not sure why remake's historically not updating p in this case is a reason to keep that behavior?

Because not every use case needs or wants to update p. Also, semver. We can't change what remake does that dramatically without a breaking release.

[AayushSabharwal]

For Optimization, you could just extend the constraint system and lift parameters into the state.

Do you mean doing this lifting during simplification, or actually turning them into unknowns? I don't quite get why we'd do the latter since it's an initialization process. Otherwise wouldn't it also make sense to lift them into states for ODEs and the like? Similarly for NonlinearProblems, initialization is not identical to SCCs since we're solving for the initial values of states that are also solved by the actual nonlinear problem, as opposed to it being a disjoint set.

[isaacsas]

If passing a quantitatively meaningless guess to a system solves this problem -- anything we pass would be no better than calling rand() or just passing 1 or 0 -- that seems suggestive that there should instead be a more general configuration option a system can be passed to tell it to opt-into whatever behavior that guess is enabling (via system metadata? parameter metadata, maybe issolvable?). It would also mean we don't have to start making assumptions about the type of user unknowns/parameters in Catalyst (because we may not have any defaults available to look at for deciding about types).

However, more generally we need the remake process in this case to work for all problem types, not just ODESystems that sit on top of integrators that opt into initialization. Otherwise we will experience the same issues with other problem types.

FWIW, whatever is being done in process_SciMLProblem right now seems to handle this correctly given that the initially constructed problem does get correct values for the conserved constants without needing the integrator initialization step or guesses for both ODEProblems and NonlinearProblems. Having an opt-into MTK system setting that tells remake to work more like that would be another non-breaking option perhaps?

[ChrisRackauckas]

FWIW, whatever is being done in process_SciMLProblem right now seems to handle this correctly given that the initially constructed problem does get correct values for the conserved constants without needing the integrator initialization step or guesses for both ODEProblems and NonlinearProblems. Having an opt-into MTK system setting that tells remake to work more like that would be another non-breaking option perhaps?

Well there's a few things that we could do. I do think it would be nice if all explicit relationships could be resolved immediately. That's always going to be slower than using a generated function though since it would require doing the symbolic substitution which is a slow interpreted process. Also, it would only apply to explicit relationships, so we'd need to error if there's any algebraic loops.

Do you mean doing this lifting during simplification, or actually turning them into unknowns? I don't quite get why we'd do the latter since it's an initialization process. Otherwise wouldn't it also make sense to lift them into states for ODEs and the like?

That's different, because the ODE is not a single implicit system but instead N+1 implicit systems.

Similarly for NonlinearProblems, initialization is not identical to SCCs since we're solving for the initial values of states that are also solved by the actual nonlinear problem, as opposed to it being a disjoint set.

It is a manual SCC decomposition. You are solving for the p's and then solving for the u's. If you put the two together, then the SCC decomposition would get one (or many) SCC for the p's, and then the next SCCs would be the u's, since those depend on the p's. You can then just represent the initial value of every u as 0 and then shift by the u0 set as parameters.

So it is equivalent. The reason to do it as two separate steps would just be an optimization over that.

[AayushSabharwal]

If passing a quantitatively meaningless guess

It isn't always meaningless; the guess can influence the result of the nonlinear solve. All zero guesses cause singularities. We require a guess for parameters for two main reasons:

• Opting into parameter initialization should be opt-in. People don't usually expect their parameter value can change when they call solve
• If they do want to opt-in, the equation may be implicit and we can't just take the guess from substituting the default

You can also opt-in to solving for a parameter by giving it a default of missing or passing a value of missing to ODEProblem. This avoids the requirement of a guess, but then whichever of the two is not missing is required to use as a guess.

It would also mean we don't have to start making assumptions about the type of user unknowns/parameters in Catalyst

What assumptions are you referring to here? The type of a parameter is always known. We validate this in MTK. You can't give a scalar parameter an array value, or vice versa, etc.

However, more generally we need the remake process in this case to work for all problem types, not just ODESystems that sit on top of integrators that opt into initialization.

remake works identically for all problem types. If it doesn't, it's likely a bug. remake doesn't change keywords it isn't passed because e.g. if you're optimizing over an ODE solve and passing a vector of initial values or directly passing a parameter object, you don't want the overhead of symbolic substitution. The initial conditions are always the same, and are handled by a nonlinear solve during integrator initialization.

FWIW, whatever is being done in process_SciMLProblem right now seems to handle this correctly given that the initially constructed problem does get correct values for the conserved constants without needing the integrator initialization step or guesses for both ODEProblems and NonlinearProblems.

process_SciMLProblem just does the substitutions. remake is similar but not identical to the ODEProblem constructor.

Having an opt-into MTK system setting that tells remake to work more like that would be another non-breaking option perhaps?

p = Dict() is already an opt-in? Also, MTK doesn't tell remake anything. If it did, that would be an implicit circular dependency and something we're actively trying to avoid ever since [email protected].

[AayushSabharwal]

It is a manual SCC decomposition. You are solving for the p's and then solving for the u's. If you put the two together, then the SCC decomposition would get one (or many) SCC for the p's, and then the next SCCs would be the u's, since those depend on the p's. You can then just represent the initial value of every u as 0 and then shift by the u0 set as parameters.

So for time-independent types, initialization should only solve for parameters?

[ChrisRackauckas]

and initial guesses.

[ChrisRackauckas]

I was just mentioning that presolving for initial guesses can be made equivalent to a parameter solve and a linear shift in the equations + SCC decomposition.

[AayushSabharwal]

I see. Yeah that makes sense. Thanks for the clarification.