Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Create function for expanding e.g. mm expressions, and use for e.g. HomotopyContinuation #713

Merged
merged 14 commits into from
Dec 2, 2023

Conversation

TorkelE
Copy link
Member

@TorkelE TorkelE commented Nov 7, 2023

This should prevent us from having to use the complicated Symbolics version we used previously. Hopefully, we can also use it under simpler circumstances later.

Tests added, if everything passes, and it looks good, it should be gtg.

One thing we can consider, is whether we should make this option default in latexify when the form=:ODE option is used. I think it is probably prettier to have these expanded.

@isaacsas
Copy link
Member

isaacsas commented Nov 7, 2023

Can we not just drop all the registration and derivatives we hard code and let symbolics pick it up by tracing, i.e. for mm we just have

mm(X, v, K) = v * X / (X + K)

and nothing else.

@TorkelE
Copy link
Member Author

TorkelE commented Nov 7, 2023

I guess, but that would be a separate issue, right?
(I am not really sure of the implications of doing so, so cannot really say much about it)

@isaacsas
Copy link
Member

isaacsas commented Nov 7, 2023

No, that is the simple way to handle this instead of adding this expansion kwarg. I don't think we should be post-processing all symbolic expressions and substituting out the registered function version for the explicit expression. This would have to include rate and affect functions for jumps, all the observed equations, symbolic expressions within events, etc, and seems like a nightmare to have to maintain going forward as we add stuff within ReactionSystems.

Instead we should just get rid of the function registration, which would then cause the function to be symbolically evaluated immediately when first used in the DSL or in a reaction. This is essentially our previous definition. Or we should leave alone and just modify the Latex printing as you suggest so it doesn't show the functions. (Could we modify how Latexify prints the functions instead, i.e. adding a custom rule for what mm should be printed as via Latexify? That might be even better since then the MTK systems we generate would print nicely too.)

@isaacsas
Copy link
Member

isaacsas commented Nov 7, 2023

Maybe we can do this as a transformation on a ReactionSystem instead. I think that will be easier. i.e.

rn_expanded = expand_rate_functions(rn)

It would still need to be applied to lots of stuff beyond just the Reactions and Equations but seems more manageable to maintain then trying to essentially do it for every possible system type we can convert to.

This is also more inline with things like expand_connections in MTK.

@TorkelE
Copy link
Member Author

TorkelE commented Nov 7, 2023

Wouldn't that require us to, when we e.g. use homotopy continuation, create a new ReactionSystem, which we then convert to a ODESystem? Just to make clear (which it maybe was already). This code doesn't change any default behaviours, it just permits me to, when I create a NonlinearSystem to use for HomotopyContinuation, tell the conversion function that in this case I want to substitute in the mm functions for the full expression. When we talked last time my understanding was that this is what you wanted (rather than the current reregister function)

If we want to change how e.g. mm function are used, as you suggest, I cannot really say much, but sounds like it could be useful.

In either case, one thing we should also consider, is that currently we export mm to the scope where Catalyst is used. Especially this one is annoying since it e.g. conflicts with the mm in Plots.Measures. If there is a solution which does not do this (and which replaces in Catalyst.mm instead) that might be good. This might even be easier if we do what you suggest.

@TorkelE
Copy link
Member Author

TorkelE commented Nov 8, 2023

Maybe we can discuss next time we meet? I think I understand how you want it set up, but if I'm to implement it I'd want to be thoroughly sure how you'd want it set up before I start working.

@TorkelE
Copy link
Member Author

TorkelE commented Nov 9, 2023

The 1.6 CI have a weird thing:
image

Otherwise this should be all good an an improvement on what we currently have.

@TorkelE
Copy link
Member Author

TorkelE commented Nov 14, 2023

This one is ready for merging as well

@isaacsas
Copy link
Member

Ok, will take another look and get back to you.

ns = convert(NonlinearSystem, rs; remove_conserved = true)
pre_varmap = [symmap_to_varmap(rs,u0)..., symmap_to_varmap(rs,ps)...]
conservationlaw_errorcheck(rs, pre_varmap)
p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, parameters(ns); defaults = ModelingToolkit.defaults(ns))
p_dict = Dict(parameters(ns) .=> p_vals)
eqs_pars_funcs = vcat(equations(ns), conservedequations(rs))
eqs_funcs = map(eq -> substitute(eq.rhs - eq.lhs, p_dict), eqs_pars_funcs)
eqs = [deregister([mm, mmr, hill, hillr, hillar], eq) for eq in eqs_funcs]
eqs = map(eq -> substitute(eq.rhs - eq.lhs, p_dict), eqs_pars_funcs)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just curious, why did you add this now?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

didn't actually add anything. The

eqs = [deregister([mm, mmr, hill, hillr, hillar], eq) for eq in eqs_funcs]

step was deleted. but the variable eqs was passed to the next stage. Hence the variable eqs_funcs was renamed to eqs directly.

Copy link
Member

@isaacsas isaacsas left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add a test with non-reaction equations to make sure they work too (given you added that functionality)?

@TorkelE
Copy link
Member Author

TorkelE commented Nov 16, 2023

Yes, I can add that, that makes sense.

@isaacsas
Copy link
Member

Do we have tests on the expansion directly? I.e. that it generates the expected symbolic equations? (Sorry, not at my computer to check easily.) it would be good to add a few of those too.

Once you’ve got the extra tests settled I’m happy for you to merge this.

@TorkelE
Copy link
Member Author

TorkelE commented Nov 16, 2023

I've added these tests, hopefully they should cover all cases:

# Tests `ReactionSystem`s.
let
    @variables t 
    @species x(t) y(t)
    @parameters k v n 
    rs1 = @reaction_network rs begin
        mm(x, v, k), 0 --> x
        mmr(x, v, k), 0 --> x
        hill(x, v, k, n), 0 --> x
        hillr(x, v, k, n), 0 --> x
        hillar(x, y, v, k, n), 0 --> x    
    end
    rs2 = @reaction_network rs begin
        v * x / (x + k), 0 --> x
        v * k / (x + k), 0 --> x
        v * (x^n) / (x^n + k^n), 0 --> x
        v * (k^n) / (x^n + k^n), 0 --> x
        v * (x^n) / (x^n + y^n + k^n), 0 --> x    
    end

    @test Catalyst.expand_registered_functions(rs1) == rs2
end

# Tests `Reaction`s.
let
    @variables t 
    @species x(t) y(t)
    @parameters k v n 
    
    r1 = @reaction mm(x, v, k), 0 --> x
    r2 = @reaction mmr(x, v, k), 0 --> x
    r3 = @reaction hill(x, v, k, n), 0 --> x
    r4 = @reaction hillr(x, v, k, n), 0 --> x
    r5 = @reaction hillar(x, y, v, k, n), 0 --> x + y
    
    @test isequal(Catalyst.expand_registered_functions(r1).rate, v * x / (x + k))
    @test isequal(Catalyst.expand_registered_functions(r2).rate, v * k / (x + k))
    @test isequal(Catalyst.expand_registered_functions(r3).rate, v * (x^n) / (x^n + k^n))
    @test isequal(Catalyst.expand_registered_functions(r4).rate, v * (k^n) / (x^n + k^n))
    @test isequal(Catalyst.expand_registered_functions(r5).rate, v * (x^n) / (x^n + y^n + k^n))
end

# Tests `Equation`s.
let
    @variables T X(T) Y(T)
    @parameters K V N 
    
    eq1 = 0 ~ mm(X, V, K)
    eq2 = 0 ~ mmr(X, V, K)
    eq3 = 0 ~ hill(X, V, K, N)
    eq4 = 0 ~ hillr(X, V, K, N)
    eq5 = 0 ~ hillar(X, Y, V, K, N)
    
    @test isequal(Catalyst.expand_registered_functions(eq1), 0 ~ V * X / (X + K))
    @test isequal(Catalyst.expand_registered_functions(eq2), 0 ~ V * K / (X + K))
    @test isequal(Catalyst.expand_registered_functions(eq3), 0 ~ V * (X^N) / (X^N + K^N))
    @test isequal(Catalyst.expand_registered_functions(eq4), 0 ~ V * (K^N) / (X^N + K^N))
    @test isequal(Catalyst.expand_registered_functions(eq5), 0 ~ V * (X^N) / (X^N + Y^N + K^N))
end

@isaacsas
Copy link
Member

isaacsas commented Dec 2, 2023

@TorkelE are you done with this? It seems like this is holding up a bunch of your other PRs right? (I had said above you could merge once you added the tests, so feel free too!)

@@ -17,6 +17,7 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
Requires = "ae029012-a4dd-5104-9daa-d747884805df"
RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, this needs a Compat

@TorkelE
Copy link
Member Author

TorkelE commented Dec 2, 2023

compat update works, merging now

@TorkelE TorkelE merged commit b179191 into master Dec 2, 2023
8 of 9 checks passed
@TorkelE TorkelE deleted the custom_function_expansion branch June 8, 2024 18:30
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants