Offset fix for ctarget matrix equations

[cadojo]

Overview

• Closes Matrix-output functions for CTarget build_function calls #735
• Adds one helper function, get_symnumber (just a copy of existing code, moved to a function)
• Adds some offset math for handling matrix equations in _build_function(::CTarget).

More Detail

Before this PR (master), calling the following code would result in an error.

julia> @variables x[1:2,1:2] y[1:4]
julia> eqs = x[:] .~ y

julia> build_function(eqs, y, target=ModelingToolkit.CTarget())
ERROR: MethodError: no method matching +(::Nothing, ::Int64)
Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...) at operators.jl:538
  +(::ChainRulesCore.One, ::Any) at /home/joe/.julia/packages/ChainRulesCore/cpHLu/src/differential_arithmetic.jl:94
  +(::ChainRulesCore.Zero, ::Any) at /home/joe/.julia/packages/ChainRulesCore/cpHLu/src/differential_arithmetic.jl:63

Upon further inspection, this issue was caused on master:build_function.jl:500.

# build_function.jl:500
i = findfirst(x->isequal(tosymbol(x isa Sym ? x : operation(x), escape=false), tosymbol(var, escape=false)),varordering)

Because the none of the output variables x appear in the state vector provided (here the state vector provided is the second argument to build_function, y), the findfirst call returns nothing.

Changes

• I used the above findfirst code in one new place, so I moved it to a new function called get_symnumber, and replaced that line (line 500) with a call to get_symnumber
• In numbered_expr(::Equation...), I changed the lhs and rhs offsets to the following logic:
  • If the output variable is not found in the provided state vector, lhs offset is equal to 1 + the index of the output vector element, and rhs offset is equal to the index of each state vector element
  • If the output variable is found in the provided state vector, then the previous logic was used
• In _build_function(::CTarget...), I added the following logic at the top:
  • If no output is in the provided state vector, use new logic (offset = index of output - 2)
  • If all outputs are in the provided state vector, use previous logic (offset = -1)

Conclusion

I'm definintely open to other implementations, but this was the least-optrusive solution I could think of. With this PR, the output of the following code block is as expected.

julia> @variables x[1:2,1:2] y[1:4]
julia> eqs = x[:] .~ y

julia> build_function(eqs, y, target=ModelingToolkit.CTarget())
"void diffeqf(double* du, double* RHS1) {\n  du[0] = RHS1[0];\n  du[1] = RHS1[1];\n  du[2] = RHS1[2];\n  du[3] = RHS1[3];\n}\n"

[ChrisRackauckas]

Can this get a test?

[cadojo]

I added the test case shown in #735 and this PR. For some reason, checks are failing for unrelated tests. Tests for derivates are now failing for 049ac92 and 8795d27, which is odd because 8795d27 tests succeeded just a few hours ago.

[ChrisRackauckas]

@YingboMa ?

[YingboMa]

This logic is too hard to understand just from the code. I don't understand what

@variables x[1:2,1:2] y[1:4]
eqs = x[:] .~ y
build_function(eqs, y, target=ModelingToolkit.CTarget())

is supposed to do.

[ChrisRackauckas]

The test failure on master got fixed though? We should update the branch when there's a passing master.

[ChrisRackauckas]

Yeah I think the issue is that the CTarget, MATLABTarget, and StanTarget are old and need some work. build_function's first argument should be the term or list of terms, not an equation. build_function(y, y, target=ModelingToolkit.CTarget()) is what it should probably be changing to. I'd be okay with accepting the patch, but I would expect this function to get overhauled.

[YingboMa]

I wouldn't want to merge this. It basically impossible to understand how +1 or some random integer is supposed to be the offset to differentiate different cases just from the code.

Also,

build_function(eqs, y, ...)

is just wrong. We shouldn't support it.

[ChrisRackauckas]

is just wrong. We shouldn't support it.

It is just wrong, these days. It's the MTK function signature from 2 years ago. It still lives on in the non-Julia targets because it hasn't been updated. And those haven't been deleted because they are still there because they are used in a bunch of tests and benchmarks. We need to completely overhaul them.

[cadojo]

Permission to close this PR in favor of the draft PR #746?

[ChrisRackauckas]

I was thinking the same thing. Yeah doing the full update of CTarget is great and I'm happy you're interested in taking it on. Thanks!