I’ll have to instrument the code again to see what’s actually happening. In this case it’s possible my starting point for k1 and k2 is poor as I’m using an explicit RK2 method for the predictor. What’s the correct way to do this?
https://diffeq.sciml.ai/stable/features/callback_library/#Manifold-Conservation-and-Projection
is our current implementation, but it’s a bit crude. We can probably do a few more Jacobian saving tricks in there. It’s a good undergrad or new contributor project so I just have left it open to improve (
). The best tome for this whole subject is Hairer’s geometric integration book which details the convergence with post-step projections, but then a bunch of schemes to avoid having to do it.
Yeah that was a GEAR thing that was preserved by LSODE and that whole series of codes, but the IRK codes have a lot more convergence rate measures that are used. But indeed, coupling adaptive time stepping with the Newton method gives wild gains, because it allows you to take a YOLO time step and then pull back if Newton isn’t converging, etc. which ultimately drops the number of calculations substantially.
Also, you might want to look at the SciMLBenchmarks, specifically:
https://benchmarks.sciml.ai/html/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.html
@SebastianM-C did a bunch of work there and showed that using really high order methods at low tolerance, and only projecting when the energy error is sufficiently high, can be one of the most efficient ways to solve such a problem. If that’s what you’re aiming for, this combination of techniques can be more efficient than a symplectic integrator.
1 Like
IIRC the Gauss-Legendre methods come with a collocation polynomial, and extrapolating that polynomial is the common predictor.
Also, note that for fully implicit Runge-Kutta methods you can exploit the structure of the A matrix by doing the Newton iteration in the eigenspace of A. Checkout Hairer II page 121 and 127 for more details.
The stage predictor for collocation RKs is the formula (8.5’) on page 120 of Hairer II.
1 Like
I mean, this isn’t really true. If you did the same thing in any other language it would be just as much of a performance issue. It’s just that most dynamic languages are so slow that it doesn’t matter and static languages don’t have any way of declaring an untyped global. So the “problem” that Julia has here is that it makes it easy to create an untyped non-constant global — because that’s a thing that people expect to work in a dynamic language — but unlike other dynamic languages, Julia is fast enough that it matters. It would be possible to fix this with some compiler work, but it’s a low priority because it’s easy to avoid once you know how, and non-constant globals are a bad idea anyway.
22 Likes
What you say makes sense. Having some background with C and Fortran as well as recent experience with Go (which conveys strong typing implicitly through a definition such as a:=1), I hadn’t thought about how some people expect to have non-typed variables in a read, evaluate, print loop. I’m also not sure dynamic typing is really needed for a REPL, as one of the first examples of an interactive language used suffixes such as $ and % to infer type.
Even though there are things that might be improved, it seems to me the most important thing for Julia right now is to stop changing. Over half the advice I find when performing a web search doesn’t work anymore on recent versions of Julia. Fortunately, this forum seems full of helpful and friendly people.
Does the const keyword essentially mean constant type rather than value?
You are lucky since Julia 1.0 was released 8 Aug 2018.
Nope, it is value (or rather the “binding” is constant, in other words, the variable will always refer to the same object).
4 Likes
You can have a global with a constant type but dynamic value by sticking it in a Ref. The only ‘downside’ is that accessing the value requires you to write a[] rather than a.
julia> const a = Ref(1)
Base.RefValue{Int64}(1)
julia> function f()
a[] = a[] + 1
a[] ^ 2
end
f (generic function with 1 method)
julia> f()
4
julia> f()
9
julia> f()
16
julia> a[]
4
This will allow the compiler to perform many optimizations always knowing that a[] isa Int, but without assuming the value is unchanging.
3 Likes
There’s an old issue to allow type annotations on global variables. This would be a great intro compiler project — not trivial by any means, but relatively straightforward. The right approach would need to be discussed with the compiler team, but I suspect the it would be to associate a type with each global binding and automatically insert a type check at each global assignment and teach the compiler that it can assume that the type of each global access has the associated type. An untyped global would then have type Any, so all globals would work the same way. Some subtleties for usability: you’ll want to at least allow redeclaring the same global with the same type so that code can be re-included interactively; you could also allow the type annotation to be made more restrictive since any code that’s already been generated with that assumption will still be correct; to be really fancy, one could add “back edges” from global binding types to methods that use those globals and increase the world age / invalidate those methods, which would allow arbitrary re-declaration of type globals at the cost of recompiling any methods that accessed the global.
17 Likes
As for using structs: Just posted an implementation in some other context.
the same code as in the link above
module Md2
Base.@kwdef mutable struct State
b1::Bool = false
b2::Bool = false
end
const st = State()
export st
end
using .Md2
julia> st.b1
false
julia> st.b1=true;
julia> st.b1
true
I my use case the performance was not relevant - but it is interesting to know how this way of encapsulation of variables into const objects compares with using Ref, array etc.
That would be an awesome feature. One of the most frustrating things in julia right now is that we don’t have a good way to define fast global parameters. There’s const, but it can’t be redefined. If we have what you suggest, we could add a keyword (or macro) typeconst which would make it an error to change the type of that variable. To define global parameters you’d just do typeconst pi = 3, then if later you want to change it to typeconst pi=3.14 you can.
as others have mentioned Ref is the solution for this, but of course this is a bit “manual” so it’d definitely be good to have both.