Tips for making a program more Julian

Hi there! I’m new to Julia and have been experimenting with it the past few months. I had been using it more as a workspace by which I could complete my school work than as a programming language. That changed recently as one of my classes is all about simulations. The assignments in class were made with MATLAB in mind, but I wished to attempt them in Julia instead.

I was wondering if anyone might have some tips for making the code below more Julian or might be able to point out any glaring mistakes I am making with for example array manipulation or my current use of global variables.

While the code does compile and run in Julia 1.6.0 and give an output, the simulation itself seems to be incorrect, so please ignore that. :see_no_evil:

Thank you for your time!
-Vero

Code
using Markdown
using InteractiveUtils
using Plots

function midpoint_method(f::Function, x, u, t, h)
    k1, y = f(x, u, t)
    k2, _ = f(x + h/2*k1, u, t + h/2)
    k3, _ = f(x - h*k1 + 2*h*k2, u, t + h)
    x_VPG = x + h*k2
    d_VPG = h/6*(k1 - 2*k2 + k3)
    return x_VPG, y, d_VPG
end

function pt1(x, u)
    Tm = 10
    xdot = -1/Tm * x + 1/Tm * u
    y = x
    return xdot, y
end

function subtractor(u1, u2)
    return u1 - u2
end

function hysteresis(u)
    h_e = 0.085
    h_a = 0.065
    if u >= h_e
        y = 1;
    elseif u <= -h_e
        y = -1;
    else
        y = 0;
    end

    if i > 1
        if hys_state[i - 1] == 1 && u <= h_a
            y = 0;
        elseif hys_state[i - 1] == -1 && u >= -h_a
            y = 0;
        end
    end
    append!(hys_state, y)
    return y
end

function sys_PRM(x, u, t)
    _, y1::Float64 = pt1(x, 0)
    y2::Float64 = subtractor(u, y1)
    y3::Float64 = hysteresis(y2)
    xdot1, _ = pt1(x, y3)
    return xdot1, [y1 y2 y3]
end

t_step = 1;
Tm = 10
t0 = 0
tf = 20
epsilon = 1e-3
h_max = 2*Tm
h_min = 12*epsilon
u_values = []
x_values = []
y_values = Vector{Float64}()
d_values = []
t_values = []
h_values = []
hys_state = []
x = 0
h = 0.2
t = 0
i = 1
while t <= tf
    if t < t_step
        global u = 0
    else
        global u = 0.17
    end

    global x, y, d = midpoint_method(sys_PRM, x, u, t, h);
    global h = 0.9*h*min(max((epsilon/max(abs(d)))^(1/3), 0.3), 2)
    global h = min(max(h, h_min), h_max)

    append!(x_values, x)
    append!(y_values, y)
    append!(d_values, d)
    append!(h_values, h)
    append!(t_values, t)

    global t = t + h
    global i = i + 1
end

y_values = reshape(y_values, (3, :))

p1 = plot(t_values, y_values[1, :])
p2 = plot(t_values, y_values[2, :])
p3 = plot(t_values, y_values[3, :])
p4 = plot(t_values, h_values)
p5 = plot(t_values, d_values)

plot(p1, p2, p3, p4, p5, layout=(5, 1))
1 Like

Don’t. Create a function for your final while loop which builds the _values arrays.

1 Like

What would be the best way to initialize all the simulation parameters then? Leave them as globals and then pass them into the function?

I would encapsulate them in a struct and pass that to the function

1 Like

Ahh ok, that makes sense. Thank you!

What are you expecting here with y = x?

Probably this is more clean:

function pt1(x, u)
    Tm = 10
    xdot = -1/Tm * x + 1/Tm * u
    return xdot, x
end

why this subtractor function? and it could be written as

subtractor(u1,u2) = u1 - u2

These are all arrays of type Any. Use Float64[] instead.

Almost certainly these annotations of type are unnecessary.

7 Likes

Also, instead of subtractor(u1,u2) = u1 - u2, you could just write -(u1,u2) and not define a function in the first place.

4 Likes

You may take a look at this video (in French with English subtitles) showing a complete development of a small simulation project. I tried to make it as Julian as possible :wink: : Julia #1/2 : Mon premier projet Julia - YouTube

2 Likes

Yeah, I believe I was going overboard with adhering to the code style that my lecturer was using, I’ll definitely implement those changes. Thank you!

I’ll check it out, thank you!

The video maybe be a little long so, in short, I would advise to make a (properly structured) Julia project even for small projects (It takes a minute with skeletons). Doing so, you will take advantage of all the Julia’secosystem tools (Revise.jl, IDE…). I find the parameters management proposed in the video (struct with @kwdef) to be very convenient.

2 Likes

Actually I just finished the video! I liked the structure of your code, I would just need to make a few adjustments, like making an update step size function, as my simulation has to account for discontinuities in the input. I will take advantage of the package structure you used in the video, as that definitely makes the code seem orderly and easy to follow.

Quick side question, I noticed you used snake case for all function names except for getvalues, is there a particular reason for this, or is it just personal preference?

1 Like

Returning a named tuple costs nothing and may add some flexibility

function pt1(x, u)
    Tm = 10
    xdot = -1/Tm * x + 1/Tm * u
    return (; xdot, x)
end
2 Likes

Just in case it is useful, this (short and in English) video https://www.youtube.com/watch?v=n8mV7eNwqzY shows with more details how to create the project. You may find the second
long video on the same channel (adding 2D case and GPU computing) and the corresponding source GitHub - triscale-innov/ScreenCastSpring.jl useful.

Quick side question, I noticed you used snake case for all function names except for getvalues, is there a particular reason for this, or is it just personal preference?

Good catch for the inconsistent function name style :wink:
I think that the Julia conventional style is to use snake case only when the normal style gets difficult to read. Makes the convention a bit subjective. IMHO, the important part is to differentiate clearly the names of types (camel case) from the variable names.

2 Likes

For what it’s worth, YASGuide, Blue style, and the JuMP style guide all suggest using the underscores. I think only the Base Julia style guide says to omit the underscores when it’s still readable. (Linking these also since they might be useful for the original question).

7 Likes

I’ve gone ahead and tried to implement as many changes as possible in order to make my code more Julian. I’m not 100% happy with my implementation, but it definitely looks cleaner and more refined than at the start. It could be that I’m over using mutable structs, but I figured it was currently the most elegant solution.

If anyone has time, I would love for you to check out the current rendition!

PRM Github Repo

The files PRM.jl and main.jl hold the new code.

Thanks once again!
-Vero

P.S. I included a block diagram of the simulation I have attempted to program. It currently resides as block-diagram.PNG in the main directory.

This structure certainly does not need to be mutable:

    mutable struct SimOutput
        x_values::Vector{Float32} # State of PT1-Block
        u_values::Vector{Float32} # Input
        y_values::Vector{Float32} # Outputs of all blocks
        d_values::Vector{Float32} # Local Discretization Error
        t_values::Vector{Float32} # Time
        h_values::Vector{Float32} # Step size
    end
    SimOutput() = SimOutput([], [], [], [], [], [])

since all of its fields are mutable themselves.

One detail is that I think that for most here, reading these [] in the constructor is unpleasant, even if in this specific case won’t hurt because the fields are properly typed in the struct. A more Julian way to write that would be:

julia> struct SimOutput{T}
         x_values::Vector{T}
         u_values::Vector{T}
       end

julia> SimOutput(T) = SimOutput{T}(T[],T[])
SimOutput

julia> SimOutput() = SimOutput(Float32)
SimOutput

with what you can use your constructor that defaults to Float32, or can change the type
of the vectors as well:

julia> SimOutput()
SimOutput{Float32}(Float32[], Float32[])

julia> SimOutput(Float64)
SimOutput{Float64}(Float64[], Float64[])

3 Likes

The main.jl file may be placed out of the src directory (in tests or even in the top level directory).
Also return xdot, [y1, y2, y3] will allocate a new array and may lead to suboptimal perf. You may pass and mutate the modified array.

1 Like

I’m a little confused how I should go about this elegantly, as I only require the output [y1, y2, y3] once per 3 function calls, as seen here:

k1, y = sim_topology(x, u, t, state.hys_state, state.i)
k2, _ = sim_topology(x + h/2*k1, u, t + h/2, state.hys_state, state.i)
k3, _ = sim_topology(x - h*k1 + 2*h*k2, u, t + h, state.hys_state, state.i)

Would it be better to pass a boolean/enum flag instead? (I’m already not too happy passing state.hys_state and state.i as parameters, but I couldn’t think of a better syntax without everything just being state.somthing.

If you know that you have exactly 3 components, you can return a tuple (y1, y2, y3) (and have a tuple in your struct). It avoids dynamic allocations.
Edit : or you can use StaticArrays of the corresponding size (3).

1 Like