bifurcationkit / BifurcationKit.jl

A Julia package to perform Bifurcation Analysis
https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable
Other
288 stars 34 forks source link

Why do you use "lenses" to set parameters? #159

Open Datseris opened 1 week ago

Datseris commented 1 week ago

Contintuing the discussion from https://github.com/JuliaDynamics/Attractors.jl/issues/131#issuecomment-2196529224 , where @rveltz said:

Why did you use lenses at all? A beginner in Julia understands what is a mutable container and the index of the container to adjust its value. Lenses is an additional concept that so far I have not found a reason to use, instead of just having a mutable parameter container, which is also what SciML recommends for the ODE problem, and it is also what it creates if you have symbolic problems via MTK

It was a suggestion of @tkf. I think the best parameter container for bifurcation analysis is a named tuple but this is debatable. All old softwares auto, matcont, etc use index to specify the parameter but with the tooling of Julia, I think we deserve better. Look at here. You want to switch to continuation wrt U0? "hang on, what is the index". Additionally, it greatly simplifies plots where the parameter axis is correctly labelled. Even better: think about the case of 2 parameters continuation. If you dont use named tuples, you get generic p1, p2 as parameter axis.

I could add an interface like update_param!(pars, index, value) and use it all around in BK. Then it would be trivial to enable the array based parameter container (that you suggest) or the immutable one. (This is actually about to be pushed since I want to use Accessors instead of Setfield)

Why not a dict, or array then?

Because the compiler will remove the par.U0 or the change @set par.U0 = newpar so it costs nothing to be immutable (it does not allocate). I think same should be true for mutable (array) container used with with @lens _[1]. Also, if you use arrays, how do you give a name to the components?

which is also what SciML recommends for the ODE problem

where?

it is also what it creates if you have symbolic problems via MTK

Yes, I saw this. We discussed the same problem with @TorkelE but he manages to nicely propagate the bifurcation parameter in his interface. The interface MTK-BK is not fully tested yet. I like the MTK way to specify the continuation parameter but I can't expect all users to use MTK to specify their problem.

Datseris commented 1 week ago

I think the best parameter container for bifurcation analysis is a named tuple but this is debatable.

There is no such thing as "best" parameter container, as it is the user that decides what they want as a parameter container, and generic code should accommodate any parameter container. Just like there is no "best dynamical system" to do continuation for. This is just the user's input.

The real question is whether this container should be mutable (what DynamicalSystems.jl and SciML does), or immutable, like a named tuple, which is what BifurcationKit.jl seems to prefer. For immutable types you need this @lens, there is no other way around it, as you constantly need to recreate the parameter container to change any of its values. But for mutable, you don't.

This is the fundamental conceptual difference, and the main reason i my opinion to use mutable parameter containers. Immutable things, by definition, are not supposed to have their values changed. Mutable ones do. At an individual level. It makes conceptual sense to have mutable parameter containers because you typically want to change the ith parameter, whatever i is. That is why I believe from a user point of view it makes more sense to have mutable parameter containers.

Also, if you use arrays, how do you give a name to the components?

You don't. You also can't do that with @lens. If you want named containers you use either a dictionary, or a custom struct, or MTK. In dynamical systems the "parameter index" can be anything that can access the parameter container. If it is a dictionary of strings Dict("this parameter" => 3), the index can be "this parameter". If the parameter container is struct A; my_parameter; end, the the parameter index would be :my_parameter (a Symbol, which is what Julia uses when you do A.my_parameter).

The same story goes with whether your parameter container was generated by ModelingToolkit.jl. Symbolic indices should "just work". They do in DynamicalSystems.jl without requiring anyone to propagate the bifurcation parameter or in fact write any glue code.

I like the MTK way to specify the continuation parameter but I can't expect all users to use MTK to specify their problem.

You don't have to expect anything nor force anything to the user. You just need a generic interface update_param!(pars, index, value) which is exactly what DynamicalSystems.jl already has in the function set_parameter!.


which is also what SciML recommends for the ODE problem

where?

Where does SciML recommend mutable parameter containers? Explicitly, nowhere. However, I have not seen any case so far where the parameter container was immutable. In all tutorials and examples they are mutable. The better question is "where in SciML have you seen immutable parameter containers?".

Datseris commented 1 week ago

A comment on using immutable structs like named tuples for performance reasons. I don't think there is any significant performance impact. Overwtiting in-place a mutable container also doesn't cost anything. So it doesn't matter if the compiler "compiles away" the @lens call, the alternative is also costless.

A cost would come if you had to copy the parameter container for whatever reason. I am not sure whether there is such a reason, at least I haven't found one in the entirety of DynamicalSystems.jl, in continuation or beyond. If there are good reasons to copy the parameter in BFKit.jl, that is still okay, because the cost of copying the parameter container would still be overwhelmingly miniscule compared to the cost of doing literaly anythning in BFKit.jl. Computing the jacobian a single time would likely cost more, and the Jacobian is typicallyu computed millions of times.

Datseris commented 1 week ago

I see a clear advantage of @lens however: that it allows nested parameter containers, like a dictionary of dictionaries where the inner dictionaries have the parameters as values. I am not sure how I would handle this in set_parameter! in dynamical systems. The good thing is, as the code is now, one can always extend the set_parameter! function to take in a @lens as a method, if @lens is indeed the best way to handle nested parameter containers.

rveltz commented 1 week ago

The real question is whether this container should be mutable (what DynamicalSystems.jl and SciML does), or immutable,

yes. The same holds for x which for now is out of place in BK except in very sensitive parts (like collocation or trapezoid). I will discuss this during JuliaCon with Chris.

An additional requirements is to be differentiable wrt to the parameter value.

If the parameter container is struct A; my_parameter; end, the the parameter index would be :my_parameter

I like this.

rveltz commented 1 week ago

is set_parameter! differentiable wrt to the parameter value say for ForwardDiff?

Datseris commented 1 week ago

is set_parameter! differentiable wrt to the parameter value say for ForwardDiff?

I don't understand this statement. Why do you need the parameter altering to be differentiable?

As far as I know it is the dynamic rule f that needs to be differentiable w.r.t. the parameter values which is always the case. In what context you need to be able to differentiate altering the parameter values? Can this operation even be expressed mathematically? Can you provide some sort of textbook that defines what it means to "differentiate altering the value of a parameter container"?

In any case the source code of set_parameter is super simple, so maybe you can check if it satisfies this requirement:

https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/main/src/core/dynamicalsystem_interface.jl#L419-L436

all if clauses are compiled away so for non-MTK models the code is reduced to either set_index! or set_property!, both of which are likely simpler than anything @lens does.

rveltz commented 1 week ago

Sorry I should have explained the context. I am just thinking out loud to see if all cases are handled. I know what differentiable means in any case ;)

I meant i need to be able to differentiate something like this

( x,p) -> prob.F(x, update_params!( prob.params, prob.idx, p))

I need to think about this when for example using forwarddiff although I know how to do it.

Datseris commented 1 week ago

I meant i need to be able to differentiate something like this

( x,p) -> prob.F(x, update_params!( prob.params, prob.idx, p))

I admit, I do not know. My understanding is that if this works for the @lens containers, I can't see a reason why it wouldn't work for mutable containers. But it also depends on what ForwardDiff is doing.

Is it not possible to first change the parameter container,

pc # parameter container
p # parameter value
set_parameter!(pc, pindex, p)

and then evaluate the function F, or its Jacobian, at state space point x with parameter container pc? Why do these have to be at the same function? even if you want the derivative of the Jacobian, the Hessian, it would still be a function of the parameter container, H(x, pc) so you could still first alter the parameter container and then estimate the Hessian.

I guess this is probably some sort of continuation method that I don't know. Where does this need come up?

EDIT: okay, I was dead wrong, the hessian has nothing to do with this, as it is the second derivative with respect to x, not a derivative with respect to p.