Solving difference equation: Part 2

In my previous post I was trying to get to reach a conclusion whether I can get more speed in solving a system of difference equation (basically a discrete dynamical system) through a hand-written code. I got useful suggestions, but today I found that the trajectory() function in DynamicalSystems runs faster. My ultimate goal was to write a code that solves the system in a biparametric space (varying two parameters in nested for loop) and store the corresponding period in a file. Here is the system and its jacobian

using DynamicalSystems, DelimitedFiles, .Threads,  BenchmarkTools
## Components of a test DiscreteDynamicalSystem

function dds_constructor(u0 = [0.5, 0.7];  r=1.0, k=2.0)
    return DiscreteDynamicalSystem(dds_rule, u0, [r, k], dds_jac)
end
## equations of motion:
function dds_rule(x, par, n)
   r, k = par;
   a = 5.0; mu = 0.5; d = 0.2;
   dx = x[1]*exp(r*(1-x[1]/k)-x[2]/(a+x[1]^2))
   dy = x[2]*exp(mu*x[1]/(a+x[1]^2)-d)
   return @SVector [dx, dy]
end
## Jacobian:
function dds_jac( x, par, n)
    r, k = par;
    a = 5.0; mu = 0.5; d = 0.2;
    J11 = exp(- x[2]/(x[1]^2 + a) - r*(x[1]/k - 1)) - x[1]*exp(- x[2]/(x[1]^2 + a) - r*(x[1]/k - 1))*(r/k - (2*x[1]*x[2])/(x[1]^2 + a)^2)
    J12 = -(x[1]*exp(- x[2]/(x[1]^2 + a) - r*(x[1]/k - 1)))/(x[1]^2 + a)
    J21 = x[2]*exp((mu*x[1])/(x[1]^2 + a) - d)*(mu/(x[1]^2 + a) - (2*mu*x[1]^2)/(x[1]^2 + a)^2)
    J22 = exp((mu*x[1])/(x[1]^2 + a) - d)
    return @SMatrix [J11  J12; J21 J22]
end

Some functions I used:

function seqper_new(x; tol=1e-3)       # function to calculate periodicity 
    @inbounds for k in 2:length(x)
        if abs(x[k] - x[1]) ≤ tol
            all(j -> abs(x[j] - x[j-k+1]) ≤ tol, k:lastindex(x)) && return k - 1
        end
    end 
    return length(x)
end
function meshh(x,y)    # create a meshgrid of x and y
    len_x = length(x);
    len_y = length(y);
    xmesh = Float64[]
    ymesh = Float64[]
    for i in 1:len_x
        for j in 1:len_y
            push!(xmesh, x[i])
            push!(ymesh, y[j])
        end
    end
    return  xmesh, ymesh 
end

Here is my function that I want to optimize. It calculates the period of each pair of xpts*ypts parameter values and store the corresponding x-values in a column of some preallocated variable.

function isoperiodic_test(ds, par_area::Matrix{Float64}, NIter::Int64, filename::String,
                         nxblock::Int64, nyblock::Int64, NTr::Int64, xpts::Int64, ypts::Int64)    

    io = open(iso_filename, "a")   # opeing file to write output
    p1_st =  par_area[1];          # beginning of first parameter 
    p1_nd =  par_area[2];          # end of first parameter
    p2_st =  par_area[3];          # beginning of second parameter
    p2_nd =  par_area[4];          # end of second parameter

    p1_block = range(p1_st, p1_nd, length=nxblock + 1);    # divide par_area in blocks 
    p2_block = range(p2_st, p2_nd, length=nyblock + 1);    # divide par_area in blocks 
    l_bipar = xpts * ypts;                                 # total number of pts in each block
    total_blocks = nxblock * nyblock                       # total number of blocks
    sol_last = Array{Float64,2}(undef, NIter + 1, l_bipar);   # stores x values of solution for l_bipar (r,k) paris

    for ii in 1:nxblock
         step1size = (p1_block[ii + 1] - p1_block[ii]) / xpts;
         par1range = range(p1_block[ii], p1_block[ii + 1] - step1size, length=xpts);
        for jj in 1:nyblock
            display("Progress: " * string(nxblock * (ii - 1) + jj) * " out of " * string(total_blocks) * " steps.")
            step2size = (p2_block[jj + 1] - p2_block[jj]) / ypts;
            par2range = range(p2_block[jj], p2_block[jj + 1] - step2size, length=ypts);   
            Threads.@threads for i = 1:xpts         
                for j = 1:ypts
                     set_parameter!(ds, [par1range[i], par2range[j]]) # change the parameter values
                     tr = trajectory(ds, NIter; Ttr=NTr);             # find the solution
                     sol_last[:,(i - 1) * xpts + j] =  tr[:,1];       # store x-values
                end
            end
            periods = seqper_new.(eachcol(sol_last), tol=0.001)       # seqper_new calculates periodicity 
            par1mesh, par2mesh = meshh(par1range, par2range)          # meshgrid of parameter values
            writedlm(io, [par1mesh par2mesh periods])                 # write data in txt file
        end 
    end
    close(io)
    return iso_filename
end
##  Test
 
parameter_area = [1.0 5.0 2.0 5.0]
nxblock = 1;
nyblock = 1;
xpts = 100;
ypts = 100;
NIter = 2000;
NTr = 50000;
init = [0.4, 0.5];
ds = dds_constructor(init);
##
iso_filename =  @time isoperiodic_test(ds, parameter_area, NIter, "isoperiodic_test_data.txt",
nxblock, nyblock, NTr, xpts, ypts)

Observation:

"Progress: 1 out of 1 steps."
  4.281922 seconds (100.33 k allocations: 628.321 MiB, 17.13% gc time)
"isoperiodic_test_data.txt_1.txt"

for 2 × 2 blocks

"Progress: 1 out of 4 steps."
"Progress: 2 out of 4 steps."
"Progress: 3 out of 4 steps."
"Progress: 4 out of 4 steps."
 16.914337 seconds (401.27 k allocations: 2.007 GiB, 10.98% gc time)
"isoperiodic_test_data.txt_1.txt"

Am I on the right track? I am still learning. Trying my best to find answers by myself but often find it difficult optimizing my own codes. Any form of help is appreciated.

If you find this difficult to read/understand let me know.

There are a few bugs in your code.

  • The code for k in 2:length(x) in the function seqper_new(x; tol=1e-3) must be corrected as for k in 2:(length(x) ÷ 2 + 1).
  • nxblock * (ii - 1) + jj must be nyblock * (ii - 1) + jj.
  • (i - 1) * xpts + j must be (i - 1) * ypts + j
  • The Threads.@threads for loop in the function isoperiodic_test is not thread-safe, so that the result changes with each execution. Threads are conflicting at the line set_parameter!(ds, [par1range[i], par2range[j]]), where the same ds will be modified by multiple threads.

Memory allocations are also unnecessarily too much.

  • The size of the array sol_last in memory is (NIter + 1) * xpts * ypts * 8 bytes = 2001 * 100 * 100 * 8 bytes = 153 MB, which is quite large, but it is discarded at the end of the function isoperiodic_test.
  • Each execution of the DynamicalSystems.trajectory function results in a wasted memory allocation, which is repeated tens of thousands of times.

Other points that need improvement:

  • The calculation results are saved to a file, and the results cannot be reused without going through the file.
  • The code is unnecessarily complicated.
  • You have created the complicated function isoperiodic_test to do most of all the work you want to do.

Furthermore, it is not necessary to use DynamicalSystems.jl in this case.

My test implementation:

using Base.Threads
using Plots
pyplot() # https://github.com/JuliaPlots/Plots.jl/issues/3560
using CSV
using DataFrames
function f(x, p)
   r, k = p
   a, mu, d = 5.0, 0.5, 0.2
   xnew1 = x[1]*exp(r*(1-x[1]/k)-x[2]/(a+x[1]^2))
   xnew2 = x[2]*exp(mu*x[1]/(a+x[1]^2)-d)
   xnew1, xnew2
end

function traj2d!(tr1, tr2, f, p, x0, niters, nwarmup)
    x = x0
    for _ in 1:nwarmup
        x = f(x, p)
    end
    @inbounds tr1[1], tr2[1] = x
    for i in 2:niters+1
        x = f(x, p)
        @inbounds tr1[i], tr2[i] = x
    end
end
function findperiod(x; tol=1e-3)
    n = length(x)
    @inbounds for k in 2:(n ÷ 2 + 1)
        if abs(x[k] - x[1]) ≤ tol
            all(j -> abs(x[j] - x[j-k+1]) ≤ tol, k:n) && return k - 1
        end
    end 
    return n
end

The function that calculates the period corresponding to the parameters (It does not save the result to a file):

function calcperiod2d(f = f, niters = 2000, nwarmup = 50000, x0 = (0.4, 0.5),
        r = (1:0.05:5)[1:end-1], k = (2:0.05:5)[1:end-1], tol=1e-3
    )
    tmp1 = [Vector{Float64}(undef, niters+1) for _ in 1:nthreads()]
    tmp2 = [Vector{Float64}(undef, niters+1) for _ in 1:nthreads()]
    period = Matrix{Int}(undef, length(r), length(k))
    @threads for j in eachindex(k)
        tid = threadid()
        tr1, tr2 = tmp1[tid], tmp2[tid]
        for i in eachindex(r)
            p = (r[i], k[j])
            traj2d!(tr1, tr2, f, p, x0, niters, nwarmup)
            @inbounds period[i, j] = findperiod(tr1; tol)
        end
    end
    period, r, k
end

Plotting function (I want to see the results immediately by plotting.):

xm1ox(x) = (x - 1)/x # x minus 1 over x : {x ≥ 1} -> {0 ≤ y < 1}
inv1my(y) = 1/(1 - y) # inverse of 1 minus 1 : {0 ≤ y < 1} -> {x ≥ 1}

function plot_period(r, k, period; kwargs...)
    cbtick = [1, 2, 3, 4, 5, 6, 10, 30]
    colorbar_ticks = (xm1ox.(cbtick), string.(cbtick))
    heatmap(r, k, xm1ox.(period)'; xlabel="r", ylabel="k", title="period",
        colorbar_ticks, kwargs...)
end

Save and load functions (I rarely use them, but many people may want to keep the results in a file.):

function save_period(r, k, period; fn = "result.csv")
    ii, jj = reim(complex.(axes(r, 1), axes(k, 1)'))
    rr, kk = reim(complex.(r, k'))
    df = DataFrame(i = vec(ii), j = vec(jj), r = vec(rr), k = vec(kk), period = vec(period))
    CSV.write(fn, df)
end

function load_period(fn = "result.csv")
    df = CSV.read(fn, DataFrame)
    imax = maximum(df.i)
    jmax = maximum(df.j)
    r = reshape(df.r, imax, jmax)[1:imax]
    k = reshape(df.k, imax, jmax)[1:imax:end]
    period = reshape(df.period, imax, jmax)
    period, r, k
end

Example:

period, r, k = @time calcperiod2d()
period, r, k = @time calcperiod2d()
period, r, k = @time calcperiod2d()
plot_period(r, k, period)
  3.109461 seconds (517.41 k allocations: 33.035 MiB, 7.54% compilation time)
  2.870107 seconds (228 allocations: 433.594 KiB)
  2.909849 seconds (221 allocations: 433.531 KiB)

result

Calculate, save, load, and plot the result:

period, r, k = calcperiod2d()
fn = save_period(r, k, period; fn = "result.csv")
period, r, k = load_period(fn)
plot_period(r, k, period)

result

Jupyter notebook: public/On Solving difference equation Part 2.ipynb at main · genkuroki/public · GitHub

I love mathematics and undoubtedly, so do the members of this community.

If you find something mathematically interesting, please let us know. :blush:

1 Like

Thank you very much @genkuroki , you are great.

  • The reason I was saving the data is to plot them in MATLAB (by calling a user-defined function from Julia) as I am not very good at plotting them in Julia.
  • I used DynamicalSystems.jl so that I can modify the code a little bit to also calculate the Lyapunov spectrum for all those parameter values and compare the iso-periodic diagram with the Lyapunov one. One way to do it is using lyapunovspectrum from DynamicalSystems.jl library, otherwise I have to write the code for finding lypunov spectrum myself.

What should I do? Should I open an issue in DynamicalSystems.jl to add another function trajectory! that takes a preallocated trajectory data and outputs the modified one. Will that be useful?

What if I want to use DynamicalSystems.jl along with @threads, can I make it threadsafe by defining the dynamical system inside the for loop where I put the new parameter values?

julia> plot_period(r, k, period)
Error showing value of type Plots.Plot{Plots.PyPlotBackend}:
ERROR: MethodError: no method matching getproperty(::Char, ::String)
Closest candidates are:
  getproperty(::DataFrames.GroupKey, ::AbstractString) at C:\Users\Hermesr\.julia\packages\DataFrames\vuMM8\src\groupeddataframe\groupeddataframe.jl:551
  getproperty(::DataFrames.DataFrameColumns, ::AbstractString) at C:\Users\Hermesr\.julia\packages\DataFrames\vuMM8\src\abstractdataframe\iteration.jl:212
  getproperty(::DataFrameRow, ::AbstractString) at C:\Users\Hermesr\.julia\packages\DataFrames\vuMM8\src\dataframerow\dataframerow.jl:300
  ...
Stacktrace:
  [1] py_set_axis_colors(sp::Plots.Subplot{Plots.PyPlotBackend}, ax::PyCall.PyObject, a::Plots.Axis)
    @ Plots C:\Users\Hermesr\.julia\packages\Plots\hyS17\src\backends\pyplot.jl:815
  [2] _before_layout_calcs(plt::Plots.Plot{Plots.PyPlotBackend})
    @ Plots C:\Users\Hermesr\.julia\packages\Plots\hyS17\src\backends\pyplot.jl:1100
  [3] prepare_output(plt::Plots.Plot{Plots.PyPlotBackend})
    @ Plots C:\Users\Hermesr\.julia\packages\Plots\hyS17\src\plot.jl:177
  [4] display(d::VSCodeServer.InlineDisplay, m::MIME{Symbol("image/svg+xml")}, x::Plots.Plot{Plots.PyPlotBackend})
    @ VSCodeServer c:\Users\Hermesr\.vscode\extensions\julialang.language-julia-1.3.32\scripts\packages\VSCodeServer\src\display.jl:0
  [5] display(d::VSCodeServer.InlineDisplay, mime::String, x::Any)
    @ Base.Multimedia .\multimedia.jl:216
  [6] display(d::VSCodeServer.InlineDisplay, x::Plots.Plot{Plots.PyPlotBackend})
    @ VSCodeServer c:\Users\Hermesr\.vscode\extensions\julialang.language-julia-1.3.32\scripts\packages\VSCodeServer\src\display.jl:109
  [7] display(x::Any)
    @ Base.Multimedia .\multimedia.jl:328
  [8] #invokelatest#2
    @ .\essentials.jl:708 [inlined]
  [9] invokelatest
    @ .\essentials.jl:706 [inlined]
 [10] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\REPL.jl:247
 [11] (::REPL.var"#40#41"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\REPL.jl:231
 [12] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\REPL.jl:462
 [13] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\REPL.jl:229
 [14] (::REPL.var"#do_respond#61"{Bool, Bool, REPL.var"#72#82"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\REPL.jl:798
 [15] #invokelatest#2
    @ .\essentials.jl:708 [inlined]
 [16] invokelatest
    @ .\essentials.jl:706 [inlined]
 [17] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\LineEdit.jl:2441
 [18] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\REPL.jl:1126
 [19] (::REPL.var"#44#49"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL .\task.jl:411

installed version of matplotlib

Z:\>pip show matplotlib
Name: matplotlib
Version: 3.4.1
Summary: Python plotting package
Home-page: https://matplotlib.org
Author: John D. Hunter, Michael Droettboom
Author-email: matplotlib-users@python.org
License: PSF

Same as Use Plots() library to graph a colored surface - #8 by genkuroki