Function depending on the global variable inside module

Make a Simulation struct that contains lambda and have functions operate on that object.

6 Likes

The performance would be similar as long as you get rid of global variables.

1 Like

I created the Jupyter and Pluto notebooks with the same contents as below:

Very Long!

Problem-Algorithm-Solver pattern

A code will be more readable and efficient if the global variables that store the parameters of a problem are passed as arguments to functions everytime.

However, many people find it troublesome to pass the parameters as arguments to the function everytime.

In order to resolve this misunderstanding, they can try the problem-algorithm-solver pattern explained below.

In order to use the (; a, b, c) = p syntax, require VERSION ≥ v"1.7.0-beta" #39285.

VERSION ≥ v"1.7.0-beta" #> true

A minimal working example of the problem-algorithm-solver pattern:

module FreeFall

"""Problem type"""
Base.@kwdef struct Problem{G, Y0, V0, TS}
    g::G = 9.80665
    y0::Y0 = 0.0
    v0::V0 = 30.0
    tspan::TS = (0.0, 8.0)
end

"""Algorithm types"""
Base.@kwdef struct EulerMethod{T}  dt::T = 0.1 end
Base.@kwdef struct ExactFormula{T} dt::T = 0.1 end
default_algorithm(::Problem) = EulerMethod()

"""Solution type"""
struct Solution{Y, V, T, P<:Problem, A} y::Y; v::V; t::T; prob::P; alg::A end

"""Solver function"""
solve(prob::Problem) = solve(prob, default_algorithm(prob))

function solve(prob::Problem, alg::ExactFormula)
    (; g, y0, v0, tspan) = prob
    (; dt) = alg
    t0, t1 = tspan
    t = range(t0, t1 + dt/2; step = dt)    
    y(t) = y0 + v0*(t - t0) - g*(t - t0)^2/2
    v(t) = v0 - g*(t - t0)
    Solution(y.(t), v.(t), t, prob, alg)
end

function solve(prob::Problem, alg::EulerMethod)
    (; g, y0, v0, tspan) = prob
    (; dt) = alg
    t0, t1 = tspan
    t = range(t0, t1 + dt/2; step = dt)    
    n = length(t)    
    y = Vector{typeof(y0)}(undef, n)
    v = Vector{typeof(v0)}(undef, n)
    y[1] = y0
    v[1] = v0
    for i in 1:n-1
        v[i+1] = v[i] - g*dt
        y[i+1] = y[i] + v[i]*dt
    end
    Solution(y, v, t, prob, alg)
end

end

This module will calculate the free-fall motion under a uniform gravity potential.

  • Problem: The parameters of problems are stored in objects of type FreeFall.Problem.

  • Algorithm: The parameters of algorithms are stored in objects of types FreeFall.EulerMethod and of FreeFall.ExactFormula.

  • Solver function: The function FreeFall.solve(prob::FreeFall.Problem, alg) returns the object of type FreeFall.Solution obtained by solving the problem prob with the algorithm alg.

Note that the (; a, b, c) = p syntax of Julia ≥ v1.7, is very useful for extracting the contents of the object that contains parameters. (You can also do the same thing with @unpack in the exellent well-known package Parameters.jl.)

Example 1: Compare a numerical solution and an exact one.

using Plots

earth = FreeFall.Problem()
sol_euler = FreeFall.solve(earth)
sol_exact = FreeFall.solve(earth, FreeFall.ExactFormula())

plot(sol_euler.t, sol_euler.y; label="Euler's method (dt = $(sol_euler.alg.dt))", ls=:auto)
plot!(sol_exact.t, sol_exact.y; label="exact solution", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft)

Since the time step dt = 0.1 is not enough small, the error of the Euler method is rather large.

Example 2: Change the algorithm parameter dt smaller.

earth = FreeFall.Problem()
sol_euler = FreeFall.solve(earth, FreeFall.EulerMethod(dt = 0.01))
sol_exact = FreeFall.solve(earth, FreeFall.ExactFormula())

plot(sol_euler.t, sol_euler.y; label="Euler's method (dt = $(sol_euler.alg.dt))", ls=:auto)
plot!(sol_exact.t, sol_exact.y; label="exact solution", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft)

Example 3: Change the problem parameter g to the gravitational acceleration on the moon.

moon = FreeFall.Problem(g = 1.62, tspan = (0.0, 40.0))
sol_euler = FreeFall.solve(moon)
sol_exact = FreeFall.solve(moon, FreeFall.ExactFormula(dt = sol_euler.alg.dt))

plot(sol_euler.t, sol_euler.y; label="Euler's method (dt = $(sol_euler.alg.dt))", ls=:auto)
plot!(sol_exact.t, sol_exact.y; label="exact solution", ls=:auto)
title!("On the Moon"; xlabel="t", legend=:bottomleft)

Extensibility of the Julia ecosystem: A different author of the FreeFall module can define , in the SomeExtension module, a new algorithm type and a new method of the FreeFall.solve function. Both a new type and a new method!

module SomeExtension

using ..FreeFall
using ..FreeFall: Problem, Solution

Base.@kwdef struct Symplectic2ndOrder{T}  dt::T = 0.1 end

function FreeFall.solve(prob::Problem, alg::Symplectic2ndOrder)
    (; g, y0, v0, tspan) = prob
    (; dt) = alg
    t0, t1 = tspan
    t = range(t0, t1 + dt/2; step = dt)    
    n = length(t)    
    y = Vector{typeof(y0)}(undef, n)
    v = Vector{typeof(v0)}(undef, n)
    y[1] = y0
    v[1] = v0
    for i in 1:n-1
        ytmp = y[i] + v[i]*dt/2
        v[i+1] = v[i] - g*dt
        y[i+1] = ytmp + v[i+1]*dt/2
    end
    Solution(y, v, t, prob, alg)
end

end

earth = FreeFall.Problem()
sol_sympl = FreeFall.solve(earth, SomeExtension.Symplectic2ndOrder(dt = 2.0))
sol_exact = FreeFall.solve(earth, FreeFall.ExactFormula())

plot(sol_sympl.t, sol_sympl.y; label="2nd order symplectic (dt = $(sol_sympl.alg.dt))", ls=:auto)
plot!(sol_exact.t, sol_exact.y; label="exact solution", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft)

The second-order symplectic method gives exact discretizations of the free-fall motions under uniform gravity potentials.

Composability of the Julia ecosystem: By combining the MonteCarloMeasurements.jl package with the modules given above, we can plot the behavior of the system for perturbations of the initial value. We did not intend for them to be used in this way!

using MonteCarloMeasurements

earth = FreeFall.Problem(y0 = 0.0 ± 0.0, v0 = 30.0 ± 1.0)
sol_euler = FreeFall.solve(earth)
sol_sympl = FreeFall.solve(earth, SomeExtension.Symplectic2ndOrder(dt = 2.0))
sol_exact = FreeFall.solve(earth, FreeFall.ExactFormula())

ylim = (-100, 60)
P = plot(sol_euler.t, sol_euler.y; label="Euler's method (dt = $(sol_euler.alg.dt))", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft, ylim)
Q = plot(sol_sympl.t, sol_sympl.y; label="2nd order symplectic (dt = $(sol_sympl.alg.dt))", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft, ylim)
R = plot(sol_exact.t, sol_exact.y; label="exact solution", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft, ylim)
plot(P, Q, R; size=(720, 600))

Edit 1: Typo. PLotsPlots
Edit 2: Add Jupyter notebook
Edit 3: Add Pluto notebooks and Julia v1.6 versions

15 Likes

Globals make your code unusable in multithreading. If someone tries to solve a physics problem with your package in parallel for different lengths and you use globals… It will Bork heavily. Use the suggested approaches above and it will be possible to multithread!

4 Likes

Thanks for your answer. Although my case is easier than the example you provided, I learn a lot from your example.

What is the purpose of using the 4 different type parameters (G, Y0, V0, TS) in the declaration of Problem structure?

struct Problem{G, Y0, V0, TS}
    g::G = 9.80665
    y0::Y0 = 0.0
    v0::V0 = 30.0
    tspan::TS = (0.0, 8.0)
end

Why not just using one type parameter instead

struct Problem{T}
    g::T = 9.80665
    y0::T = 0.0
    v0::T = 30.0
    tspan::Tuple{T, T} = (0.0, 8.0)
end
1 Like

I created the Jupyter notebook with this change:

Most results are the same as the original. However, only the last example results in an error:

MethodError: no method matching Main.FreeFall.Problem(::Float64, ::Particles{Float64, 2000}, ::Particles{Float64, 2000}, ::Tuple{Float64, Float64})
Closest candidates are:
  Main.FreeFall.Problem(::T, ::T, ::T, ::Tuple{T, T}) where T at In[2]:13

In the last example of using MonteCarloMeasurements, the types of the fields y0 and v0 are set to MonteCarloMeasurements.Particles{Float64, 2000} instead of Float64.

Fixed version:

Base.@kwdef struct Problem{T, U}
    g::T = 9.80665
    y0::U = 0.0
    v0::U = 30.0
    tspan::Tuple{T, T} = (0.0, 8.0)
end

Thus, overly restricting the type combinations may inhibit the composability, the strengths of the Julia ecosystem.

When in doubt, it is safe to mechanically separate the type parameters of all potentially different types of fields.

Doing so has the advantage of reducing the amount of something to worry about and makes it easier to concentrate on mathematics and physics. :wink:

2 Likes

@genkuroki, thanks for this outstanding post, a great service to promote Julia language IMHO.

As an exercise, added to your code the different physical units using Unitful and UnitfulRecipes. It worked everywhere at the first attempt except for the MonteCarloMeasurements plots. Do you know of any recipe for that?

NB: will be glad to post the preliminary edited code if useful.

By saying multithreading, do you mean things like openMP?

However, if use MPI instead, will global variables still have such problems you mentioned?

Is it always recommended to parameterize the type of a struct? For instance, if a struct has 10 fields, then do we need to specify 10 type parameters for each field?

Sorry for a stupid question, what does the … mean in

using ..FreeFall

I know if I write a module called FreeFall.jl, I can do

include("FreeFall.jl")
using .FreeFall

But what does the … do in using …FreeFall?

Thanks!

Yes.

First, for type stability. (cf. Avoid fields with abstract type)

Second, for composability with various packages that cannot be predicted in advance. (Example: Composability with MonteCarloMeasurements.jl)

If you really don’t want to write type parameters, you can use ConcreteStructs.jl. The macro ConcreteStructs.@concrete automatically add type parameters to all fields in struct.

See [ANN] ConcreteStructs.jl - Cut the boilerplate when concretely parameterizing structs

Small working example. (Note that we have an issue someone should solve. Play nicely with other struct def macros · Issue #4 · jonniedie/ConcreteStructs.jl · GitHub)

using ConcreteStructs

Base.@kwdef @concrete struct Problem g; y0; v0; tspan end
Problem(;
    g = 9.80665,
    y0 = 0.0,
    v0 = 30.0,
    tspan = (0.0, 8.0)
) = Problem(g, y0, v0, tspan)

prob1 = Problem()

Output: Problem{Float64, Float64, Float64, Tuple{Float64, Float64}}(9.80665, 0.0, 30.0, (0.0, 8.0))

prob2 = Problem(v0 = 100.0)

Output: Problem{Float64, Float64, Float64, Tuple{Float64, Float64}}(9.80665, 0.0, 100.0, (0.0, 8.0))

using MonteCarloMeasurements
prob3 = Problem(y0 = 0.0 ± 0.0, v0 = 30.0 ± 1.0)

Output: Problem{Float64, Particles{Float64, 2000}, Particles{Float64, 2000}, Tuple{Float64, Float64}}(9.80665, 0.0, 30.0 ± 1.0, (0.0, 8.0))

1 Like

Often you don’t want to, and instead want to parametrize many of your vars with the same parameter. For example,

Base.@kwdef struct Problem{T}
    g::T = 9.80665
    y0::T = 0.0
    v0::T= 30.0
    tspan::Tuple{T,T} = (0.0, 8.0)
end
2 Likes

Also, only if you need them to be parametric. They can also have static types:

struct A
  x::Int 
  y::Float64
end

in which case you don’t need any parameter whatsoever.

2 Likes

no I mean things like

@threads for i in 1:n .... or @spawn begin .... end

Yes, that is basically “like OpenMP:” shared-memory threads, with decorator directives to make it easy. As opposed to distributed-memory programming like MPI.

Thanks. The second point above seems to be a design choice. Sometimes I feel that it is quite cumbersome to specify the argument type in a function if that type itself is a parametric one.

Consider again your above example. In solve(prob::Problem, alg::EulerMethod), no constraints are exerted on Problem. However, I usually have some constraints in mind during design. For example, I will not expect that the type parameter G to be a String (at least in my design). Then, shall I enforce this constraint in the signature of solve(prob::Problem, alg::EulerMethod) or simply throw an error during execution?

Thank you!
Yes, I was thinking it must be something related with shared memory stuff that can cause issue for global variable. For MPI, it seems using global variable is safe because each core are independent.

I think both options can be fine. Sometimes I prefer to throw an error on execution, because the error message can be more user friendly.

2 Likes

Thanks @lmiq .

Did you mean we check the type on the first line of the function body (e.g., using something like @assert) or just let it run and encounter finally an error spontaneously? :grinning_face_with_smiling_eyes:

Yes, one or other option. The use of an assertion allow you to throw a customized error message.

(personally I would like to have an option to throw an error without launching the complete stack trace, meaning, I would like to have something that behaves like this:

julia> function f(x)
         if typeof(x) != Int
           println("x must be an Int")
           return nothing
         end
         x
       end
f (generic function with 1 method)

julia> f(1.0)
x must be an Int

Instead of:

julia> function f(x)
         @assert typeof(x) == Int "x must be an Int"
         x
       end
f (generic function with 1 method)

julia> f(1.0)
ERROR: AssertionError: x must be an Int
Stacktrace:
 [1] f(x::Float64)
   @ Main ./REPL[6]:2
 [2] top-level scope
   @ REPL[7]:1

In this case the stack trace is not large, but in some cases it can be overwhelming for a user.

1 Like