Function depending on the global variable inside module

Hello,
I want to define a global variable inside a module. This global variable can be set by the user and the behaviour (the returned value) of the functions inside the module depends on this global variable.
For example:

module M
  global AMP
  set_amp(x) = (AMP = x)
  func1(x; amp) = amp * x
  func1(x) = func1(x; amp=AMP)
end

I want to achieve something like that

M.set_amp(5.0)
M.func1(2.0)
#  the return value is 10.0

However, I got the error instead
ERROR: UndefVarError: AMP not defined
Stacktrace:
[1] func1(x::Float64)
@ Main.M ~/test.jl:5
[2] top-level scope
@ REPL[3]:1

Is there any approach to do the similar thing? What is the recommended way? Thanks!

Function has local scope of variables. So you can make it realize global variables via

function set_amp(x)
    global AMP = x
end
4 Likes

Thanks for the answer. Does this approach hinder the performance? I am wondering if there is other better approach to control the behavior of functions inside a module.

Yes. See the performance tips: Performance Tips · The Julia Language

Yes. Pass data as parameters.

It’s not just better from a performance perspective, but also from a code-organization perspective — having your functions depend on global parameters means that it’s harder to tell what your code is doing because it could be affected by code someplace else that changed the parameter.

6 Likes

You can just define the function with a default parameter

func1(x; amp=DEFAULT_VALUE)

and pass a different value to it when needed

2 Likes

Thanks for the answer!

I want to develop a physics simulation package. Every physical quantity is normalized by a user-chosen characteristic length, says lambda.

module M
  func1(t; lambda) = do_something
  func2(t; lambda) = do_something
  func3(t; lambda) = do_something
  func4(t; lambda) = do_something
end

func1, func2, func3, func4 are part of simulation procedure. Since there will be only one characteristic length in the entire simulation, I am trying to provide a way to set the characteristic length before the simulation, instead of passing it every time when calling the simulation procedure. Would it be possible to achieve that?

Thanks for the answer!
Not sure if this approach has the similar performance issue mentioned by @stevengj ?

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

14 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