From d38963d749b73474faa80307cd6ebdb4549e4435 Mon Sep 17 00:00:00 2001 From: kmdeck Date: Tue, 6 Sep 2022 20:15:26 -0700 Subject: [PATCH 1/4] standard map wip --- examples/henon_heiles.jl | 21 ++++ examples/henon_heiles_integration_test.jl | 124 ++++++++++++++++++++++ examples/standard_map.jl | 23 ++++ examples/standard_map_integration_test.jl | 118 ++++++++++++++++++++ 4 files changed, 286 insertions(+) create mode 100644 examples/henon_heiles.jl create mode 100644 examples/henon_heiles_integration_test.jl create mode 100644 examples/standard_map.jl create mode 100644 examples/standard_map_integration_test.jl diff --git a/examples/henon_heiles.jl b/examples/henon_heiles.jl new file mode 100644 index 0000000..6dfeb0a --- /dev/null +++ b/examples/henon_heiles.jl @@ -0,0 +1,21 @@ +using Distributed +@everywhere using DifferentialEquations: ODEProblem, solve, ImplicitMidpoint + +struct HenonHeiles end + +@everywhere function deterministic_system(model::HenonHeiles) + function tendency!(dudt,u,p,t) + dudt[1] = u[3] + dudt[2] = u[4] + dudt[3] = -u[1]-2.0*u[1]*u[2] + dudt[4] = -u[2]-(u[1]^2-u[2]^2) + end + return tendency! +end + +@everywhere function evolve_deterministic_system(model::HenonHeiles, u0, tspan, dt, alg_kwargs) + tendency! = deterministic_system(model) + prob = ODEProblem(tendency!, u0, tspan, nothing) + sol = solve(prob, ImplicitMidpoint(), dt = dt, saveat = tspan[1]:dt:tspan[2]; alg_kwargs...); + return sol.u +end diff --git a/examples/henon_heiles_integration_test.jl b/examples/henon_heiles_integration_test.jl new file mode 100644 index 0000000..ee635cc --- /dev/null +++ b/examples/henon_heiles_integration_test.jl @@ -0,0 +1,124 @@ +using Revise +using RareEvents +using Statistics +using DelimitedFiles +using LinearAlgebra +using SpecialFunctions +using Plots + +examples_dir = joinpath(pkgdir(RareEvents), "examples") +include(joinpath(examples_dir, "henon_heiles.jl")) +n_processes = Sys.CPU_THREADS +addprocs(n_processes) + +# Set up model and wrapper function that evolves a single trajectory +dt = 1.0 +alg_kwargs = (); +FT = Float64 +model = HenonHeiles() +@everywhere evolve_wrapper(tspan) = (u0) -> evolve_deterministic_system(model, u0, tspan, dt, alg_kwargs) + + +# Set up rare event algorithm +τ = 1.0 +function set_ic_via_y(E, y; py = 0.0, x = 0.0) + twiceV = (x^2 + y^2 + 2 * x^2 * y - 2 / 3 * y^3) + px = sqrt(2.0 * E - py^2 - twiceV) + return [x,y, px, py] +end; +E = 0.125; +yvals = -0.35:0.001:0.35; +nensemble = length(yvals) +u0 = [set_ic_via_y(E, y) for y in yvals] +k = 0.0 +ϵ = 0.000 +metric(y) = y[2] +t_end = 1000.0 +t_start = 0.0 +tspan = (t_start, t_end) +# averaging window +T = 50 +steps_per_window = Int64.(round(T/dt)) + + +# Run GKLT - bootstrap iterations - and compute necessary statistics per iteration +iters = 5 +em = zeros(nensemble*iters) +em = reshape(em, (nensemble,iters)) +r = similar(em) +r_paper = similar(em) +σr = similar(em) +bin = 0.03 +a_range = Array(0.0:bin:0.6) +counts = zeros(length(a_range)) +algorithm_cost = nensemble*(t_end-t_start)*iters*dt +for iter in 1:iters + println(iter) + sim = RareEventSampler{Float64}(dt, τ, tspan,u0, evolve_wrapper, nensemble, metric, k, ϵ, iter) + + run!(sim); + a_m = zeros(nensemble) + lr_vector = zeros(nensemble) + λ = sum(log.(sim.weight_norm))/(t_end-t_start) + + for i in 1:nensemble + trajectory = [metric(sim.ensemble[i][j]) for j in 1:length(sim.ensemble[i])] + lr_vector[i] = likelihood_ratio(sim, trajectory, λ, t_end-t_start) + # This just counts the number of (treated as independent) block means per bin, over all trajectories, over all iterations + counts .+= count_events_per_bin(block_mean(trajectory, steps_per_window), lr_vector[i],a_range) + # This is for estimating return times via the Ragone paper method; via `return_curve`. + a_m[i] = maximum(block_mean(trajectory, steps_per_window)) + + end + em[:,iter], r[:,iter], r_paper[:,iter], σr[:, iter]= return_curve(a_m, t_end-t_start, lr_vector); + +end +rmprocs([p for p in procs() if p != myid()]) + + +# Evolve the solution directly; compute statistics +tspan_direct = (0.0, 3.e6) +direct_solution = evolve_deterministic_system(model, u0[end], tspan_direct, dt, (maxiters = 1e8,)); +direct_cost = (tspan_direct[2]-tspan_direct[1])*dt +u = [metric(direct_solution[k]) for k in 1:length(direct_solution)] +A = block_mean(u, steps_per_window) +m = 100 +em_direct, r_direct, r_paper_direct, σr_direct = return_curve(A, FT(T), FT.(ones(length(A)))) + +# Fit GEV +m = 100 +blocks = block_maxima(A, m) +params, nll = fit_gev(blocks, [std(blocks), mean(blocks), 0.1]) +cdf_A= (gev_cdf.(a_range, params[1], params[2], params[3])).^(1.0/m) +pdf_A= (cdf_A[2:end] .- cdf_A[1:end-1])/bin + + +# Make plots +plot1 = plot() +mean_a = zeros(length(a_range)-1) +mean_rtn = zeros(length(a_range)-1) +std_rtn = zeros(length(a_range)-1) +for i in 1:(length(a_range)-1) + mask = (em[:] .< a_range[i+1]) .& ( em[:] .>= a_range[i]) + if sum(mask) >0 + mean_a[i] = mean(em[:][mask]) + mean_rtn[i] = log10.(mean(r_paper[:][mask])) + std_rtn[i] = std(r_paper[:][mask])/mean(r_paper[:][mask])./log(10.0) + end +end +nonzero = (mean_a .!= 0.0) .& (isnan.(std_rtn) .== 0) +final_a = mean_a[nonzero] +final_r = mean_rtn[nonzero] +final_σr = std_rtn[nonzero] +plot!(plot1,final_a, final_r, ribbon = final_σr, label = "GKLT Algorithm; using maxes") +denominator = 1.0 .- cumsum(counts)/sum(counts) +safe_denominator = denominator .> 10.0*eps(FT) +plot!(plot1,a_range[safe_denominator], log10.(50.0 ./denominator[safe_denominator]), label = "GKLT Algorithm; using counts") +plot!(plot1,em_direct, log10.(r_paper_direct), ribbon = σr_direct ./ r_paper_direct ./ log(10.0), label = "Direct Integration") +plot!(plot1, a_range, log10.(50.0 ./ (1.0 .- cdf_A)), label = "GEV") +plot!(plot1, a_range, log10.(50.0 ./ (1.0 .- 0.5 .*(1 .+erf.(a_range ./2 .^0.5 ./std(A))))), label = "Gaussian") +plot!(plot1,xlabel = "Event magnitude") +plot!(plot1,ylabel = "Log10(Return Time)") +plot!(plot1, legend = :topleft) +plot!(plot1, xrange = [0.0,1.0]) +savefig(plot1, "return_curve_ou.png") diff --git a/examples/standard_map.jl b/examples/standard_map.jl new file mode 100644 index 0000000..04ab479 --- /dev/null +++ b/examples/standard_map.jl @@ -0,0 +1,23 @@ +using Distributed + +struct StandardMap{FT<: AbstractFloat} + K::FT +end + +@everywhere function mapping_system(model::StandardMap) + function update!(u) + u[2] = mod(u[2] + model.K*sin(u[1]),2π) + u[1] = mod(u[1] + u[2],2π) + end + return update! +end + +@everywhere function evolve_mapping_system(model::StandardMap, u0, tspan) + update! = mapping_system(model) + nsteps = Int(tspan[2]-tspan[1]) + solution = map(0:nsteps) do (i) + update!(u0) + copy(u0) + end + return solution +end diff --git a/examples/standard_map_integration_test.jl b/examples/standard_map_integration_test.jl new file mode 100644 index 0000000..b297aed --- /dev/null +++ b/examples/standard_map_integration_test.jl @@ -0,0 +1,118 @@ +using Revise +using RareEvents +using Statistics +using DelimitedFiles +using LinearAlgebra +using SpecialFunctions +using Plots + +examples_dir = joinpath(pkgdir(RareEvents), "examples") +include(joinpath(examples_dir, "standard_map.jl")) +n_processes = Sys.CPU_THREADS +addprocs(n_processes) + +# Set up model and wrapper function that evolves a single trajectory +dt = 1.0 +alg_kwargs = (); +FT = Float64 +model = StandardMap{Float64}(4.0) +@everywhere evolve_wrapper(tspan) = (u0) -> evolve_mapping_system(model, u0, tspan) + + +# Set up rare event algorithm +τ = 10.0 +yvals = 0.0:0.1:2π; +nensemble = length(yvals) +u0 = [[y, 0.0] for y in yvals] +k = -0.2 +ϵ = 0.001 +metric(y) = abs(y[1]-π) +t_end = 500.0 +t_start = 0.0 +tspan = (t_start, t_end) +# averaging window +T = 20 +steps_per_window = Int64.(round(T/dt)) + + +# Run GKLT - bootstrap iterations - and compute necessary statistics per iteration +iters = 1 +em = zeros(nensemble*iters) +em = reshape(em, (nensemble,iters)) +r = similar(em) +r_paper = similar(em) +σr = similar(em) +bin = 0.03 +a_range = Array(π:bin:2π) +counts = zeros(length(a_range)) +algorithm_cost = nensemble*(t_end-t_start)*iters*dt +for iter in 1:iters + println(iter) + sim = RareEventSampler{Float64}(dt, τ, tspan,u0, evolve_wrapper, nensemble, metric, k, ϵ, iter) + + run!(sim); + a_m = zeros(nensemble) + lr_vector = zeros(nensemble) + λ = sum(log.(sim.weight_norm))/(t_end-t_start) + + for i in 1:nensemble + trajectory = [metric(sim.ensemble[i][j]) for j in 1:length(sim.ensemble[i])] + lr_vector[i] = likelihood_ratio(sim, trajectory, λ, t_end-t_start) + # This just counts the number of (treated as independent) block means per bin, over all trajectories, over all iterations + counts .+= count_events_per_bin(block_mean(trajectory, steps_per_window), lr_vector[i],a_range) + # This is for estimating return times via the Ragone paper method; via `return_curve`. + a_m[i] = maximum(block_mean(trajectory, steps_per_window)) + + end + em[:,iter], r[:,iter], r_paper[:,iter], σr[:, iter]= return_curve(a_m, t_end-t_start, lr_vector); + +end +rmprocs([p for p in procs() if p != myid()]) + + +# Evolve the solution directly; compute statistics +tspan_direct = (0.0, algorithm_cost) +direct_solution = evolve_mapping_system(model, u0[end-1], tspan_direct) +direct_cost = (tspan_direct[2]-tspan_direct[1])*dt +u = [metric(direct_solution[k]) for k in 1:length(direct_solution)] +A = block_mean(u, steps_per_window) +m = 100 +em_direct, r_direct, r_paper_direct, σr_direct = return_curve(A, FT(T), FT.(ones(length(A)))) + +# Fit GEV +m = 100 +blocks = block_maxima(A, m) +params, nll = fit_gev(blocks, [std(blocks), mean(blocks), 0.1]) +cdf_A= (gev_cdf.(a_range, params[1], params[2], params[3])).^(1.0/m) +pdf_A= (cdf_A[2:end] .- cdf_A[1:end-1])/bin + + +# Make plots +plot1 = plot() +mean_a = zeros(length(a_range)-1) +mean_rtn = zeros(length(a_range)-1) +std_rtn = zeros(length(a_range)-1) +for i in 1:(length(a_range)-1) + mask = (em[:] .< a_range[i+1]) .& ( em[:] .>= a_range[i]) + if sum(mask) >0 + mean_a[i] = mean(em[:][mask]) + mean_rtn[i] = log10.(mean(r_paper[:][mask])) + std_rtn[i] = std(r_paper[:][mask])/mean(r_paper[:][mask])./log(10.0) + end +end +nonzero = (mean_a .!= 0.0) .& (isnan.(std_rtn) .== 0) +final_a = mean_a[nonzero] +final_r = mean_rtn[nonzero] +final_σr = std_rtn[nonzero] +plot!(plot1,final_a, final_r, ribbon = final_σr, label = "GKLT Algorithm; using maxes") +denominator = 1.0 .- cumsum(counts)/sum(counts) +safe_denominator = denominator .> 10.0*eps(FT) +plot!(plot1,a_range[safe_denominator], log10.(50.0 ./denominator[safe_denominator]), label = "GKLT Algorithm; using counts") +plot!(plot1,em_direct, log10.(r_paper_direct), ribbon = σr_direct ./ r_paper_direct ./ log(10.0), label = "Direct Integration") +plot!(plot1, a_range, log10.(50.0 ./ (1.0 .- cdf_A)), label = "GEV") +plot!(plot1, a_range, log10.(50.0 ./ (1.0 .- 0.5 .*(1 .+erf.(a_range ./2 .^0.5 ./std(A))))), label = "Gaussian") +plot!(plot1,xlabel = "Event magnitude") +plot!(plot1,ylabel = "Log10(Return Time)") +plot!(plot1, legend = :topleft) +plot!(plot1, xrange = [0.0,1.0]) +savefig(plot1, "return_curve_ou.png") From d87023850e398faf6529b06bcfd154daf4e5daa2 Mon Sep 17 00:00:00 2001 From: kmdeck Date: Tue, 6 Sep 2022 20:39:21 -0700 Subject: [PATCH 2/4] finding regular orbits --- examples/standard_map_integration_test.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/standard_map_integration_test.jl b/examples/standard_map_integration_test.jl index b297aed..e8c0c0d 100644 --- a/examples/standard_map_integration_test.jl +++ b/examples/standard_map_integration_test.jl @@ -8,8 +8,8 @@ using Plots examples_dir = joinpath(pkgdir(RareEvents), "examples") include(joinpath(examples_dir, "standard_map.jl")) -n_processes = Sys.CPU_THREADS -addprocs(n_processes) +#n_processes = Sys.CPU_THREADS +#addprocs(n_processes) # Set up model and wrapper function that evolves a single trajectory dt = 1.0 @@ -24,10 +24,10 @@ model = StandardMap{Float64}(4.0) yvals = 0.0:0.1:2π; nensemble = length(yvals) u0 = [[y, 0.0] for y in yvals] -k = -0.2 +k = -0.5 ϵ = 0.001 -metric(y) = abs(y[1]-π) -t_end = 500.0 +metric(y) = abs(y[2]-π) + abs(y[1] -π) +t_end = 300.0 t_start = 0.0 tspan = (t_start, t_end) # averaging window @@ -67,7 +67,7 @@ for iter in 1:iters em[:,iter], r[:,iter], r_paper[:,iter], σr[:, iter]= return_curve(a_m, t_end-t_start, lr_vector); end -rmprocs([p for p in procs() if p != myid()]) +#rmprocs([p for p in procs() if p != myid()]) # Evolve the solution directly; compute statistics From bdf0dded7c9edd4db0819c2baf654eee908b6989 Mon Sep 17 00:00:00 2001 From: kmdeck Date: Wed, 7 Sep 2022 07:59:28 -0700 Subject: [PATCH 3/4] lyapunov_weighted_Dynamics --- examples/standard_map.jl | 12 +-- examples/standard_map_integration_test.jl | 122 +++++----------------- 2 files changed, 33 insertions(+), 101 deletions(-) diff --git a/examples/standard_map.jl b/examples/standard_map.jl index 04ab479..137db50 100644 --- a/examples/standard_map.jl +++ b/examples/standard_map.jl @@ -1,18 +1,18 @@ -using Distributed - struct StandardMap{FT<: AbstractFloat} K::FT end -@everywhere function mapping_system(model::StandardMap) +function mapping_system(model::StandardMap) function update!(u) - u[2] = mod(u[2] + model.K*sin(u[1]),2π) - u[1] = mod(u[1] + u[2],2π) + u[4] = u[4] + model.K*cos(u[1])*u[3] # δI + u[3] = u[3] + u[4] # δθ + u[2] = mod(u[2] + model.K*sin(u[1]),2π) # I + u[1] = mod(u[1] + u[2],2π) # θ end return update! end -@everywhere function evolve_mapping_system(model::StandardMap, u0, tspan) +function evolve_mapping_system(model::StandardMap, u0, tspan) update! = mapping_system(model) nsteps = Int(tspan[2]-tspan[1]) solution = map(0:nsteps) do (i) diff --git a/examples/standard_map_integration_test.jl b/examples/standard_map_integration_test.jl index e8c0c0d..c3e6e62 100644 --- a/examples/standard_map_integration_test.jl +++ b/examples/standard_map_integration_test.jl @@ -2,117 +2,49 @@ using Revise using RareEvents using Statistics using DelimitedFiles -using LinearAlgebra -using SpecialFunctions using Plots examples_dir = joinpath(pkgdir(RareEvents), "examples") include(joinpath(examples_dir, "standard_map.jl")) -#n_processes = Sys.CPU_THREADS -#addprocs(n_processes) -# Set up model and wrapper function that evolves a single trajectory dt = 1.0 alg_kwargs = (); FT = Float64 -model = StandardMap{Float64}(4.0) -@everywhere evolve_wrapper(tspan) = (u0) -> evolve_mapping_system(model, u0, tspan) +K_sm =2.0 +model = StandardMap{Float64}(K_sm) +evolve_wrapper(tspan) = (u0) -> evolve_mapping_system(model, u0, tspan) # Set up rare event algorithm τ = 10.0 -yvals = 0.0:0.1:2π; -nensemble = length(yvals) -u0 = [[y, 0.0] for y in yvals] -k = -0.5 -ϵ = 0.001 -metric(y) = abs(y[2]-π) + abs(y[1] -π) -t_end = 300.0 +nensemble = 150 +u0 = [[rand()*2π, rand()*2π, 1e-5, 1e-5] for _ in 1:nensemble] +kvals = [0.0,-0.05,-0.1] +ϵ = 0.0001 +metric(y; K= K_sm) = (y[4]*y[3]*2.0*(1.0+K*cos(y[1])))/(y[3]^2.0+y[4]^2.0) +t_end = 150.0 t_start = 0.0 tspan = (t_start, t_end) -# averaging window -T = 20 -steps_per_window = Int64.(round(T/dt)) - +plot1 = plot() +extract_I(u) = u[2] +extract_θ(u) = u[1] +colors = [:cyan, :purple, :red] +ms = [1.0, 1.0, 1.0] +labels = ["k=0", "k = -0.05", "k = -0.1"] +# Run GKLT +for j in 1:2 + k = kvals[j] -# Run GKLT - bootstrap iterations - and compute necessary statistics per iteration -iters = 1 -em = zeros(nensemble*iters) -em = reshape(em, (nensemble,iters)) -r = similar(em) -r_paper = similar(em) -σr = similar(em) -bin = 0.03 -a_range = Array(π:bin:2π) -counts = zeros(length(a_range)) -algorithm_cost = nensemble*(t_end-t_start)*iters*dt -for iter in 1:iters - println(iter) sim = RareEventSampler{Float64}(dt, τ, tspan,u0, evolve_wrapper, nensemble, metric, k, ϵ, iter) - run!(sim); - a_m = zeros(nensemble) - lr_vector = zeros(nensemble) - λ = sum(log.(sim.weight_norm))/(t_end-t_start) - - for i in 1:nensemble - trajectory = [metric(sim.ensemble[i][j]) for j in 1:length(sim.ensemble[i])] - lr_vector[i] = likelihood_ratio(sim, trajectory, λ, t_end-t_start) - # This just counts the number of (treated as independent) block means per bin, over all trajectories, over all iterations - counts .+= count_events_per_bin(block_mean(trajectory, steps_per_window), lr_vector[i],a_range) - # This is for estimating return times via the Ragone paper method; via `return_curve`. - a_m[i] = maximum(block_mean(trajectory, steps_per_window)) - - end - em[:,iter], r[:,iter], r_paper[:,iter], σr[:, iter]= return_curve(a_m, t_end-t_start, lr_vector); - -end -#rmprocs([p for p in procs() if p != myid()]) - - -# Evolve the solution directly; compute statistics -tspan_direct = (0.0, algorithm_cost) -direct_solution = evolve_mapping_system(model, u0[end-1], tspan_direct) -direct_cost = (tspan_direct[2]-tspan_direct[1])*dt -u = [metric(direct_solution[k]) for k in 1:length(direct_solution)] -A = block_mean(u, steps_per_window) -m = 100 -em_direct, r_direct, r_paper_direct, σr_direct = return_curve(A, FT(T), FT.(ones(length(A)))) - -# Fit GEV -m = 100 -blocks = block_maxima(A, m) -params, nll = fit_gev(blocks, [std(blocks), mean(blocks), 0.1]) -cdf_A= (gev_cdf.(a_range, params[1], params[2], params[3])).^(1.0/m) -pdf_A= (cdf_A[2:end] .- cdf_A[1:end-1])/bin - - -# Make plots -plot1 = plot() -mean_a = zeros(length(a_range)-1) -mean_rtn = zeros(length(a_range)-1) -std_rtn = zeros(length(a_range)-1) -for i in 1:(length(a_range)-1) - mask = (em[:] .< a_range[i+1]) .& ( em[:] .>= a_range[i]) - if sum(mask) >0 - mean_a[i] = mean(em[:][mask]) - mean_rtn[i] = log10.(mean(r_paper[:][mask])) - std_rtn[i] = std(r_paper[:][mask])/mean(r_paper[:][mask])./log(10.0) + flatten_θ = extract_θ.(sim.ensemble[1]) + flatten_I = extract_I.(sim.ensemble[1]) + for i in 2:nensemble + flatten_θ = vcat(flatten_θ, extract_θ.(sim.ensemble[i])) + flatten_I = vcat(flatten_I, extract_I.(sim.ensemble[i])) end + scatter!(plot1, flatten_θ, flatten_I, color = colors[j], ms = ms[j], label = labels[j], markerstrokewidth=0) + end -nonzero = (mean_a .!= 0.0) .& (isnan.(std_rtn) .== 0) -final_a = mean_a[nonzero] -final_r = mean_rtn[nonzero] -final_σr = std_rtn[nonzero] -plot!(plot1,final_a, final_r, ribbon = final_σr, label = "GKLT Algorithm; using maxes") -denominator = 1.0 .- cumsum(counts)/sum(counts) -safe_denominator = denominator .> 10.0*eps(FT) -plot!(plot1,a_range[safe_denominator], log10.(50.0 ./denominator[safe_denominator]), label = "GKLT Algorithm; using counts") -plot!(plot1,em_direct, log10.(r_paper_direct), ribbon = σr_direct ./ r_paper_direct ./ log(10.0), label = "Direct Integration") -plot!(plot1, a_range, log10.(50.0 ./ (1.0 .- cdf_A)), label = "GEV") -plot!(plot1, a_range, log10.(50.0 ./ (1.0 .- 0.5 .*(1 .+erf.(a_range ./2 .^0.5 ./std(A))))), label = "Gaussian") -plot!(plot1,xlabel = "Event magnitude") -plot!(plot1,ylabel = "Log10(Return Time)") -plot!(plot1, legend = :topleft) -plot!(plot1, xrange = [0.0,1.0]) -savefig(plot1, "return_curve_ou.png") +plot!(plot1, xlabel = "θ", ylabel = "I", xticks = ([0,π,2π],["0","π","2π"]), yticks = ([0, π, 2π],["0","π","2π"])) +savefig(plot1, "standard_map.png") From 2967a1d25cfa6aab6f3f90289d3cfeef701d0a3f Mon Sep 17 00:00:00 2001 From: kmdeck Date: Wed, 31 May 2023 08:36:29 -0700 Subject: [PATCH 4/4] adding old file --- examples/standard_map_integration_test.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/examples/standard_map_integration_test.jl b/examples/standard_map_integration_test.jl index c3e6e62..99598ec 100644 --- a/examples/standard_map_integration_test.jl +++ b/examples/standard_map_integration_test.jl @@ -17,7 +17,7 @@ evolve_wrapper(tspan) = (u0) -> evolve_mapping_system(model, u0, tspan) # Set up rare event algorithm τ = 10.0 -nensemble = 150 +nensemble = 300 u0 = [[rand()*2π, rand()*2π, 1e-5, 1e-5] for _ in 1:nensemble] kvals = [0.0,-0.05,-0.1] ϵ = 0.0001 @@ -26,6 +26,7 @@ t_end = 150.0 t_start = 0.0 tspan = (t_start, t_end) plot1 = plot() +plot2 = plot() extract_I(u) = u[2] extract_θ(u) = u[1] colors = [:cyan, :purple, :red] @@ -35,16 +36,20 @@ labels = ["k=0", "k = -0.05", "k = -0.1"] for j in 1:2 k = kvals[j] - sim = RareEventSampler{Float64}(dt, τ, tspan,u0, evolve_wrapper, nensemble, metric, k, ϵ, iter) + sim = RareEventSampler{Float64}(dt, τ, tspan,u0, evolve_wrapper, nensemble, metric, k, ϵ, j) run!(sim); flatten_θ = extract_θ.(sim.ensemble[1]) flatten_I = extract_I.(sim.ensemble[1]) + flatten_λ = -sum(metric.(sim.ensemble[1])*dt)/(tspan[2]-tspan[1]) for i in 2:nensemble flatten_θ = vcat(flatten_θ, extract_θ.(sim.ensemble[i])) flatten_I = vcat(flatten_I, extract_I.(sim.ensemble[i])) + flatten_λ = vcat(flatten_λ, sum(metric.(sim.ensemble[i])*dt)/(tspan[2]-tspan[1])) end scatter!(plot1, flatten_θ, flatten_I, color = colors[j], ms = ms[j], label = labels[j], markerstrokewidth=0) - + plot!(plot2, flatten_λ, seriestype=:stephist, label = labels[j], bins = 10) end plot!(plot1, xlabel = "θ", ylabel = "I", xticks = ([0,π,2π],["0","π","2π"]), yticks = ([0, π, 2π],["0","π","2π"])) +plot!(plot2, xlabel = "Finite time estimate of Lyapunov exponent λ", ylabel = "Relative Frequency") savefig(plot1, "standard_map.png") +savefig(plot2, "lyapunov_histograms.png")