Local lyapunov exponents map about chaostools.jl HOT 3 CLOSED

function compute_lyapunov(prob::ODEProblem; d0=1e-9, threshold=10^4*d0, dt = 0.1,
    diff_eq_kwargs = Dict(:abstol=>d0, :reltol=>d0))

    threshold <= d0 && throw(ArgumentError("Threshold must be bigger than d0!"))

    if haskey(diff_eq_kwargs, :solver)
        solver = diff_eq_kwargs[:solver]
    else
        println("Using DPRKN12 as default solver")
        solver = DPRKN12()
    end
    # Initialize
    st1 = prob.u0
    integ1 = init(prob, solver; diff_eq_kwargs..., save_first=false, save_everystep=false)
    integ1.opts.advance_to_tstop = true
    prob.u0 = st1 .+ d0
    integ2 = init(prob, solver; diff_eq_kwargs..., save_first=false, save_everystep=false)
    integ2.opts.advance_to_tstop = true
    prob.u0 = st1
    lyapunov(integ1, integ2, prob.tspan[2]; d0=d0, threshold=threshold, dt=dt)
end

function compute_λs(q0list, p0list, tmax, d0, dt, tr, kwargs, prefix)
    n = size(q0list, 1)    # number of initial conditions
    λs = SharedArray{Float64}(n)
    pmap(i->(prob = defProb(q0list[i,:], p0list[i,:], (0., tmax));
            λs[i] = compute_lyapunov(prob, d0=d0, dt=dt, threshold=tr,
                                    diff_eq_kwargs=kwargs)),
        Progress(n),
        1:n)

    save("$prefix/lyapunov.jld", "λs", λs, "d0", d0, "dt", dt, "tr", tr,
        "tmax", tmax, "n", n);
    return λs
end

addprocs(4);
using JLD
using ArgParse
using Plots, LaTeXStrings
using ProgressMeter
using PmapProgressMeter
@everywhere begin
    using OrdinaryDiffEq, DiffEqCallbacks
    using ParallelDataTransfer
end
function main()
    # Hamiltonian parameters
    A, B, D = readdlm("param.dat")
    E_list, tmax, d0, dt, tr, solver, kwargs, defProb = input_param()
    # Broadcast parameters to all workers
    sendto(workers(), A=A, B=B, D=D, defProb=defProb)

    for E in E_list
        sendto(workers(), E=E)
        prefix = "../output/B$B D$D E$E"
        if !isdir(prefix)
            mkpath(prefix)
        end
        if isfile("$prefix/z0.jld")
            q0list, p0list = load("$prefix/z0.jld", "q0list", "p0list")
        else
            error("$prefix/z0.jld not found! Generate the initial conditions.")
        end
        λs = compute_λs(q0list, p0list, tmax, d0, dt, tr, kwargs, prefix)
    end
end