fix: ensure lower order derivative SCCs occur before higher order ones

[AayushSabharwal]

No description provided.

[oscardssmith]

why is this a necessary constraint?

[AayushSabharwal]

Semantically, because higher order derivatives depend on lower order ones. Practically, because it breaks BLT ordering. Structural passes only deal in the most differentiated variables, but for reassembly (this piece of code) we need all levels of differentiation in the SCCs. So this function is responsible for inserting them in. E.g. var_sccs will typically only contain an entry for D(D(D(x))). So this pass finds the SCC containing D(D(D(x))), says "we need to insert D(D(x)) and D(x) here". Depending on variable ordering, it might issue the insertion commands in either order, which is the problem this PR tackles.

We generate equations SCC-by-SCC, in BLT sorted order, in generate_system_equations!. The dictionary total_sub keeps track of all differential equations we've generated so far in this process, to substitute into subsequent equations. So e.g. if we generate

D(x) ~ 2x + t

then total_sub will contain D(x) => 2x + t. A subsequent SCC might then have the equation D(x) + D(y) ~ 3 matched to D(y). So ordinarily it would generate D(y) ~ 3 - D(x), but that's not how MTK works. total_sub is substituted in so that we get D(y) ~ 3 - 2x - t. If we don't have SCCs correctly sorted, total_sub won't contain the right values by the time we generate the D(y) equation. This often manifests itself as an error "Unexpected Differential(t, 1) in system with independent variable nothing" when creating an ODEProblem from the system. This is usually because there is an observed equation containing D(x) on the RHS which didn't get substituted away and ended up in the initialization system, which of course doesn't have an independent variable and complains if it sees a Differential operator.

In this case, the bug manifested as an error in this typeassert. What it means is that the inline linear SCC pass generated an A and b, where an equation eq is matched by tearing to be solved for variable v, but the corresponding entry A[eq, v] is zero. A little digging revealed this is a 1x1 SCC of the form 0 ~ D(D(x)) - x_tt which should have become D(x_t) ~ x_tt but didn't. I set Claude on figuring out why, and (eventually) it made these changes.

It's often very tempting to preprocess all equations and populate total_sub first, or populate total_sub while generating equations and substitute them all at the end. However, every single time I've found a bug where total_sub didn't contain the right entry, it is due to bad SCC ordering. To this end (and because Claude loves to try both of the aforementioned "solutions") I added this comment. Basically, total_sub is never wrong. Some structural information is wrong, which causes total_sub to be populated incorrectly. Despite this, though, it took some prodding to get Claude to abandon bad ideas.