JumpSystems for constant and variable rate jumps

[isaacsas]

@ChrisRackauckas a few issues here to address before proceeding:

1. This needs a version of build_function that generates a function only taking an integrator (for affects as we discussed).
2. I'm planning that convert will be valid for producing a JumpSet unless you have another suggestion. Will just directly create JumpProblems from JumpSystems and bypass JumpSets.
3. MassActionJumps only allow for rates that are numbers so that is also problematic. I set up the conversion code to error if a user tries to convert a MassActionJump with a rate that has a type of Operation or Variable. (Since Variable is itself a Number I couldn't check if the type was a Number...) Otherwise we have to do some revising of DiffEqJump to somehow extract rates from parameter vectors (which means the mass action jumps must know their index or key in the parameter data structure). I'm actually not quite sure how to best handle this in converting a ReactionSystem over... Plan for 2. will get around this...

[isaacsas]

For 3., I guess we can just create a new JumpProblem constructor here that requires a parameter map like the new ODEProblem constructor. That would allows us to generate MassActionJumps that have real numbers for rates at that point instead of in the conversion phase.

[ChrisRackauckas]

it might be easier to do this the other way around, i.e.

function generate_affect_function(js, affect)
    f = build_function(...)
    function (integrator)
       u = integrator.u
       f(u,u,p,t)
    end
end

kind of thing.

[isaacsas]

Do we know for sure that the compiler would optimize out the extra function call?

Another option would be to modify build_function to create a different header for the generated functions based on a kwarg (say integratorarg=false by default). If true we use something like (for oop):

ex_head = :(integrator -> begin 
($(fargs.args...),) = (integrator.u,integrator.p,integrator.t)
end

with a corresponding header for the normal (u,p,t) and (du,u,p,t) cases extracted from what is currently there. Then we merge these headers with everything in the current begin to end block within ex_iip and ex_oop.

Or I could have build_function generate code, intercept it, and modify the header on the generated expressions. Then we could call mk_function within here.

[ChrisRackauckas]

Or I could have build_function generate code, intercept it, and modify the header on the generated expressions. Then we could call mk_function within here.

That'll be hard to support in the program staging, since it would then need to separately call to GeneralizedGenerated.

Do we know for sure that the compiler would optimize out the extra function call?

Thinking about it more, no, I don't think it would inline well, so we should add the header support.

[isaacsas]

Ok, I have a first pass at adding a header. It seems that generated constant and variable rate jumps are consistent with hand specified rate/affect functions. Mass action jumps and problem creation are still to do.

[isaacsas]

@ChrisRackauckas So it looks like build_function doesn't resolve from the states and eqs which unknowns the equations that are passed in correspond to? I'm finding for the two-reaction SIR example in the docs that the affect! function that is generated for the I -> R reaction sets S and I instead of I and R. Below is what is generated for S+I -> 2I and I -> R. Notice the indices used on the in-place output vector in both cases correspond to S and I. Is there a way to specify which unknowns the equations correspond to updating?

(function = (##MTKIntegrator#433;) -> begin
    begin
        (var"##MTIIPVar#430", var"##MTKArg#426", var"##MTKArg#427", var"##MTKArg#428") = ((var"##MTKIntegrator#433").u, (var"##MTKIntegrator#433").u, (var"##MTKIntegrator#433").p, (var"##MTKIntegrator#433").t)
        begin
            $(Expr(:inbounds, true))
            local var"#80#val" = begin
                        let (S, I, R, β, γ, t) = (var"##MTKArg#426"[1], var"##MTKArg#426"[2], var"##MTKArg#426"[3], var"##MTKArg#427"[1], var"##MTKArg#427"[2], var"##MTKArg#428")
                            begin
                                var"##MTIIPVar#430"[1] = (S = (ModelingToolkit).:-(S, 1))
                                var"##MTIIPVar#430"[2] = (I = (ModelingToolkit).:+(I, 1))
                            end
                        end
                    end
            $(Expr(:inbounds, :((ModelingToolkit).pop)))
            var"#80#val"
        end
        (ModelingToolkit).nothing
    end
end, function = (##MTKIntegrator#446;) -> begin
    begin
        (var"##MTIIPVar#443", var"##MTKArg#439", var"##MTKArg#440", var"##MTKArg#441") = ((var"##MTKIntegrator#446").u, (var"##MTKIntegrator#446").u, (var"##MTKIntegrator#446").p, (var"##MTKIntegrator#446").t)
        begin
            $(Expr(:inbounds, true))
            local var"#84#val" = begin
                        let (S, I, R, β, γ, t) = (var"##MTKArg#439"[1], var"##MTKArg#439"[2], var"##MTKArg#439"[3], var"##MTKArg#440"[1], var"##MTKArg#440"[2], var"##MTKArg#441")
                            begin
                                var"##MTIIPVar#443"[1] = (I = (ModelingToolkit).:-(I, 1))
                                var"##MTIIPVar#443"[2] = (R = (ModelingToolkit).:+(R, 1))
                            end
                        end
                    end
            $(Expr(:inbounds, :((ModelingToolkit).pop)))
            var"#84#val"
        end
        (ModelingToolkit).nothing
    end
end)

[ChrisRackauckas]

Is there a way to specify which unknowns the equations correspond to updating?

We should maybe add an output_idxs argument. This could be done by views, but we could avoid the possible performance issue by just generating better code, so we might as well.

[isaacsas]

Yeah, I was just hacking up something to that effect...

…

On Sun, May 3, 2020 at 7:05 PM Christopher Rackauckas < ***@***.***> wrote: Is there a way to specify which unknowns the equations correspond to updating? We should maybe add an output_idxs argument. This could be done by views, but we could avoid the possible performance issue by just generating better code, so we might as well. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#317 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACHTJT7CHI635P5PPIBWKN3RPX2FBANCNFSM4MLDNYBA> .

[isaacsas]

@ChrisRackauckas With this jumps seem to be working now. If this design looks ok to you I can add tests, and if you want add the oidx call for the other equation types of rhss I can add them too (right now it is only applied for the default rhss case).

[ChrisRackauckas]

probably fine to just have it here.

[isaacsas]

Ok, this should be good to go. I'll plan to add MassActionJumps next, and then can work on documenting this all. Hopefully next weekend.