From 389531abe334572b93561c9d710751e63c4f6479 Mon Sep 17 00:00:00 2001 From: mschauer Date: Thu, 30 Jan 2020 23:21:52 +0100 Subject: [PATCH 01/29] SIR model --- scripts/inference/sir.jl | 116 ++++++++++++++++++++++++++++ src/BridgeSDEInference.jl | 1 + src/BridgeSDEInference_for_tests.jl | 2 + src/examples/sir.jl | 87 +++++++++++++++++++++ 4 files changed, 206 insertions(+) create mode 100644 scripts/inference/sir.jl create mode 100644 src/examples/sir.jl diff --git a/scripts/inference/sir.jl b/scripts/inference/sir.jl new file mode 100644 index 0000000..1212516 --- /dev/null +++ b/scripts/inference/sir.jl @@ -0,0 +1,116 @@ +SRC_DIR = joinpath(Base.source_dir(), "..", "..", "src") +OUT_DIR = joinpath(Base.source_dir(), "..", "..", "output") +mkpath(OUT_DIR) +using Makie: lines, lines! + +include(joinpath(SRC_DIR, "BridgeSDEInference_for_tests.jl")) + +using StaticArrays +using Distributions +using Random +# Let's generate the data +# ----------------------- +using Bridge +#import BridgeSDEInference: CIR, CIRaux +DIR = "auxiliary" +include(joinpath(SRC_DIR, DIR, "data_simulation_fns.jl")) +include(joinpath(SRC_DIR, DIR, "utility_functions.jl")) +Random.seed!(4) +pop = 50_000_000 +K = 1.0 +θˣ = [0.37, 0.05, 0.1, 0.01, K] +Pˣ = SIR(θˣ...) + +x0, dt, T = ℝ{2}(1/pop, 0.), 1/5000, 30.0 +tt = 0.0:dt:T + +XX, _ = simulate_segment(ℝ{2}(1.0, 0.0), x0, Pˣ, tt) +last(XX)[2]*pop + +#lines(XX.tt, K .- sum.(XX.yy)) +lines(XX.tt, pop*first.(XX.yy)) +lines!(XX.tt, pop*last.(XX.yy)) + +θ_init = copy(θˣ) +Pˣ = SIR(θ_init...) + +length(XX.tt) +skip = 20000 + +Σdiagel = 10/pop +Σ = SMatrix{2,2}(1.0I)*Σdiagel +L = @SMatrix[1.0 0.0; 0.0 1.0] + +obs_time, obs_vals = XX.tt[1:skip:end], [rand(Gaussian(L*x, Σ)) for x in XX.yy[1:skip:end]] + +days = [0.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0] +cases = [1.0 0.0; 62.0 0.0; 121.0 0.1; 198.0 3.0; 291.0 6.0; 440.0 9.0; 571.0 17.0; 830.0 25.0; 1287.0 41.0; 1975.0 56.0; 2744.0 80.0; 4515.0 106.0; 5974.0 132.0; 7771.0 170.0] +if true + obs_time = days + obs_vals = (1/pop)*reinterpret(SVector{2,Float64}, (cases)') +end + +P̃ = [SIRAux(θ_init..., t₀, u, T, v) for (t₀, T, u, v) + in zip(obs_time[1:end-1], obs_time[2:end], obs_vals[1:end-1], obs_vals[2:end])] + +model_setup = DiffusionSetup(Pˣ, P̃, PartObs()) +set_observations!(model_setup, [L for _ in P̃], [Σ for _ in P̃], obs_vals, obs_time) # uses default fpt +set_imputation_grid!(model_setup, 1/1000) +set_x0_prior!(model_setup, + GsnStartingPt(x0, @SMatrix [0.01 0.0; + 0.0 0.01;]), + x0) +set_auxiliary!(model_setup; skip_for_save=10^0, + adaptive_prop=NoAdaptation()) +initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt(100)) +#:step, :scale, :min, :max, :trgt, :offset +readj = (100, 0.001, 0.001, 999.9, 0.4, 50) + +mcmc_setup = MCMCSetup( + Imputation(NoBlocking(), 0.99, Vern7()), + ParamUpdate(ConjugateUpdt(), [1,2], θ_init, nothing, + MvNormal(fill(0.0, 2), diagm(0=>fill(1000.0, 2))), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1,2]))), + ParamUpdate(MetropolisHastingsUpdt(), 1, θ_init, + UniformRandomWalk(0.1, true), ImproperPrior(), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 1)), readj), + ParamUpdate(MetropolisHastingsUpdt(), 2, θ_init, + UniformRandomWalk(0.01, true), ImproperPrior(), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 2)), readj), + ParamUpdate(MetropolisHastingsUpdt(), 3, θ_init, + UniformRandomWalk(0.1, true), ImproperPrior(), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 3)), readj), + ParamUpdate(MetropolisHastingsUpdt(), 4, θ_init, + UniformRandomWalk(0.1, true), ImproperPrior(), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 4)), readj), + ) + +schedule = MCMCSchedule(4*10^4, [[1],[2],[5]], #[[1],[2], [5]], + (save=1*10^3, verbose=10^2, warm_up=100, + readjust=(x->x%100==0), fuse=(x->false))) + +Random.seed!(4) +out = mcmc(mcmc_setup, schedule, model_setup) +ws = out[2] + +θs = ws.θ_chain +±(a, b) = a - b, a + b +beta, gamma = [median(getindex.(θs, i)) for i in [1,2]] .± [std(getindex.(θs, i)) for i in [1,2]] + +@show beta +@show gamma + +mean(getindex.(θs, 1)./getindex.(θs, 2)) +std(getindex.(θs, 1)./getindex.(θs, 2)) + +lines(getindex.(θs, 1)) +lines(getindex.(θs, 2)) +lines(getindex.(θs, 3)) + +include(joinpath(SRC_DIR, DIR, "plotting_fns.jl")) +plot_chains(ws; truth=θˣ) +#= +plot_paths(out[1], out[2], schedule; obs=(times=obs_time[2:end], + vals=[[v[1] for v in obs_vals[2:end]], + [v[2] for v in obs_vals[2:end]]], indices=[2,3])) +=# diff --git a/src/BridgeSDEInference.jl b/src/BridgeSDEInference.jl index f0bf8d0..8b9a3c8 100644 --- a/src/BridgeSDEInference.jl +++ b/src/BridgeSDEInference.jl @@ -103,6 +103,7 @@ include(joinpath(_DIR, "lorenz_system_const_vola.jl")) include(joinpath(_DIR, "prokaryotic_autoregulatory_gene_network.jl")) include(joinpath(_DIR, "Jansen_and_Rit_simple.jl")) include(joinpath(_DIR, "lotka_volterra.jl")) +include(joinpath(_DIR, "sir.jl")) _DIR = "mcmc" include(joinpath(_DIR, "priors.jl")) diff --git a/src/BridgeSDEInference_for_tests.jl b/src/BridgeSDEInference_for_tests.jl index 91e9f74..2e31662 100644 --- a/src/BridgeSDEInference_for_tests.jl +++ b/src/BridgeSDEInference_for_tests.jl @@ -36,6 +36,8 @@ include(joinpath(_DIR, "lorenz_system.jl")) include(joinpath(_DIR, "lorenz_system_const_vola.jl")) include(joinpath(_DIR, "prokaryotic_autoregulatory_gene_network.jl")) include(joinpath(_DIR, "lotka_volterra.jl")) +include(joinpath(_DIR, "sir.jl")) + _DIR = "mcmc" include(joinpath(_DIR, "priors.jl")) diff --git a/src/examples/sir.jl b/src/examples/sir.jl new file mode 100644 index 0000000..5bd863e --- /dev/null +++ b/src/examples/sir.jl @@ -0,0 +1,87 @@ +using Bridge +using StaticArrays +import Bridge: b, σ, B, β, a, constdiff +const ℝ = SVector{N,T} where {N,T} +import Base.display + +struct SIR{T} <: ContinuousTimeProcess{SVector{2,T}} + α::T + β::T + σ1::T + σ2::T + k::T +end + +b(t, u, P::SIR) = @SVector [P.α*(P.k - u[1] - u[2])*u[1] - P.β*u[1], P.β*u[1]] +σ(t, u, P::SIR) = @SMatrix Float64[ + -P.σ1*(P.k - u[1] - u[2])*u[1] -P.σ2*u[1] + 0.0 P.σ2*u[1] + ] +a(t, u, P::SIR) = σ(t, u, P)*σ(t, u, P)' +constdiff(::SIR) = false +clone(P::SIR, θ) = SIR(θ...) +params(P::SIR) = [P.α, P.β, P.σ1, P.σ2, P.k] +#domain(P::SIR) = LowerBoundedDomain((0.0, 0.0), (1,2)) +domain(P::SIR) = LowerBoundedDomain((0.0,), (1,)) + + +# <--------------------------------------------- +# this is optional, needed for conjugate updates +phi(::Val{0}, t, u, P::SIR) = (zero(u[1]), zero(u[2])) +#[P.α*(P.k - u[1] - u[2])*u[1] - P.β*u[1], P.β*u[1]] +phi(::Val{1}, t, u, P::SIR) = ((P.k - u[1] - u[2])*u[1], zero(u[2])) +phi(::Val{2}, t, u, P::SIR) = (-u[1], u[1]) +phi(::Val{3}, t, x, P::SIR) = (zero(x[1]), zero(x[2])) +phi(::Val{4}, t, x, P::SIR) = (zero(x[1]), zero(x[2])) +phi(::Val{5}, t, x, P::SIR) = (zero(x[1]), zero(x[2])) + + + +nonhypo(P::SIR, x) = x +@inline hypo_a_inv(P::SIR, t, x) = inv(a(t, x, P)) +num_non_hypo(P::Type{<:SIR}) = 2 + + +struct SIRAux{T,S1,S2} <: ContinuousTimeProcess{SVector{2,T}} + α::T + β::T + σ1::T + σ2::T + k::T + t::Float64 + u::S1 + T::Float64 + v::S2 +end + + +function B(t, P::SIRAux) +# b(t, u, P::SIR) = @SVector [P.α*(P.k - u[1] - u[2])*u[1] - P.β*u[1], P.β*u[1]] + 1000*@SMatrix [P.α -P.β; + P.β -0.0] +end + +# mean = ℝ{2}(P.γ/P.δ, P.α/P.β) +function β(t, P::SIRAux) + ℝ{2}(0.0, 0.0) +end + +function σ(t, P::SIRAux) + @SMatrix Float64[ + -P.σ1 -P.σ2 + 0.0 P.σ2 + ] +end +σ(t, x, P::SIRAux) = σ(t, P) + +depends_on_params(::SIRAux) = (3, 4, 5) + +constdiff(::SIRAux) = true +b(t, x, P::SIRAux) = B(t,P) * x + β(t,P) +a(t, P::SIRAux) = σ(t,P) * σ(t, P)' + + +clone(P::SIRAux, θ) = SIRAux(θ..., P.t, P.u, P.T, P.v) + +clone(P::SIRAux, θ, v) = SIRAux(θ..., P.t, v, P.T, v) +params(P::SIRAux) = [P.α, P.β, P.σ1, P.σ2, P.k] From 75af874db3ad88a2e9b908aa42ab25ea58fd1518 Mon Sep 17 00:00:00 2001 From: mschauer Date: Thu, 30 Jan 2020 23:37:20 +0100 Subject: [PATCH 02/29] Tuning --- scripts/inference/sir.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/scripts/inference/sir.jl b/scripts/inference/sir.jl index 1212516..342062f 100644 --- a/scripts/inference/sir.jl +++ b/scripts/inference/sir.jl @@ -18,7 +18,8 @@ include(joinpath(SRC_DIR, DIR, "utility_functions.jl")) Random.seed!(4) pop = 50_000_000 K = 1.0 -θˣ = [0.37, 0.05, 0.1, 0.01, K] +θˣ = [0.37, 0.05, 0.05, 0.01, K] + Pˣ = SIR(θˣ...) x0, dt, T = ℝ{2}(1/pop, 0.), 1/5000, 30.0 @@ -37,8 +38,8 @@ Pˣ = SIR(θ_init...) length(XX.tt) skip = 20000 -Σdiagel = 10/pop -Σ = SMatrix{2,2}(1.0I)*Σdiagel +#Σdiagel = +Σ = @SMatrix[5.0 0.0; 0.0 1.0]/pop L = @SMatrix[1.0 0.0; 0.0 1.0] obs_time, obs_vals = XX.tt[1:skip:end], [rand(Gaussian(L*x, Σ)) for x in XX.yy[1:skip:end]] @@ -85,7 +86,7 @@ mcmc_setup = MCMCSetup( UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 4)), readj), ) -schedule = MCMCSchedule(4*10^4, [[1],[2],[5]], #[[1],[2], [5]], +schedule = MCMCSchedule(2*10^4, [[1],[2]], #[[1],[2], [5]], (save=1*10^3, verbose=10^2, warm_up=100, readjust=(x->x%100==0), fuse=(x->false))) @@ -95,13 +96,13 @@ ws = out[2] θs = ws.θ_chain ±(a, b) = a - b, a + b -beta, gamma = [median(getindex.(θs, i)) for i in [1,2]] .± [std(getindex.(θs, i)) for i in [1,2]] +beta, gamma, s1 = [median(getindex.(θs, i)) for i in [1,2,3]] .± [std(getindex.(θs, i)) for i in [1,2, 3]] +R0 = mean(getindex.(θs, 1)./getindex.(θs, 2)) ± std(getindex.(θs, 1)./getindex.(θs, 2)) @show beta @show gamma - -mean(getindex.(θs, 1)./getindex.(θs, 2)) -std(getindex.(θs, 1)./getindex.(θs, 2)) +@show s1 +@show R0 lines(getindex.(θs, 1)) lines(getindex.(θs, 2)) From 676a21ab84248221b6d2b267dd6e8980ae95192b Mon Sep 17 00:00:00 2001 From: mschauer Date: Fri, 31 Jan 2020 04:14:32 +0100 Subject: [PATCH 03/29] Fixes SIR --- scripts/inference/sir.jl | 29 ++++++++++++++++------------- src/examples/sir.jl | 12 ++++++------ src/general/readjustments.jl | 6 ++++++ 3 files changed, 28 insertions(+), 19 deletions(-) diff --git a/scripts/inference/sir.jl b/scripts/inference/sir.jl index 342062f..58d4d52 100644 --- a/scripts/inference/sir.jl +++ b/scripts/inference/sir.jl @@ -15,14 +15,15 @@ using Bridge DIR = "auxiliary" include(joinpath(SRC_DIR, DIR, "data_simulation_fns.jl")) include(joinpath(SRC_DIR, DIR, "utility_functions.jl")) -Random.seed!(4) +Random.seed!(2) pop = 50_000_000 K = 1.0 -θˣ = [0.37, 0.05, 0.05, 0.01, K] +#θˣ = [0.37, 0.05, 0.05, 0.01, K] +θˣ = [0.37, 0.05, 0.3/sqrt(pop), 1/sqrt(pop), K] Pˣ = SIR(θˣ...) -x0, dt, T = ℝ{2}(1/pop, 0.), 1/5000, 30.0 +x0, dt, T = ℝ{2}(1/pop, 0.), 1/10000, 30.0 tt = 0.0:dt:T XX, _ = simulate_segment(ℝ{2}(1.0, 0.0), x0, Pˣ, tt) @@ -39,18 +40,19 @@ length(XX.tt) skip = 20000 #Σdiagel = -Σ = @SMatrix[5.0 0.0; 0.0 1.0]/pop +Σ = @SMatrix[1.0 0.0; 0.0 0.5]/pop L = @SMatrix[1.0 0.0; 0.0 1.0] obs_time, obs_vals = XX.tt[1:skip:end], [rand(Gaussian(L*x, Σ)) for x in XX.yy[1:skip:end]] days = [0.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0] -cases = [1.0 0.0; 62.0 0.0; 121.0 0.1; 198.0 3.0; 291.0 6.0; 440.0 9.0; 571.0 17.0; 830.0 25.0; 1287.0 41.0; 1975.0 56.0; 2744.0 80.0; 4515.0 106.0; 5974.0 132.0; 7771.0 170.0] +cases = [1.0 0.0; 62.0 0.0; 121.0 0.0; 198.0 3.0; 291.0 6.0; 440.0 9.0; 571.0 17.0; 830.0 25.0; 1287.0 41.0; 1975.0 56.0; 2744.0 80.0; 4515.0 106.0; 5974.0 132.0; 7771.0 170.0] + if true obs_time = days - obs_vals = (1/pop)*reinterpret(SVector{2,Float64}, (cases)') + obs_vals = 1/pop*reinterpret(SVector{2,Float64}, [1 -1; 0 1]*cases') # first coordinate is I + R end - +obs_vals[end] P̃ = [SIRAux(θ_init..., t₀, u, T, v) for (t₀, T, u, v) in zip(obs_time[1:end-1], obs_time[2:end], obs_vals[1:end-1], obs_vals[2:end])] @@ -66,9 +68,10 @@ set_auxiliary!(model_setup; skip_for_save=10^0, initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt(100)) #:step, :scale, :min, :max, :trgt, :offset readj = (100, 0.001, 0.001, 999.9, 0.4, 50) +readj2 = (100, 0.01, 0.05, 0.2, 0.7, 50) mcmc_setup = MCMCSetup( - Imputation(NoBlocking(), 0.99, Vern7()), + Imputation(NoBlocking(), 0.999, Vern7()), ParamUpdate(ConjugateUpdt(), [1,2], θ_init, nothing, MvNormal(fill(0.0, 2), diagm(0=>fill(1000.0, 2))), UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1,2]))), @@ -79,14 +82,14 @@ mcmc_setup = MCMCSetup( UniformRandomWalk(0.01, true), ImproperPrior(), UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 2)), readj), ParamUpdate(MetropolisHastingsUpdt(), 3, θ_init, - UniformRandomWalk(0.1, true), ImproperPrior(), - UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 3)), readj), + UniformRandomWalk(0.1, true), ImproperPosPrior(), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 3)), readj2), ParamUpdate(MetropolisHastingsUpdt(), 4, θ_init, - UniformRandomWalk(0.1, true), ImproperPrior(), - UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 4)), readj), + UniformRandomWalk(0.1, true), ImproperPosPrior(), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 4)), readj2), ) -schedule = MCMCSchedule(2*10^4, [[1],[2]], #[[1],[2], [5]], +schedule = MCMCSchedule(1*10^4, [[1],[2],[5]], #[[1],[2], [5]], (save=1*10^3, verbose=10^2, warm_up=100, readjust=(x->x%100==0), fuse=(x->false))) diff --git a/src/examples/sir.jl b/src/examples/sir.jl index 5bd863e..7223869 100644 --- a/src/examples/sir.jl +++ b/src/examples/sir.jl @@ -3,7 +3,7 @@ using StaticArrays import Bridge: b, σ, B, β, a, constdiff const ℝ = SVector{N,T} where {N,T} import Base.display - +sq(x) = sqrt(max(x, 2e-10)) struct SIR{T} <: ContinuousTimeProcess{SVector{2,T}} α::T β::T @@ -14,8 +14,8 @@ end b(t, u, P::SIR) = @SVector [P.α*(P.k - u[1] - u[2])*u[1] - P.β*u[1], P.β*u[1]] σ(t, u, P::SIR) = @SMatrix Float64[ - -P.σ1*(P.k - u[1] - u[2])*u[1] -P.σ2*u[1] - 0.0 P.σ2*u[1] + -P.σ1*sq((P.k - u[1] - u[2])*u[1]) -P.σ2*sq(u[1]) + 0.0 P.σ2*sq.(u[1]) ] a(t, u, P::SIR) = σ(t, u, P)*σ(t, u, P)' constdiff(::SIR) = false @@ -57,8 +57,8 @@ end function B(t, P::SIRAux) # b(t, u, P::SIR) = @SVector [P.α*(P.k - u[1] - u[2])*u[1] - P.β*u[1], P.β*u[1]] - 1000*@SMatrix [P.α -P.β; - P.β -0.0] + @SMatrix [(P.α) -(P.β); + (P.β) -0.0] end # mean = ℝ{2}(P.γ/P.δ, P.α/P.β) @@ -67,7 +67,7 @@ function β(t, P::SIRAux) end function σ(t, P::SIRAux) - @SMatrix Float64[ + sq(0.0001)*@SMatrix Float64[ -P.σ1 -P.σ2 0.0 P.σ2 ] diff --git a/src/general/readjustments.jl b/src/general/readjustments.jl index e53bf42..5848535 100644 --- a/src/general/readjustments.jl +++ b/src/general/readjustments.jl @@ -6,8 +6,14 @@ function named_readjust(p) (step=p[1], scale=p[2], min=p[3], max=p[4], trgt=p[5], offset=p[6]) end +""" +δ decreases roughly proportional to scale/sqrt(iteration) +""" compute_δ(p, mcmc_iter) = p.scale/sqrt(max(1.0, mcmc_iter/p.step-p.offset)) +""" +ϵ is moved by δ to adapt to target acceptance rate +""" function compute_ϵ(ϵ_old, p, a_r, δ, flip=1.0, f=identity, finv=identity) ϵ = finv(f(ϵ_old) + flip*(2*(a_r > p.trgt)-1)*δ) ϵ = max(min(ϵ, p.max), p.min) # trim excessive updates From c1d01a5d27c9f6ad6778c781219f9bdcb69fbfe7 Mon Sep 17 00:00:00 2001 From: mschauer Date: Sat, 1 Feb 2020 10:27:09 +0100 Subject: [PATCH 04/29] Experiment --- scripts/inference/sir.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/inference/sir.jl b/scripts/inference/sir.jl index 58d4d52..f0f2e8e 100644 --- a/scripts/inference/sir.jl +++ b/scripts/inference/sir.jl @@ -89,7 +89,7 @@ mcmc_setup = MCMCSetup( UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 4)), readj2), ) -schedule = MCMCSchedule(1*10^4, [[1],[2],[5]], #[[1],[2], [5]], +schedule = MCMCSchedule(10*10^4, [[1],[2],[5]], #[[1],[2], [5]], (save=1*10^3, verbose=10^2, warm_up=100, readjust=(x->x%100==0), fuse=(x->false))) From 93fa09bdf12fa271cd1a0d0dbdb4f12bdd500101 Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 11 Mar 2020 11:02:42 +0100 Subject: [PATCH 05/29] js particle animation --- scripts/nextjournal/domserv.jl | 170 +++++++++++++++++++++++++++++++ scripts/nextjournal/particle.css | 5 + 2 files changed, 175 insertions(+) create mode 100644 scripts/nextjournal/domserv.jl create mode 100644 scripts/nextjournal/particle.css diff --git a/scripts/nextjournal/domserv.jl b/scripts/nextjournal/domserv.jl new file mode 100644 index 0000000..23dc9f8 --- /dev/null +++ b/scripts/nextjournal/domserv.jl @@ -0,0 +1,170 @@ +using JSServe, WGLMakie, AbstractPlotting +using JSServe: JSServe.DOM, @js_str, onjs +global three, scene + +using StaticArrays +using Colors +using Random +rebirth(α, R) = x -> (rand() > α ? x : (2rand(typeof(x)) .- 1).*R) +const 𝕏 = SVector + + + +using Hyperscript, Markdown +using JSServe, Observables +using JSServe: Application, Session, evaljs, linkjs, div, active_sessions, Asset +using JSServe: @js_str, onjs, Button, TextField, Slider, JSString, Dependency, with_session +using JSServe.DOM + + +#= +function dom_handler3(session, request) + R = 𝕏(1.5,6.0) + k = 50 + n = 500 + + + #c = 1 + ms = 0.08 + #i = 1; + + + x0 = zero(𝕏{2,Float32}) + limits = FRect(-R[1], -R[2], 2R[1], 2R[2]) + p = Scene(resolution=(800,800), limits=limits, backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22)) + + pts = [rebirth(1.0, R)(x0) for i in 1:n] + xraw = [copy(pts) for i in 1:k] + x = [Node(xraw[i]) for i in 1:k] + col = [Node(RGBA{Float32}(0.0, 0.0, 0.0, 0.0)) for i in 1:k] + + for i in randperm(k) + scatter!(p, x[i], color = col[i], markersize = ms, #=show_axis = true,limits=limits,=# glowwidth = 0.005, glowcolor = :white) + end + #update_cam!(p, limits) + axis = p[Axis] + axis[:grid, :linewidth] = (1, 1) + axis[:grid, :linecolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.5),RGBA{Float32}(0.5, 0.7, 1.0, 0.5)) + axis[:names][:textsize] = (0.0,0.0) + axis[:ticks, :textcolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.8),RGBA{Float32}(0.5, 0.7, 1.0, 0.8)) + r = range(-R[1], 2R[1], length=100) + + scene = p + three, canvas = WGLMakie.three_display(session, scene) + return DOM.div(canvas) +end + +JSServe.with_session() do session, request + dom_handler2(session, request) +end +=# + +function dom_handler(session, request) + global three, scene + + # slider and field for sigma + sliders = JSServe.Slider(0.01:0.01:1) + nrs = JSServe.NumberInput(0.0) + linkjs(session, sliders.value, nrs.value) + + # time wheel ;-) + button = JSServe.Slider(1:109) + + # init + R = 𝕏(1.5,6.0) + limits = FRect(-R[1], -R[2], 2R[1], 2R[2]) + n = 800 + K = 80 + dt = 0.001 + sqrtdt = sqrt(dt) + + particlecss = Asset(joinpath(@__DIR__,"particle.css")) + + # init javascript + evaljs(session, js""" + console.log("Hello"); + iter = 1; + si = 0.0; + //setInterval( + """) + + #css("#slider", var"background-image"="linear-gradient(blue, green, blue)") + + scene = scatter(repeat(2randn(n), outer=K), repeat(2randn(n),outer=K), color = fill(:white, n*K), + backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = 0.03, + resolution=(600,600), limits = limits, + ) + axis = scene[Axis] + axis[:grid, :linewidth] = (0.3, 0.3) + axis[:grid, :linecolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.3),RGBA{Float32}(0.5, 0.7, 1.0, 0.3)) + axis[:names][:textsize] = (0.0,0.0) + axis[:ticks, :textcolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.5),RGBA{Float32}(0.5, 0.7, 1.0, 0.5)) + + + splot = scene[end] + three, canvas = WGLMakie.three_display(session, scene) + js_scene = WGLMakie.to_jsscene(three, scene) + mesh = js_scene.getObjectByName(string(objectid(splot))) + + onjs(session, sliders.value, js"""function (value){ + si = value; + }""") + + onjs(session, button.value, js"""function (value){ + function randn_bm() { + var u = 0, v = 0; + while(u === 0) u = Math.random(); //Converting [0,1) to (0,1) + while(v === 0) v = Math.random(); + return Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v ); + } + var mu = 0.2; + var mesh = $(mesh); + var K = $(K); + var n = $(n); + var dt = $(dt); + console.log(iter++); + var sqrtdt = $(sqrtdt); + + k = iter%K; + var positions = mesh.geometry.attributes.offset.array; + var color = mesh.geometry.attributes.color.array; + console.log(color.length); + for ( var i = 0; i < n; i++ ) { + inew = k*2*n + 2*i; + iold = ((K + k - 1)%K)*2*n + 2*i; + positions[inew] = (1 - mu*dt)*positions[iold] - 3*dt*positions[iold+1] + si*sqrtdt*randn_bm(); // x + positions[inew+1] = (1 - mu*dt)*positions[iold+1] + 3*dt*positions[iold] + si*sqrtdt*randn_bm(); + color[k*4*n + 4*i] = 1.0; + color[k*4*n + 4*i + 1] = 1.0; + color[k*4*n + 4*i + 2] = 1.0; + color[k*4*n + 4*i + 3] = 1.0; + + } + for ( var k = 0; k < K; k++ ) { + for ( var i = 0; i < n; i++ ) { + color[k*4*n + 4*i + 3] = 0.98*color[k*4*n + 4*i + 3]; + } + } + mesh.geometry.attributes.color.needsUpdate = true; + mesh.geometry.attributes.offset.needsUpdate = true; + + }""") + dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters", DOM.div(sliders, id="slider")), + DOM.p(nrs), DOM.p(button)) +# JSServe.onload(session, dom, js""" +# iter = 1; +# """) + println("running...") + dom +end + + +app = JSServe.Application( + dom_handler, + get(ENV, "WEBIO_SERVER_HOST_URL", "127.0.0.1"), + parse(Int, get(ENV, "WEBIO_HTTP_PORT", "8081")), + verbose = false +) +cl() = (close(app), "stopped") +println("Done.") +cl() diff --git a/scripts/nextjournal/particle.css b/scripts/nextjournal/particle.css new file mode 100644 index 0000000..e985711 --- /dev/null +++ b/scripts/nextjournal/particle.css @@ -0,0 +1,5 @@ +#slider { + background-color: #778899; + width: 4cm; + background-image: linear-gradient(to right, #8888FF, #FFFF88, #8888FF); +} From aef9ae96e52563af8e4565b351d7317472667a8c Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 11 Mar 2020 11:37:25 +0100 Subject: [PATCH 06/29] Restart particles --- scripts/nextjournal/domserv.jl | 51 ++++++---------------------------- 1 file changed, 9 insertions(+), 42 deletions(-) diff --git a/scripts/nextjournal/domserv.jl b/scripts/nextjournal/domserv.jl index 23dc9f8..a31959a 100644 --- a/scripts/nextjournal/domserv.jl +++ b/scripts/nextjournal/domserv.jl @@ -17,48 +17,6 @@ using JSServe: @js_str, onjs, Button, TextField, Slider, JSString, Dependency, w using JSServe.DOM -#= -function dom_handler3(session, request) - R = 𝕏(1.5,6.0) - k = 50 - n = 500 - - - #c = 1 - ms = 0.08 - #i = 1; - - - x0 = zero(𝕏{2,Float32}) - limits = FRect(-R[1], -R[2], 2R[1], 2R[2]) - p = Scene(resolution=(800,800), limits=limits, backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22)) - - pts = [rebirth(1.0, R)(x0) for i in 1:n] - xraw = [copy(pts) for i in 1:k] - x = [Node(xraw[i]) for i in 1:k] - col = [Node(RGBA{Float32}(0.0, 0.0, 0.0, 0.0)) for i in 1:k] - - for i in randperm(k) - scatter!(p, x[i], color = col[i], markersize = ms, #=show_axis = true,limits=limits,=# glowwidth = 0.005, glowcolor = :white) - end - #update_cam!(p, limits) - axis = p[Axis] - axis[:grid, :linewidth] = (1, 1) - axis[:grid, :linecolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.5),RGBA{Float32}(0.5, 0.7, 1.0, 0.5)) - axis[:names][:textsize] = (0.0,0.0) - axis[:ticks, :textcolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.8),RGBA{Float32}(0.5, 0.7, 1.0, 0.8)) - r = range(-R[1], 2R[1], length=100) - - scene = p - three, canvas = WGLMakie.three_display(session, scene) - return DOM.div(canvas) -end - -JSServe.with_session() do session, request - dom_handler2(session, request) -end -=# - function dom_handler(session, request) global three, scene @@ -72,6 +30,7 @@ function dom_handler(session, request) # init R = 𝕏(1.5,6.0) + R1, R2 = R limits = FRect(-R[1], -R[2], 2R[1], 2R[2]) n = 800 K = 80 @@ -85,6 +44,8 @@ function dom_handler(session, request) console.log("Hello"); iter = 1; si = 0.0; + R1 = $(R1); + R2 = $(R2); //setInterval( """) @@ -92,6 +53,7 @@ function dom_handler(session, request) scene = scatter(repeat(2randn(n), outer=K), repeat(2randn(n),outer=K), color = fill(:white, n*K), backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = 0.03, + glowwidth = 0.005, glowcolor = :white, resolution=(600,600), limits = limits, ) axis = scene[Axis] @@ -138,6 +100,11 @@ function dom_handler(session, request) color[k*4*n + 4*i + 1] = 1.0; color[k*4*n + 4*i + 2] = 1.0; color[k*4*n + 4*i + 3] = 1.0; + if (Math.random() < 0.01) + { + positions[inew] = (2*Math.random()-1)*R1; + positions[inew+1] = (2*Math.random()-1)*R2; + } } for ( var k = 0; k < K; k++ ) { From 4f1fcfa7e2af3307a26a9a99d6f1830ea0f47536 Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Wed, 11 Mar 2020 15:23:53 +0100 Subject: [PATCH 07/29] interface with FitzHugh-Nagumo dynamics --- scripts/nextjournal/domserv_fhn.jl | 141 +++++++++++++++++++++++++++++ 1 file changed, 141 insertions(+) create mode 100644 scripts/nextjournal/domserv_fhn.jl diff --git a/scripts/nextjournal/domserv_fhn.jl b/scripts/nextjournal/domserv_fhn.jl new file mode 100644 index 0000000..20d86fa --- /dev/null +++ b/scripts/nextjournal/domserv_fhn.jl @@ -0,0 +1,141 @@ +using JSServe, WGLMakie, AbstractPlotting +using JSServe: JSServe.DOM, @js_str, onjs +global three, scene + +using StaticArrays +using Colors +using Random +rebirth(α, R) = x -> (rand() > α ? x : (2rand(typeof(x)) .- 1).*R) +const 𝕏 = SVector + + + +using Hyperscript, Markdown +using JSServe, Observables +using JSServe: Application, Session, evaljs, linkjs, div, active_sessions, Asset +using JSServe: @js_str, onjs, Button, TextField, Slider, JSString, Dependency, with_session +using JSServe.DOM + + +function dom_handler(session, request) + global three, scene + + # slider and field for sigma + sliders = JSServe.Slider(0.01:0.01:1) + nrs = JSServe.NumberInput(0.0) + linkjs(session, sliders.value, nrs.value) + + # time wheel ;-) + button = JSServe.Slider(1:109) + + # init + R = 𝕏(1.5,6.0) + R1, R2 = R + limits = FRect(-R[1], -R[2], 2R[1], 2R[2]) + n = 800 + K = 80 + dt = 0.001 + sqrtdt = sqrt(dt) + + particlecss = Asset(joinpath(@__DIR__,"particle.css")) + + # init javascript + evaljs(session, js""" + console.log("Hello"); + iter = 1; + eps = 0.1; + s = -0.8; + gamma = 1.5; + beta = 0.0; + si = 0.0; + R1 = $(R1); + R2 = $(R2); + //setInterval( + """) + + #css("#slider", var"background-image"="linear-gradient(blue, green, blue)") + + scene = scatter(repeat(2randn(n), outer=K), repeat(2randn(n),outer=K), color = fill(:white, n*K), + backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = 0.03, + glowwidth = 0.005, glowcolor = :white, + resolution=(600,600), limits = limits, + ) + axis = scene[Axis] + axis[:grid, :linewidth] = (0.3, 0.3) + axis[:grid, :linecolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.3),RGBA{Float32}(0.5, 0.7, 1.0, 0.3)) + axis[:names][:textsize] = (0.0,0.0) + axis[:ticks, :textcolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.5),RGBA{Float32}(0.5, 0.7, 1.0, 0.5)) + + + splot = scene[end] + three, canvas = WGLMakie.three_display(session, scene) + js_scene = WGLMakie.to_jsscene(three, scene) + mesh = js_scene.getObjectByName(string(objectid(splot))) + + onjs(session, sliders.value, js"""function (value){ + si = value; + }""") + + onjs(session, button.value, js"""function (value){ + function randn_bm() { + var u = 0, v = 0; + while(u === 0) u = Math.random(); //Converting [0,1) to (0,1) + while(v === 0) v = Math.random(); + return Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v ); + } + var mesh = $(mesh); + var K = $(K); + var n = $(n); + var dt = $(dt); + console.log(iter++); + var sqrtdt = $(sqrtdt); + + k = iter%K; + var positions = mesh.geometry.attributes.offset.array; + var color = mesh.geometry.attributes.color.array; + console.log(color.length); + for ( var i = 0; i < n; i++ ) { + inew = k*2*n + 2*i; + iold = ((K + k - 1)%K)*2*n + 2*i; + positions[inew] = positions[iold] + dt/eps*((1 - positions[iold]*positions[iold])*positions[iold] - positions[iold+1] - s); // x + positions[inew+1] = positions[iold+1] + dt*(-positions[iold+1] + gamma*positions[iold] + beta) + si*sqrtdt*randn_bm(); + color[k*4*n + 4*i] = 1.0; + color[k*4*n + 4*i + 1] = 1.0; + color[k*4*n + 4*i + 2] = 1.0; + color[k*4*n + 4*i + 3] = 1.0; + if (Math.random() < 0.01) + { + positions[inew] = (2*Math.random()-1)*R1; + positions[inew+1] = (2*Math.random()-1)*R2; + } + + } + for ( var k = 0; k < K; k++ ) { + for ( var i = 0; i < n; i++ ) { + color[k*4*n + 4*i + 3] = 0.98*color[k*4*n + 4*i + 3]; + } + } + mesh.geometry.attributes.color.needsUpdate = true; + mesh.geometry.attributes.offset.needsUpdate = true; + + }""") + dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters", DOM.div(sliders, id="slider")), + DOM.p(nrs), DOM.p(button)) +# JSServe.onload(session, dom, js""" +# iter = 1; +# """) + println("running...") + dom +end + + +app = JSServe.Application( + dom_handler, + get(ENV, "WEBIO_SERVER_HOST_URL", "127.0.0.1"), + parse(Int, get(ENV, "WEBIO_HTTP_PORT", "8081")), + verbose = false +) + +#cl() = (close(app), "stopped") +#println("Done.") +#cl() From 5c91968d63fa3b1848574ccf4b02a426778fa5c2 Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Wed, 11 Mar 2020 15:27:22 +0100 Subject: [PATCH 08/29] Fix sir diffusion and auxiliary diffusion --- src/examples/sir.jl | 43 ++++++++++++++++++++++++++++++------------- 1 file changed, 30 insertions(+), 13 deletions(-) diff --git a/src/examples/sir.jl b/src/examples/sir.jl index 7223869..25ee3c2 100644 --- a/src/examples/sir.jl +++ b/src/examples/sir.jl @@ -9,18 +9,17 @@ struct SIR{T} <: ContinuousTimeProcess{SVector{2,T}} β::T σ1::T σ2::T - k::T end -b(t, u, P::SIR) = @SVector [P.α*(P.k - u[1] - u[2])*u[1] - P.β*u[1], P.β*u[1]] +b(t, u, P::SIR) = @SVector [P.α*(1 - u[1] - u[2])*u[1] - P.β*u[1], P.β*u[1]] σ(t, u, P::SIR) = @SMatrix Float64[ - -P.σ1*sq((P.k - u[1] - u[2])*u[1]) -P.σ2*sq(u[1]) + -P.σ1*sq((1 - u[1] - u[2])*u[1]) -P.σ2*sq(u[1]) 0.0 P.σ2*sq.(u[1]) ] a(t, u, P::SIR) = σ(t, u, P)*σ(t, u, P)' constdiff(::SIR) = false clone(P::SIR, θ) = SIR(θ...) -params(P::SIR) = [P.α, P.β, P.σ1, P.σ2, P.k] +params(P::SIR) = [P.α, P.β, P.σ1, P.σ2] #domain(P::SIR) = LowerBoundedDomain((0.0, 0.0), (1,2)) domain(P::SIR) = LowerBoundedDomain((0.0,), (1,)) @@ -47,7 +46,6 @@ struct SIRAux{T,S1,S2} <: ContinuousTimeProcess{SVector{2,T}} β::T σ1::T σ2::T - k::T t::Float64 u::S1 T::Float64 @@ -55,23 +53,42 @@ struct SIRAux{T,S1,S2} <: ContinuousTimeProcess{SVector{2,T}} end +# function B(t, P::SIRAux) +# # b(t, u, P::SIR) = @SVector [P.α*(P.k - u[1] - u[2])*u[1] - P.β*u[1], P.β*u[1]] +# @SMatrix [(P.α) -(P.β); +# (P.β) -0.0] +# end function B(t, P::SIRAux) -# b(t, u, P::SIR) = @SVector [P.α*(P.k - u[1] - u[2])*u[1] - P.β*u[1], P.β*u[1]] - @SMatrix [(P.α) -(P.β); - (P.β) -0.0] +# b(t, u, P::SIR) = @SVector [P.α*(1 - u[1] - u[2])*u[1] - P.β*u[1], P.β*u[1]] + @SMatrix [(P.α*(1 - P.v[1] - P.v[2]) - (P.β)) 0.0; + (P.β) 0.0] end + + + + # mean = ℝ{2}(P.γ/P.δ, P.α/P.β) function β(t, P::SIRAux) ℝ{2}(0.0, 0.0) end +# function σ(t, P::SIRAux) +# sq(0.0001)*@SMatrix Float64[ +# -P.σ1 -P.σ2 +# 0.0 P.σ2 +# ] +# end + function σ(t, P::SIRAux) - sq(0.0001)*@SMatrix Float64[ - -P.σ1 -P.σ2 - 0.0 P.σ2 - ] + @SMatrix Float64[ + (-P.σ1*(1 - P.v[2] - P.v[1])*P.v[1]) -P.σ2*P.v[1]; + 0.0 P.σ2*P.v[1]] end + + + + σ(t, x, P::SIRAux) = σ(t, P) depends_on_params(::SIRAux) = (3, 4, 5) @@ -84,4 +101,4 @@ a(t, P::SIRAux) = σ(t,P) * σ(t, P)' clone(P::SIRAux, θ) = SIRAux(θ..., P.t, P.u, P.T, P.v) clone(P::SIRAux, θ, v) = SIRAux(θ..., P.t, v, P.T, v) -params(P::SIRAux) = [P.α, P.β, P.σ1, P.σ2, P.k] +params(P::SIRAux) = [P.α, P.β, P.σ1, P.σ2] From 15c36a0a857a9e00f9f16d2891e9e51cf52af7b3 Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Wed, 11 Mar 2020 15:33:50 +0100 Subject: [PATCH 09/29] smoothing --- scripts/inference/sir.jl | 51 ++++++++++++++++++++++------------------ 1 file changed, 28 insertions(+), 23 deletions(-) diff --git a/scripts/inference/sir.jl b/scripts/inference/sir.jl index f0f2e8e..1cddbc8 100644 --- a/scripts/inference/sir.jl +++ b/scripts/inference/sir.jl @@ -17,21 +17,21 @@ include(joinpath(SRC_DIR, DIR, "data_simulation_fns.jl")) include(joinpath(SRC_DIR, DIR, "utility_functions.jl")) Random.seed!(2) pop = 50_000_000 -K = 1.0 -#θˣ = [0.37, 0.05, 0.05, 0.01, K] -θˣ = [0.37, 0.05, 0.3/sqrt(pop), 1/sqrt(pop), K] +#θˣ = [0.37, 0.05, 0.05, 0.01] +θˣ = [0.37, 0.05, 0.3/sqrt(pop), 1.0/sqrt(pop)] Pˣ = SIR(θˣ...) x0, dt, T = ℝ{2}(1/pop, 0.), 1/10000, 30.0 tt = 0.0:dt:T +Random.seed!(1) XX, _ = simulate_segment(ℝ{2}(1.0, 0.0), x0, Pˣ, tt) last(XX)[2]*pop #lines(XX.tt, K .- sum.(XX.yy)) -lines(XX.tt, pop*first.(XX.yy)) -lines!(XX.tt, pop*last.(XX.yy)) +lines(XX.tt, first.(XX.yy), color = :red) +lines!(XX.tt,last.(XX.yy), color = :blue) θ_init = copy(θˣ) Pˣ = SIR(θ_init...) @@ -72,30 +72,35 @@ readj2 = (100, 0.01, 0.05, 0.2, 0.7, 50) mcmc_setup = MCMCSetup( Imputation(NoBlocking(), 0.999, Vern7()), - ParamUpdate(ConjugateUpdt(), [1,2], θ_init, nothing, - MvNormal(fill(0.0, 2), diagm(0=>fill(1000.0, 2))), - UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1,2]))), - ParamUpdate(MetropolisHastingsUpdt(), 1, θ_init, - UniformRandomWalk(0.1, true), ImproperPrior(), - UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 1)), readj), - ParamUpdate(MetropolisHastingsUpdt(), 2, θ_init, - UniformRandomWalk(0.01, true), ImproperPrior(), - UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 2)), readj), - ParamUpdate(MetropolisHastingsUpdt(), 3, θ_init, - UniformRandomWalk(0.1, true), ImproperPosPrior(), - UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 3)), readj2), - ParamUpdate(MetropolisHastingsUpdt(), 4, θ_init, - UniformRandomWalk(0.1, true), ImproperPosPrior(), - UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 4)), readj2), + #ParamUpdate(ConjugateUpdt(), [1,2], θ_init, nothing, + # MvNormal(fill(0.0, 2), diagm(0=>fill(1000.0, 2))), + # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1,2]))), + # ParamUpdate(MetropolisHastingsUpdt(), 1, θ_init, + # UniformRandomWalk(0.1, true), ImproperPrior(), + # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 1)), readj), + # ParamUpdate(MetropolisHastingsUpdt(), 2, θ_init, + # UniformRandomWalk(0.01, true), ImproperPrior(), + # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 2)), readj), + # ParamUpdate(MetropolisHastingsUpdt(), 3, θ_init, + # UniformRandomWalk(0.1, true), ImproperPosPrior(), + # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 3)), readj2), + # ParamUpdate(MetropolisHastingsUpdt(), 4, θ_init, + # UniformRandomWalk(0.1, true), ImproperPosPrior(), + # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 4)), readj2), ) -schedule = MCMCSchedule(10*10^4, [[1],[2],[5]], #[[1],[2], [5]], - (save=1*10^3, verbose=10^2, warm_up=100, +schedule = MCMCSchedule(10^2, [[1]],#[5]], #[[1],[2], [5]], + (save=10^1, verbose=10^2, warm_up=100, readjust=(x->x%100==0), fuse=(x->false))) Random.seed!(4) out = mcmc(mcmc_setup, schedule, model_setup) -ws = out[2] +error("STOP HERE") + + +using Makie +ws = out[1] + θs = ws.θ_chain ±(a, b) = a - b, a + b From d9cd9c33bdfeff280754ddb78680b3a230b2f3e7 Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 11 Mar 2020 15:53:44 +0100 Subject: [PATCH 10/29] Run automatically --- scripts/nextjournal/domserv.jl | 115 +++++++++++++++++---------------- 1 file changed, 59 insertions(+), 56 deletions(-) diff --git a/scripts/nextjournal/domserv.jl b/scripts/nextjournal/domserv.jl index a31959a..143a0a8 100644 --- a/scripts/nextjournal/domserv.jl +++ b/scripts/nextjournal/domserv.jl @@ -39,15 +39,7 @@ function dom_handler(session, request) particlecss = Asset(joinpath(@__DIR__,"particle.css")) - # init javascript - evaljs(session, js""" - console.log("Hello"); - iter = 1; - si = 0.0; - R1 = $(R1); - R2 = $(R2); - //setInterval( - """) + #css("#slider", var"background-image"="linear-gradient(blue, green, blue)") @@ -68,56 +60,67 @@ function dom_handler(session, request) js_scene = WGLMakie.to_jsscene(three, scene) mesh = js_scene.getObjectByName(string(objectid(splot))) - onjs(session, sliders.value, js"""function (value){ - si = value; - }""") - - onjs(session, button.value, js"""function (value){ - function randn_bm() { - var u = 0, v = 0; - while(u === 0) u = Math.random(); //Converting [0,1) to (0,1) - while(v === 0) v = Math.random(); - return Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v ); - } - var mu = 0.2; - var mesh = $(mesh); - var K = $(K); - var n = $(n); - var dt = $(dt); - console.log(iter++); - var sqrtdt = $(sqrtdt); - - k = iter%K; - var positions = mesh.geometry.attributes.offset.array; - var color = mesh.geometry.attributes.color.array; - console.log(color.length); - for ( var i = 0; i < n; i++ ) { - inew = k*2*n + 2*i; - iold = ((K + k - 1)%K)*2*n + 2*i; - positions[inew] = (1 - mu*dt)*positions[iold] - 3*dt*positions[iold+1] + si*sqrtdt*randn_bm(); // x - positions[inew+1] = (1 - mu*dt)*positions[iold+1] + 3*dt*positions[iold] + si*sqrtdt*randn_bm(); - color[k*4*n + 4*i] = 1.0; - color[k*4*n + 4*i + 1] = 1.0; - color[k*4*n + 4*i + 2] = 1.0; - color[k*4*n + 4*i + 3] = 1.0; - if (Math.random() < 0.01) - { - positions[inew] = (2*Math.random()-1)*R1; - positions[inew+1] = (2*Math.random()-1)*R2; - } + # init javascript + evaljs(session, js""" + console.log("Hello"); + iter = 1; + si = 0.0; + R1 = $(R1); + R2 = $(R2); + setInterval( + function (){ + function randn_bm() { + var u = 0, v = 0; + while(u === 0) u = Math.random(); //Converting [0,1) to (0,1) + while(v === 0) v = Math.random(); + return Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v ); + } + var mu = 0.2; + var mesh = $(mesh); + var K = $(K); + var n = $(n); + var dt = $(dt); + console.log(iter++); + var sqrtdt = $(sqrtdt); + + k = iter%K; + var positions = mesh.geometry.attributes.offset.array; + var color = mesh.geometry.attributes.color.array; + console.log(color.length); + for ( var i = 0; i < n; i++ ) { + inew = k*2*n + 2*i; + iold = ((K + k - 1)%K)*2*n + 2*i; + positions[inew] = (1 - mu*dt)*positions[iold] - 3*dt*positions[iold+1] + si*sqrtdt*randn_bm(); // x + positions[inew+1] = (1 - mu*dt)*positions[iold+1] + 3*dt*positions[iold] + si*sqrtdt*randn_bm(); + color[k*4*n + 4*i] = 1.0; + color[k*4*n + 4*i + 1] = 1.0; + color[k*4*n + 4*i + 2] = 1.0; + color[k*4*n + 4*i + 3] = 1.0; + if (Math.random() < 0.01) + { + positions[inew] = (2*Math.random()-1)*R1; + positions[inew+1] = (2*Math.random()-1)*R2; + } + + } + for ( var k = 0; k < K; k++ ) { + for ( var i = 0; i < n; i++ ) { + color[k*4*n + 4*i + 3] = 0.98*color[k*4*n + 4*i + 3]; + } + } + mesh.geometry.attributes.color.needsUpdate = true; + mesh.geometry.attributes.offset.needsUpdate = true; - } - for ( var k = 0; k < K; k++ ) { - for ( var i = 0; i < n; i++ ) { - color[k*4*n + 4*i + 3] = 0.98*color[k*4*n + 4*i + 3]; } - } - mesh.geometry.attributes.color.needsUpdate = true; - mesh.geometry.attributes.offset.needsUpdate = true; + , 50); + """) + onjs(session, sliders.value, js"""function (value){ + si = value; }""") - dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters", DOM.div(sliders, id="slider")), - DOM.p(nrs), DOM.p(button)) + + dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters"), DOM.div(sliders, id="slider"), + DOM.p(nrs)) # JSServe.onload(session, dom, js""" # iter = 1; # """) @@ -134,4 +137,4 @@ app = JSServe.Application( ) cl() = (close(app), "stopped") println("Done.") -cl() +#cl() From 36fac3add578fffa9a59cc1b32ec59fe55192397 Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Wed, 11 Mar 2020 15:55:43 +0100 Subject: [PATCH 11/29] try multiple sliders --- scripts/nextjournal/domserv_fhn.jl | 21 +++++++++++++++------ scripts/nextjournal/particle.css | 8 +++++++- 2 files changed, 22 insertions(+), 7 deletions(-) diff --git a/scripts/nextjournal/domserv_fhn.jl b/scripts/nextjournal/domserv_fhn.jl index 20d86fa..efe88aa 100644 --- a/scripts/nextjournal/domserv_fhn.jl +++ b/scripts/nextjournal/domserv_fhn.jl @@ -21,9 +21,15 @@ function dom_handler(session, request) global three, scene # slider and field for sigma - sliders = JSServe.Slider(0.01:0.01:1) - nrs = JSServe.NumberInput(0.0) - linkjs(session, sliders.value, nrs.value) + slider1 = JSServe.Slider(0.01:0.01:1) + nrs1 = JSServe.NumberInput(0.0) + linkjs(session, slider1.value, nrs1.value) + + # slider and field for beta + slider2 = JSServe.Slider(0.01:0.01:10) + nrs2 = JSServe.NumberInput(0.0) + linkjs(session, slider2.value, nrs2.value) + # time wheel ;-) button = JSServe.Slider(1:109) @@ -72,9 +78,12 @@ function dom_handler(session, request) js_scene = WGLMakie.to_jsscene(three, scene) mesh = js_scene.getObjectByName(string(objectid(splot))) - onjs(session, sliders.value, js"""function (value){ + onjs(session, slider1.value, js"""function (value){ si = value; }""") + onjs(session, slider2.value, js"""function (value){ + beta = value; + }""") onjs(session, button.value, js"""function (value){ function randn_bm() { @@ -119,8 +128,8 @@ function dom_handler(session, request) mesh.geometry.attributes.offset.needsUpdate = true; }""") - dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters", DOM.div(sliders, id="slider")), - DOM.p(nrs), DOM.p(button)) + dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters", DOM.div(slider1, id="slider1"), DOM.div(slider2, id="slider2")), + DOM.p(nrs1), DOM.p(nrs2), DOM.p(button)) # JSServe.onload(session, dom, js""" # iter = 1; # """) diff --git a/scripts/nextjournal/particle.css b/scripts/nextjournal/particle.css index e985711..af01d28 100644 --- a/scripts/nextjournal/particle.css +++ b/scripts/nextjournal/particle.css @@ -1,4 +1,10 @@ -#slider { +#slider1 { + background-color: #778899; + width: 4cm; + background-image: linear-gradient(to right, #8888FF, #FFFF88, #8888FF); +} + +#slider2 { background-color: #778899; width: 4cm; background-image: linear-gradient(to right, #8888FF, #FFFF88, #8888FF); From 5452f6a16feb8f4779827e991fdca80be057cb2d Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Wed, 11 Mar 2020 16:18:31 +0100 Subject: [PATCH 12/29] update --- scripts/nextjournal/domserv_fhn.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/scripts/nextjournal/domserv_fhn.jl b/scripts/nextjournal/domserv_fhn.jl index efe88aa..b97b675 100644 --- a/scripts/nextjournal/domserv_fhn.jl +++ b/scripts/nextjournal/domserv_fhn.jl @@ -9,7 +9,6 @@ rebirth(α, R) = x -> (rand() > α ? x : (2rand(typeof(x)) .- 1).*R) const 𝕏 = SVector - using Hyperscript, Markdown using JSServe, Observables using JSServe: Application, Session, evaljs, linkjs, div, active_sessions, Asset From b901196ba6d0a1405a3bbbfd85cebe9d410d1e5b Mon Sep 17 00:00:00 2001 From: mschauer Date: Thu, 12 Mar 2020 13:29:12 +0100 Subject: [PATCH 13/29] Experiment with kline --- scripts/nextjournal/domserv.jl | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/scripts/nextjournal/domserv.jl b/scripts/nextjournal/domserv.jl index 143a0a8..dc54b0a 100644 --- a/scripts/nextjournal/domserv.jl +++ b/scripts/nextjournal/domserv.jl @@ -38,13 +38,9 @@ function dom_handler(session, request) sqrtdt = sqrt(dt) particlecss = Asset(joinpath(@__DIR__,"particle.css")) - - - - #css("#slider", var"background-image"="linear-gradient(blue, green, blue)") - - scene = scatter(repeat(2randn(n), outer=K), repeat(2randn(n),outer=K), color = fill(:white, n*K), - backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = 0.03, + ms = 0.03 + global scene = scatter(repeat(2randn(n), outer=K), repeat(2randn(n),outer=K), color = fill(:white, n*K), + backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = ms, glowwidth = 0.005, glowcolor = :white, resolution=(600,600), limits = limits, ) @@ -56,9 +52,13 @@ function dom_handler(session, request) splot = scene[end] + scatter!(scene, -R1:0.01:R1, sin.(-R1:0.01:R1), color = RGBA{Float32}(0.5, 0.7, 1.0, 0.8), markersize=ms) + kplot = scene[end] + three, canvas = WGLMakie.three_display(session, scene) js_scene = WGLMakie.to_jsscene(three, scene) mesh = js_scene.getObjectByName(string(objectid(splot))) + mesh2 = js_scene.getObjectByName(string(objectid(kplot))) # init javascript evaljs(session, js""" @@ -117,6 +117,16 @@ function dom_handler(session, request) onjs(session, sliders.value, js"""function (value){ si = value; + var mesh = $(mesh2); + var positions = mesh.geometry.attributes.offset.array; + var color = mesh.geometry.attributes.color.array; + + for ( var i = 0, l = positions.length; i < l; i += 2 ) { + positions[i+1] = si*Math.sin(positions[i]); + } + mesh.geometry.attributes.offset.needsUpdate = true; + //mesh.geometry.attributes.color.needsUpdate = true; + }""") dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters"), DOM.div(sliders, id="slider"), @@ -137,4 +147,5 @@ app = JSServe.Application( ) cl() = (close(app), "stopped") println("Done.") -#cl() +# +cl() From d540cfe0c0757c7a0d5e1875586ea809ba0f9df4 Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Tue, 17 Mar 2020 14:59:31 +0100 Subject: [PATCH 14/29] Fixes. Parameter estimation not working yet --- scripts/inference/sir.jl | 46 +++++++++++++++++++--------------------- 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/scripts/inference/sir.jl b/scripts/inference/sir.jl index 1cddbc8..9ad3f74 100644 --- a/scripts/inference/sir.jl +++ b/scripts/inference/sir.jl @@ -18,11 +18,12 @@ include(joinpath(SRC_DIR, DIR, "utility_functions.jl")) Random.seed!(2) pop = 50_000_000 #θˣ = [0.37, 0.05, 0.05, 0.01] +#α, β, σ1, σ2 θˣ = [0.37, 0.05, 0.3/sqrt(pop), 1.0/sqrt(pop)] Pˣ = SIR(θˣ...) -x0, dt, T = ℝ{2}(1/pop, 0.), 1/10000, 30.0 +x0, dt, T = ℝ{2}(1/pop, 0.), 1/10000, 10.0 tt = 0.0:dt:T Random.seed!(1) @@ -30,14 +31,14 @@ XX, _ = simulate_segment(ℝ{2}(1.0, 0.0), x0, Pˣ, tt) last(XX)[2]*pop #lines(XX.tt, K .- sum.(XX.yy)) -lines(XX.tt, first.(XX.yy), color = :red) -lines!(XX.tt,last.(XX.yy), color = :blue) +lines(XX.tt, first.(XX.yy), color = :red) #infected +lines!(XX.tt,last.(XX.yy), color = :blue) #recovered + dead people -θ_init = copy(θˣ) +θ_init = [0.5, 0.1, 0.3/sqrt(pop), 1.0/sqrt(pop)] Pˣ = SIR(θ_init...) length(XX.tt) -skip = 20000 +skip = 2000 #Σdiagel = Σ = @SMatrix[1.0 0.0; 0.0 0.5]/pop @@ -70,27 +71,17 @@ initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt(100)) readj = (100, 0.001, 0.001, 999.9, 0.4, 50) readj2 = (100, 0.01, 0.05, 0.2, 0.7, 50) + mcmc_setup = MCMCSetup( - Imputation(NoBlocking(), 0.999, Vern7()), - #ParamUpdate(ConjugateUpdt(), [1,2], θ_init, nothing, - # MvNormal(fill(0.0, 2), diagm(0=>fill(1000.0, 2))), - # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1,2]))), - # ParamUpdate(MetropolisHastingsUpdt(), 1, θ_init, - # UniformRandomWalk(0.1, true), ImproperPrior(), - # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 1)), readj), - # ParamUpdate(MetropolisHastingsUpdt(), 2, θ_init, - # UniformRandomWalk(0.01, true), ImproperPrior(), - # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 2)), readj), - # ParamUpdate(MetropolisHastingsUpdt(), 3, θ_init, - # UniformRandomWalk(0.1, true), ImproperPosPrior(), - # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 3)), readj2), - # ParamUpdate(MetropolisHastingsUpdt(), 4, θ_init, - # UniformRandomWalk(0.1, true), ImproperPosPrior(), - # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 4)), readj2), + Imputation(NoBlocking(), 0.99, Vern7()), + ParamUpdate(ConjugateUpdt(), [1,2], θ_init, nothing, + MvNormal(fill(0.05, 2), diagm(0=>fill(1000.0, 2))), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1,2]))), + ) -schedule = MCMCSchedule(10^2, [[1]],#[5]], #[[1],[2], [5]], - (save=10^1, verbose=10^2, warm_up=100, +schedule = MCMCSchedule(10^4, [[1],[2]], #[[1],[2], [5]], + (save=10^2, verbose=10^2, warm_up=100, readjust=(x->x%100==0), fuse=(x->false))) Random.seed!(4) @@ -99,7 +90,7 @@ error("STOP HERE") using Makie -ws = out[1] +ws = out[2] θs = ws.θ_chain @@ -107,6 +98,10 @@ ws = out[1] beta, gamma, s1 = [median(getindex.(θs, i)) for i in [1,2,3]] .± [std(getindex.(θs, i)) for i in [1,2, 3]] R0 = mean(getindex.(θs, 1)./getindex.(θs, 2)) ± std(getindex.(θs, 1)./getindex.(θs, 2)) +θs = ws.θ_chain + + + @show beta @show gamma @show s1 @@ -114,6 +109,9 @@ R0 = mean(getindex.(θs, 1)./getindex.(θs, 2)) ± std(getindex.(θs, 1)./getind lines(getindex.(θs, 1)) lines(getindex.(θs, 2)) + + + lines(getindex.(θs, 3)) include(joinpath(SRC_DIR, DIR, "plotting_fns.jl")) From 97b160f0564e2d3f62e2814bc82607d89e249c10 Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Tue, 17 Mar 2020 14:59:56 +0100 Subject: [PATCH 15/29] Fixes. Parameter estimation not working yet --- scripts/inference/sir.jl | 46 +++++++++++++++++++--------------------- src/examples/sir.jl | 4 ++-- 2 files changed, 24 insertions(+), 26 deletions(-) diff --git a/scripts/inference/sir.jl b/scripts/inference/sir.jl index 1cddbc8..9ad3f74 100644 --- a/scripts/inference/sir.jl +++ b/scripts/inference/sir.jl @@ -18,11 +18,12 @@ include(joinpath(SRC_DIR, DIR, "utility_functions.jl")) Random.seed!(2) pop = 50_000_000 #θˣ = [0.37, 0.05, 0.05, 0.01] +#α, β, σ1, σ2 θˣ = [0.37, 0.05, 0.3/sqrt(pop), 1.0/sqrt(pop)] Pˣ = SIR(θˣ...) -x0, dt, T = ℝ{2}(1/pop, 0.), 1/10000, 30.0 +x0, dt, T = ℝ{2}(1/pop, 0.), 1/10000, 10.0 tt = 0.0:dt:T Random.seed!(1) @@ -30,14 +31,14 @@ XX, _ = simulate_segment(ℝ{2}(1.0, 0.0), x0, Pˣ, tt) last(XX)[2]*pop #lines(XX.tt, K .- sum.(XX.yy)) -lines(XX.tt, first.(XX.yy), color = :red) -lines!(XX.tt,last.(XX.yy), color = :blue) +lines(XX.tt, first.(XX.yy), color = :red) #infected +lines!(XX.tt,last.(XX.yy), color = :blue) #recovered + dead people -θ_init = copy(θˣ) +θ_init = [0.5, 0.1, 0.3/sqrt(pop), 1.0/sqrt(pop)] Pˣ = SIR(θ_init...) length(XX.tt) -skip = 20000 +skip = 2000 #Σdiagel = Σ = @SMatrix[1.0 0.0; 0.0 0.5]/pop @@ -70,27 +71,17 @@ initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt(100)) readj = (100, 0.001, 0.001, 999.9, 0.4, 50) readj2 = (100, 0.01, 0.05, 0.2, 0.7, 50) + mcmc_setup = MCMCSetup( - Imputation(NoBlocking(), 0.999, Vern7()), - #ParamUpdate(ConjugateUpdt(), [1,2], θ_init, nothing, - # MvNormal(fill(0.0, 2), diagm(0=>fill(1000.0, 2))), - # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1,2]))), - # ParamUpdate(MetropolisHastingsUpdt(), 1, θ_init, - # UniformRandomWalk(0.1, true), ImproperPrior(), - # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 1)), readj), - # ParamUpdate(MetropolisHastingsUpdt(), 2, θ_init, - # UniformRandomWalk(0.01, true), ImproperPrior(), - # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 2)), readj), - # ParamUpdate(MetropolisHastingsUpdt(), 3, θ_init, - # UniformRandomWalk(0.1, true), ImproperPosPrior(), - # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 3)), readj2), - # ParamUpdate(MetropolisHastingsUpdt(), 4, θ_init, - # UniformRandomWalk(0.1, true), ImproperPosPrior(), - # UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 4)), readj2), + Imputation(NoBlocking(), 0.99, Vern7()), + ParamUpdate(ConjugateUpdt(), [1,2], θ_init, nothing, + MvNormal(fill(0.05, 2), diagm(0=>fill(1000.0, 2))), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1,2]))), + ) -schedule = MCMCSchedule(10^2, [[1]],#[5]], #[[1],[2], [5]], - (save=10^1, verbose=10^2, warm_up=100, +schedule = MCMCSchedule(10^4, [[1],[2]], #[[1],[2], [5]], + (save=10^2, verbose=10^2, warm_up=100, readjust=(x->x%100==0), fuse=(x->false))) Random.seed!(4) @@ -99,7 +90,7 @@ error("STOP HERE") using Makie -ws = out[1] +ws = out[2] θs = ws.θ_chain @@ -107,6 +98,10 @@ ws = out[1] beta, gamma, s1 = [median(getindex.(θs, i)) for i in [1,2,3]] .± [std(getindex.(θs, i)) for i in [1,2, 3]] R0 = mean(getindex.(θs, 1)./getindex.(θs, 2)) ± std(getindex.(θs, 1)./getindex.(θs, 2)) +θs = ws.θ_chain + + + @show beta @show gamma @show s1 @@ -114,6 +109,9 @@ R0 = mean(getindex.(θs, 1)./getindex.(θs, 2)) ± std(getindex.(θs, 1)./getind lines(getindex.(θs, 1)) lines(getindex.(θs, 2)) + + + lines(getindex.(θs, 3)) include(joinpath(SRC_DIR, DIR, "plotting_fns.jl")) diff --git a/src/examples/sir.jl b/src/examples/sir.jl index 25ee3c2..5cefbff 100644 --- a/src/examples/sir.jl +++ b/src/examples/sir.jl @@ -27,8 +27,8 @@ domain(P::SIR) = LowerBoundedDomain((0.0,), (1,)) # <--------------------------------------------- # this is optional, needed for conjugate updates phi(::Val{0}, t, u, P::SIR) = (zero(u[1]), zero(u[2])) -#[P.α*(P.k - u[1] - u[2])*u[1] - P.β*u[1], P.β*u[1]] -phi(::Val{1}, t, u, P::SIR) = ((P.k - u[1] - u[2])*u[1], zero(u[2])) +#[P.α*(1 - u[1] - u[2])*u[1] - P.β*u[1], P.β*u[1]] +phi(::Val{1}, t, u, P::SIR) = ((1 - u[1] - u[2])*u[1], zero(u[2])) phi(::Val{2}, t, u, P::SIR) = (-u[1], u[1]) phi(::Val{3}, t, x, P::SIR) = (zero(x[1]), zero(x[2])) phi(::Val{4}, t, x, P::SIR) = (zero(x[1]), zero(x[2])) From edb77434991229369241fd7b1d5005f41b41780d Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Fri, 20 Mar 2020 21:58:29 +0100 Subject: [PATCH 16/29] add sliders for each parameter and automatic update of the particles --- scripts/nextjournal/domserv_fhn.jl | 167 ++++++++++++++++++----------- scripts/nextjournal/particle.css | 19 ++++ 2 files changed, 124 insertions(+), 62 deletions(-) diff --git a/scripts/nextjournal/domserv_fhn.jl b/scripts/nextjournal/domserv_fhn.jl index b97b675..bee83c9 100644 --- a/scripts/nextjournal/domserv_fhn.jl +++ b/scripts/nextjournal/domserv_fhn.jl @@ -9,6 +9,7 @@ rebirth(α, R) = x -> (rand() > α ? x : (2rand(typeof(x)) .- 1).*R) const 𝕏 = SVector + using Hyperscript, Markdown using JSServe, Observables using JSServe: Application, Session, evaljs, linkjs, div, active_sessions, Asset @@ -29,6 +30,20 @@ function dom_handler(session, request) nrs2 = JSServe.NumberInput(0.0) linkjs(session, slider2.value, nrs2.value) + # slider and field for gamma + slider3 = JSServe.Slider(0.01:0.01:10) + nrs3 = JSServe.NumberInput(0.0) + linkjs(session, slider3.value, nrs3.value) + + # slider and field for s + slider4 = JSServe.Slider(0.01:0.01:10) + nrs4 = JSServe.NumberInput(0.0) + linkjs(session, slider4.value, nrs4.value) + + # slider and field for esp + slider5 = JSServe.Slider(0.01:0.01:10) + nrs5 = JSServe.NumberInput(0.0) + linkjs(session, slider5.value, nrs5.value) # time wheel ;-) button = JSServe.Slider(1:109) @@ -43,25 +58,9 @@ function dom_handler(session, request) sqrtdt = sqrt(dt) particlecss = Asset(joinpath(@__DIR__,"particle.css")) - - # init javascript - evaljs(session, js""" - console.log("Hello"); - iter = 1; - eps = 0.1; - s = -0.8; - gamma = 1.5; - beta = 0.0; - si = 0.0; - R1 = $(R1); - R2 = $(R2); - //setInterval( - """) - - #css("#slider", var"background-image"="linear-gradient(blue, green, blue)") - - scene = scatter(repeat(2randn(n), outer=K), repeat(2randn(n),outer=K), color = fill(:white, n*K), - backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = 0.03, + ms = 0.03 + global scene = scatter(repeat(2randn(n), outer=K), repeat(2randn(n),outer=K), color = fill(:white, n*K), + backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = ms, glowwidth = 0.005, glowcolor = :white, resolution=(600,600), limits = limits, ) @@ -73,62 +72,106 @@ function dom_handler(session, request) splot = scene[end] + scatter!(scene, -R1:0.01:R1, sin.(-R1:0.01:R1), color = RGBA{Float32}(0.5, 0.7, 1.0, 0.8), markersize=ms) + kplot = scene[end] + three, canvas = WGLMakie.three_display(session, scene) js_scene = WGLMakie.to_jsscene(three, scene) mesh = js_scene.getObjectByName(string(objectid(splot))) + mesh2 = js_scene.getObjectByName(string(objectid(kplot))) + + # init javascript + evaljs(session, js""" + console.log("Hello"); + console.log("Hello"); + iter = 1; + eps = 0.1; + s = -0.8; + gamma = 1.5; + beta = 0.0; + si = 0.0; + R1 = $(R1); + R2 = $(R2); + setInterval( + function (){ + function randn_bm() { + var u = 0, v = 0; + while(u === 0) u = Math.random(); //Converting [0,1) to (0,1) + while(v === 0) v = Math.random(); + return Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v ); + } + var mu = 0.2; + var mesh = $(mesh); + var K = $(K); + var n = $(n); + var dt = $(dt); + console.log(iter++); + var sqrtdt = $(sqrtdt); + + k = iter%K; + var positions = mesh.geometry.attributes.offset.array; + var color = mesh.geometry.attributes.color.array; + console.log(color.length); + for ( var i = 0; i < n; i++ ) { + inew = k*2*n + 2*i; + iold = ((K + k - 1)%K)*2*n + 2*i; + positions[inew] = positions[iold] + dt/eps*((1 - positions[iold]*positions[iold])*positions[iold] - positions[iold+1] - s); // x + positions[inew+1] = positions[iold+1] + dt*(-positions[iold+1] + gamma*positions[iold] + beta) + si*sqrtdt*randn_bm(); + color[k*4*n + 4*i] = 1.0; + color[k*4*n + 4*i + 1] = 1.0; + color[k*4*n + 4*i + 2] = 1.0; + color[k*4*n + 4*i + 3] = 1.0; + if (Math.random() < 0.01) + { + positions[inew] = (2*Math.random()-1)*R1; + positions[inew+1] = (2*Math.random()-1)*R2; + } + + } + for ( var k = 0; k < K; k++ ) { + for ( var i = 0; i < n; i++ ) { + color[k*4*n + 4*i + 3] = 0.98*color[k*4*n + 4*i + 3]; + } + } + mesh.geometry.attributes.color.needsUpdate = true; + mesh.geometry.attributes.offset.needsUpdate = true; + + } + , 50); + """) - onjs(session, slider1.value, js"""function (value){ - si = value; - }""") onjs(session, slider2.value, js"""function (value){ beta = value; }""") - onjs(session, button.value, js"""function (value){ - function randn_bm() { - var u = 0, v = 0; - while(u === 0) u = Math.random(); //Converting [0,1) to (0,1) - while(v === 0) v = Math.random(); - return Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v ); - } - var mesh = $(mesh); - var K = $(K); - var n = $(n); - var dt = $(dt); - console.log(iter++); - var sqrtdt = $(sqrtdt); - - k = iter%K; + onjs(session, slider3.value, js"""function (value){ + gamma = value; + }""") + + onjs(session, slider4.value, js"""function (value){ + s = value; + }""") + + onjs(session, slider5.value, js"""function (value){ + eps = value; + }""") + onjs(session, slider1.value, js"""function (value){ + si = value; + var mesh = $(mesh2); var positions = mesh.geometry.attributes.offset.array; var color = mesh.geometry.attributes.color.array; - console.log(color.length); - for ( var i = 0; i < n; i++ ) { - inew = k*2*n + 2*i; - iold = ((K + k - 1)%K)*2*n + 2*i; - positions[inew] = positions[iold] + dt/eps*((1 - positions[iold]*positions[iold])*positions[iold] - positions[iold+1] - s); // x - positions[inew+1] = positions[iold+1] + dt*(-positions[iold+1] + gamma*positions[iold] + beta) + si*sqrtdt*randn_bm(); - color[k*4*n + 4*i] = 1.0; - color[k*4*n + 4*i + 1] = 1.0; - color[k*4*n + 4*i + 2] = 1.0; - color[k*4*n + 4*i + 3] = 1.0; - if (Math.random() < 0.01) - { - positions[inew] = (2*Math.random()-1)*R1; - positions[inew+1] = (2*Math.random()-1)*R2; - } - } - for ( var k = 0; k < K; k++ ) { - for ( var i = 0; i < n; i++ ) { - color[k*4*n + 4*i + 3] = 0.98*color[k*4*n + 4*i + 3]; + for ( var i = 0, l = positions.length; i < l; i += 2 ) { + positions[i+1] = si*Math.sin(positions[i]); } - } - mesh.geometry.attributes.color.needsUpdate = true; mesh.geometry.attributes.offset.needsUpdate = true; + //mesh.geometry.attributes.color.needsUpdate = true; }""") - dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters", DOM.div(slider1, id="slider1"), DOM.div(slider2, id="slider2")), - DOM.p(nrs1), DOM.p(nrs2), DOM.p(button)) + + dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters", DOM.div(slider1, id="slider1"), DOM.div(slider2, id="slider2"), + DOM.div(slider3, id="slider3"),DOM.div(slider4, id="slider4"), DOM.div(slider5, id="slider5")), + DOM.p(nrs1), DOM.p(nrs2), DOM.p(nrs3), DOM.p(nrs4), DOM.p(nrs5), DOM.p(button)) # JSServe.onload(session, dom, js""" # iter = 1; # """) @@ -143,7 +186,7 @@ app = JSServe.Application( parse(Int, get(ENV, "WEBIO_HTTP_PORT", "8081")), verbose = false ) - -#cl() = (close(app), "stopped") +cl() = (close(app), "stopped") #println("Done.") +# #cl() diff --git a/scripts/nextjournal/particle.css b/scripts/nextjournal/particle.css index af01d28..6eafc55 100644 --- a/scripts/nextjournal/particle.css +++ b/scripts/nextjournal/particle.css @@ -9,3 +9,22 @@ width: 4cm; background-image: linear-gradient(to right, #8888FF, #FFFF88, #8888FF); } + + +#slider3 { + background-color: #778899; + width: 4cm; + background-image: linear-gradient(to right, #8888FF, #FFFF88, #8888FF); +} + +#slider4 { + background-color: #778899; + width: 4cm; + background-image: linear-gradient(to right, #8888FF, #FFFF88, #8888FF); +} + +#slider5 { + background-color: #778899; + width: 4cm; + background-image: linear-gradient(to right, #8888FF, #FFFF88, #8888FF); +} From c25d3d4af7c8af4fa5b393a25040320704fc1fee Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Sat, 21 Mar 2020 09:59:58 +0100 Subject: [PATCH 17/29] add visible nullclines and small fixes --- scripts/nextjournal/domserv_fhn.jl | 37 ++++++++++++++++++------------ 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/scripts/nextjournal/domserv_fhn.jl b/scripts/nextjournal/domserv_fhn.jl index bee83c9..f0f0809 100644 --- a/scripts/nextjournal/domserv_fhn.jl +++ b/scripts/nextjournal/domserv_fhn.jl @@ -31,18 +31,18 @@ function dom_handler(session, request) linkjs(session, slider2.value, nrs2.value) # slider and field for gamma - slider3 = JSServe.Slider(0.01:0.01:10) - nrs3 = JSServe.NumberInput(0.0) + slider3 = JSServe.Slider(1.51:0.01:10) + nrs3 = JSServe.NumberInput(1.5) linkjs(session, slider3.value, nrs3.value) # slider and field for s - slider4 = JSServe.Slider(0.01:0.01:10) + slider4 = JSServe.Slider(0.00:0.01:10) nrs4 = JSServe.NumberInput(0.0) linkjs(session, slider4.value, nrs4.value) # slider and field for esp - slider5 = JSServe.Slider(0.01:0.01:10) - nrs5 = JSServe.NumberInput(0.0) + slider5 = JSServe.Slider(0.11:0.01:10) + nrs5 = JSServe.NumberInput(0.1) linkjs(session, slider5.value, nrs5.value) # time wheel ;-) @@ -56,11 +56,16 @@ function dom_handler(session, request) K = 80 dt = 0.001 sqrtdt = sqrt(dt) - + eps = 0.1; + s = -0.0; + gamma = 1.5; + beta = 0.0; + si = 0.0; particlecss = Asset(joinpath(@__DIR__,"particle.css")) - ms = 0.03 + ms1 = 0.04 + ms2 = 0.08 global scene = scatter(repeat(2randn(n), outer=K), repeat(2randn(n),outer=K), color = fill(:white, n*K), - backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = ms, + backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = ms1, glowwidth = 0.005, glowcolor = :white, resolution=(600,600), limits = limits, ) @@ -72,8 +77,10 @@ function dom_handler(session, request) splot = scene[end] - scatter!(scene, -R1:0.01:R1, sin.(-R1:0.01:R1), color = RGBA{Float32}(0.5, 0.7, 1.0, 0.8), markersize=ms) - kplot = scene[end] + scatter!(scene, -R1:0.01:R1, (-R1:0.01:R1) .- (-R1:0.01:R1).^3 + s, color = RGBA{Float32}(255, 0.0, 4.0, 1.0), markersize=ms2) + kplot1 = scene[end] + scatter!(scene, -R1:0.01:R1, gamma*(-R1:0.01:R1) .+ beta , color = RGBA{Float32}(255, 0.0, 4.0, 1.0), markersize=ms2) + kplot2 = scene[end] three, canvas = WGLMakie.three_display(session, scene) js_scene = WGLMakie.to_jsscene(three, scene) @@ -85,11 +92,11 @@ function dom_handler(session, request) console.log("Hello"); console.log("Hello"); iter = 1; - eps = 0.1; - s = -0.8; - gamma = 1.5; - beta = 0.0; - si = 0.0; + eps = $(eps); + s = $(s); + gamma = $(gamma); + beta = $(beta); + si = $(si); R1 = $(R1); R2 = $(R2); setInterval( From 090737a12244a78c2d7bf497adabbe83f22c70b8 Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Sat, 21 Mar 2020 12:23:23 +0100 Subject: [PATCH 18/29] updatekline not working --- scripts/nextjournal/domserv.jl | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/scripts/nextjournal/domserv.jl b/scripts/nextjournal/domserv.jl index dc54b0a..eec8558 100644 --- a/scripts/nextjournal/domserv.jl +++ b/scripts/nextjournal/domserv.jl @@ -52,7 +52,7 @@ function dom_handler(session, request) splot = scene[end] - scatter!(scene, -R1:0.01:R1, sin.(-R1:0.01:R1), color = RGBA{Float32}(0.5, 0.7, 1.0, 0.8), markersize=ms) + scatter!(scene, -R1:0.01:R1, sin.(-R1:0.01:R1), color = RGBA{Float32}(255, 0.0, 4.0, 1.0), markersize=ms) kplot = scene[end] three, canvas = WGLMakie.three_display(session, scene) @@ -67,6 +67,17 @@ function dom_handler(session, request) si = 0.0; R1 = $(R1); R2 = $(R2); + function updatekline(value){ + si = value; + var mesh = $(mesh2); + var positions = mesh.geometry.attributes.offset.array; + var color = mesh.geometry.attributes.color.array; + for ( var i = 0, l = positions.length; i < l; i += 2 ) { + positions[i+1] = si*Math.sin(positions[i]); + } + mesh.geometry.attributes.offset.needsUpdate = true; + //mesh.geometry.attributes.color.needsUpdate = true; + } setInterval( function (){ function randn_bm() { @@ -114,19 +125,8 @@ function dom_handler(session, request) } , 50); """) - onjs(session, sliders.value, js"""function (value){ - si = value; - var mesh = $(mesh2); - var positions = mesh.geometry.attributes.offset.array; - var color = mesh.geometry.attributes.color.array; - - for ( var i = 0, l = positions.length; i < l; i += 2 ) { - positions[i+1] = si*Math.sin(positions[i]); - } - mesh.geometry.attributes.offset.needsUpdate = true; - //mesh.geometry.attributes.color.needsUpdate = true; - + updatekline(value); }""") dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters"), DOM.div(sliders, id="slider"), @@ -146,6 +146,6 @@ app = JSServe.Application( verbose = false ) cl() = (close(app), "stopped") -println("Done.") +#println("Done.") # cl() From 8cbef36e0fe8c8dfe83c8e8ea5d39eb8f3c974ca Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Sat, 21 Mar 2020 12:39:38 +0100 Subject: [PATCH 19/29] samll fix, still nullkline does not move --- scripts/nextjournal/domserv.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/nextjournal/domserv.jl b/scripts/nextjournal/domserv.jl index eec8558..0cdc37d 100644 --- a/scripts/nextjournal/domserv.jl +++ b/scripts/nextjournal/domserv.jl @@ -67,7 +67,7 @@ function dom_handler(session, request) si = 0.0; R1 = $(R1); R2 = $(R2); - function updatekline(value){ + updatekline = function (value){ si = value; var mesh = $(mesh2); var positions = mesh.geometry.attributes.offset.array; From 8e309901432aee1b11a3562073bd4f1532c9d343 Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Sat, 21 Mar 2020 16:25:34 +0100 Subject: [PATCH 20/29] interface fully functioning --- scripts/nextjournal/domserv_fhn.jl | 63 ++++++++++++++++++++---------- 1 file changed, 43 insertions(+), 20 deletions(-) diff --git a/scripts/nextjournal/domserv_fhn.jl b/scripts/nextjournal/domserv_fhn.jl index f0f0809..19f80d3 100644 --- a/scripts/nextjournal/domserv_fhn.jl +++ b/scripts/nextjournal/domserv_fhn.jl @@ -16,6 +16,8 @@ using JSServe: Application, Session, evaljs, linkjs, div, active_sessions, Asset using JSServe: @js_str, onjs, Button, TextField, Slider, JSString, Dependency, with_session using JSServe.DOM +#Import dataframe + function dom_handler(session, request) global three, scene @@ -57,7 +59,7 @@ function dom_handler(session, request) dt = 0.001 sqrtdt = sqrt(dt) eps = 0.1; - s = -0.0; + s = 0.0; gamma = 1.5; beta = 0.0; si = 0.0; @@ -77,7 +79,7 @@ function dom_handler(session, request) splot = scene[end] - scatter!(scene, -R1:0.01:R1, (-R1:0.01:R1) .- (-R1:0.01:R1).^3 + s, color = RGBA{Float32}(255, 0.0, 4.0, 1.0), markersize=ms2) + scatter!(scene, -R1:0.01:R1, (-R1:0.01:R1) .- (-R1:0.01:R1).^3 .+ s, color = RGBA{Float32}(255, 0.0, 4.0, 1.0), markersize=ms2) kplot1 = scene[end] scatter!(scene, -R1:0.01:R1, gamma*(-R1:0.01:R1) .+ beta , color = RGBA{Float32}(255, 0.0, 4.0, 1.0), markersize=ms2) kplot2 = scene[end] @@ -85,8 +87,8 @@ function dom_handler(session, request) three, canvas = WGLMakie.three_display(session, scene) js_scene = WGLMakie.to_jsscene(three, scene) mesh = js_scene.getObjectByName(string(objectid(splot))) - mesh2 = js_scene.getObjectByName(string(objectid(kplot))) - + mesh1 = js_scene.getObjectByName(string(objectid(kplot1))) + mesh2 = js_scene.getObjectByName(string(objectid(kplot2))) # init javascript evaljs(session, js""" console.log("Hello"); @@ -99,6 +101,36 @@ function dom_handler(session, request) si = $(si); R1 = $(R1); R2 = $(R2); + updateklinebeta = function (value){ + beta = value; + var mesh = $(mesh2); + var positions = mesh.geometry.attributes.offset.array; + for ( var i = 0, l = positions.length; i < l; i += 2 ) { + positions[i+1] = beta + positions[i]*gamma; + } + mesh.geometry.attributes.offset.needsUpdate = true; + //mesh.geometry.attributes.color.needsUpdate = true; + } + updateklinegamma = function (value){ + gamma = value; + var mesh = $(mesh2); + var positions = mesh.geometry.attributes.offset.array; + for ( var i = 0, l = positions.length; i < l; i += 2 ) { + positions[i+1] = beta + positions[i]*gamma; + } + mesh.geometry.attributes.offset.needsUpdate = true; + //mesh.geometry.attributes.color.needsUpdate = true; + } + updateklines = function (value){ + s = value; + var mesh = $(mesh1); + var positions = mesh.geometry.attributes.offset.array; + for ( var i = 0, l = positions.length; i < l; i += 2 ) { + positions[i+1] = positions[i] - positions[i]*positions[i]*positions[i] + s; + } + mesh.geometry.attributes.offset.needsUpdate = true; + //mesh.geometry.attributes.color.needsUpdate = true; + } setInterval( function (){ function randn_bm() { @@ -122,7 +154,7 @@ function dom_handler(session, request) for ( var i = 0; i < n; i++ ) { inew = k*2*n + 2*i; iold = ((K + k - 1)%K)*2*n + 2*i; - positions[inew] = positions[iold] + dt/eps*((1 - positions[iold]*positions[iold])*positions[iold] - positions[iold+1] - s); // x + positions[inew] = positions[iold] + dt/eps*((1 - positions[iold]*positions[iold])*positions[iold] - positions[iold+1] + s); // x positions[inew+1] = positions[iold+1] + dt*(-positions[iold+1] + gamma*positions[iold] + beta) + si*sqrtdt*randn_bm(); color[k*4*n + 4*i] = 1.0; color[k*4*n + 4*i + 1] = 1.0; @@ -148,32 +180,22 @@ function dom_handler(session, request) """) onjs(session, slider2.value, js"""function (value){ - beta = value; + updateklinebeta(value); }""") onjs(session, slider3.value, js"""function (value){ - gamma = value; + updateklinegamma(value); }""") onjs(session, slider4.value, js"""function (value){ - s = value; + updateklines(value); }""") onjs(session, slider5.value, js"""function (value){ eps = value; }""") onjs(session, slider1.value, js"""function (value){ - si = value; - var mesh = $(mesh2); - var positions = mesh.geometry.attributes.offset.array; - var color = mesh.geometry.attributes.color.array; - - for ( var i = 0, l = positions.length; i < l; i += 2 ) { - positions[i+1] = si*Math.sin(positions[i]); - } - mesh.geometry.attributes.offset.needsUpdate = true; - //mesh.geometry.attributes.color.needsUpdate = true; - + sigma = value; }""") dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters", DOM.div(slider1, id="slider1"), DOM.div(slider2, id="slider2"), @@ -193,7 +215,8 @@ app = JSServe.Application( parse(Int, get(ENV, "WEBIO_HTTP_PORT", "8081")), verbose = false ) + cl() = (close(app), "stopped") #println("Done.") # -#cl() +cl() From d29254bb142511887c3503ab3da77a78c8676c24 Mon Sep 17 00:00:00 2001 From: mschauer Date: Sat, 21 Mar 2020 17:19:03 +0100 Subject: [PATCH 21/29] Small fixes --- scripts/nextjournal/domserv_fhn.jl | 31 ++++++++++-------------------- 1 file changed, 10 insertions(+), 21 deletions(-) diff --git a/scripts/nextjournal/domserv_fhn.jl b/scripts/nextjournal/domserv_fhn.jl index 19f80d3..890d035 100644 --- a/scripts/nextjournal/domserv_fhn.jl +++ b/scripts/nextjournal/domserv_fhn.jl @@ -2,11 +2,9 @@ using JSServe, WGLMakie, AbstractPlotting using JSServe: JSServe.DOM, @js_str, onjs global three, scene -using StaticArrays using Colors using Random -rebirth(α, R) = x -> (rand() > α ? x : (2rand(typeof(x)) .- 1).*R) -const 𝕏 = SVector +using WGLMakie: scatter, scatter! @@ -42,16 +40,13 @@ function dom_handler(session, request) nrs4 = JSServe.NumberInput(0.0) linkjs(session, slider4.value, nrs4.value) - # slider and field for esp + # slider and field for eps slider5 = JSServe.Slider(0.11:0.01:10) nrs5 = JSServe.NumberInput(0.1) linkjs(session, slider5.value, nrs5.value) - # time wheel ;-) - button = JSServe.Slider(1:109) - # init - R = 𝕏(1.5,6.0) + R = (1.5,6.0) R1, R2 = R limits = FRect(-R[1], -R[2], 2R[1], 2R[2]) n = 800 @@ -65,7 +60,7 @@ function dom_handler(session, request) si = 0.0; particlecss = Asset(joinpath(@__DIR__,"particle.css")) ms1 = 0.04 - ms2 = 0.08 + ms2 = 0.04 global scene = scatter(repeat(2randn(n), outer=K), repeat(2randn(n),outer=K), color = fill(:white, n*K), backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = ms1, glowwidth = 0.005, glowcolor = :white, @@ -79,9 +74,9 @@ function dom_handler(session, request) splot = scene[end] - scatter!(scene, -R1:0.01:R1, (-R1:0.01:R1) .- (-R1:0.01:R1).^3 .+ s, color = RGBA{Float32}(255, 0.0, 4.0, 1.0), markersize=ms2) + scatter!(scene, -R1:0.01:R1, (-R1:0.01:R1) .- (-R1:0.01:R1).^3 .+ s, color = RGBA{Float32}(0.5, 0.7, 1.0, 0.8), markersize=ms2) kplot1 = scene[end] - scatter!(scene, -R1:0.01:R1, gamma*(-R1:0.01:R1) .+ beta , color = RGBA{Float32}(255, 0.0, 4.0, 1.0), markersize=ms2) + scatter!(scene, -R1:0.01:R1, gamma*(-R1:0.01:R1) .+ beta , color = RGBA{Float32}(0.5, 0.7, 1.0, 0.8), markersize=ms2) kplot2 = scene[end] three, canvas = WGLMakie.three_display(session, scene) @@ -146,7 +141,6 @@ function dom_handler(session, request) var dt = $(dt); console.log(iter++); var sqrtdt = $(sqrtdt); - k = iter%K; var positions = mesh.geometry.attributes.offset.array; var color = mesh.geometry.attributes.color.array; @@ -165,7 +159,6 @@ function dom_handler(session, request) positions[inew] = (2*Math.random()-1)*R1; positions[inew+1] = (2*Math.random()-1)*R2; } - } for ( var k = 0; k < K; k++ ) { for ( var i = 0; i < n; i++ ) { @@ -174,7 +167,6 @@ function dom_handler(session, request) } mesh.geometry.attributes.color.needsUpdate = true; mesh.geometry.attributes.offset.needsUpdate = true; - } , 50); """) @@ -195,15 +187,12 @@ function dom_handler(session, request) eps = value; }""") onjs(session, slider1.value, js"""function (value){ - sigma = value; + si = value; }""") - dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters", DOM.div(slider1, id="slider1"), DOM.div(slider2, id="slider2"), - DOM.div(slider3, id="slider3"),DOM.div(slider4, id="slider4"), DOM.div(slider5, id="slider5")), - DOM.p(nrs1), DOM.p(nrs2), DOM.p(nrs3), DOM.p(nrs4), DOM.p(nrs5), DOM.p(button)) -# JSServe.onload(session, dom, js""" -# iter = 1; -# """) + dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters"), "sigma", DOM.div(slider1, id="slider1"), DOM.p(nrs1), "beta", DOM.div(slider2, id="slider2"), DOM.p(nrs2), "gamma", + DOM.div(slider3, id="slider3"), DOM.p(nrs3), "s", DOM.div(slider4, id="slider4"), DOM.p(nrs4), "eps", DOM.div(slider5, id="slider5"), DOM.p(nrs5)) + println("running...") dom end From bc20498ba6c5031a3b0e65e85ccf9ca4fec95163 Mon Sep 17 00:00:00 2001 From: mschauer Date: Sun, 22 Mar 2020 14:18:00 +0100 Subject: [PATCH 22/29] Bisschen huebscher --- scripts/nextjournal/domserv_fhn.jl | 51 +++++++++++++++++------------- 1 file changed, 29 insertions(+), 22 deletions(-) diff --git a/scripts/nextjournal/domserv_fhn.jl b/scripts/nextjournal/domserv_fhn.jl index 890d035..8cbddf7 100644 --- a/scripts/nextjournal/domserv_fhn.jl +++ b/scripts/nextjournal/domserv_fhn.jl @@ -20,48 +20,50 @@ using JSServe.DOM function dom_handler(session, request) global three, scene + eps = 0.1; + s = 0.0; + gamma = 1.5; + beta = 0.0; + si = 0.0; + # slider and field for sigma - slider1 = JSServe.Slider(0.01:0.01:1) + slider1 = JSServe.Slider(0.01:0.01:1, si) nrs1 = JSServe.NumberInput(0.0) linkjs(session, slider1.value, nrs1.value) # slider and field for beta - slider2 = JSServe.Slider(0.01:0.01:10) + slider2 = JSServe.Slider(0.01:0.01:10, beta) nrs2 = JSServe.NumberInput(0.0) linkjs(session, slider2.value, nrs2.value) # slider and field for gamma - slider3 = JSServe.Slider(1.51:0.01:10) + slider3 = JSServe.Slider(1.51:0.01:10, gamma) nrs3 = JSServe.NumberInput(1.5) linkjs(session, slider3.value, nrs3.value) # slider and field for s - slider4 = JSServe.Slider(0.00:0.01:10) + slider4 = JSServe.Slider(0.00:0.01:10, s) nrs4 = JSServe.NumberInput(0.0) linkjs(session, slider4.value, nrs4.value) # slider and field for eps - slider5 = JSServe.Slider(0.11:0.01:10) + slider5 = JSServe.Slider(0.11:0.01:10, eps) nrs5 = JSServe.NumberInput(0.1) linkjs(session, slider5.value, nrs5.value) # init - R = (1.5,6.0) + R = (1.5,3.0) R1, R2 = R limits = FRect(-R[1], -R[2], 2R[1], 2R[2]) - n = 800 - K = 80 - dt = 0.001 + n = 500 + K = 200 + dt = 0.0005 sqrtdt = sqrt(dt) - eps = 0.1; - s = 0.0; - gamma = 1.5; - beta = 0.0; - si = 0.0; + particlecss = Asset(joinpath(@__DIR__,"particle.css")) - ms1 = 0.04 - ms2 = 0.04 - global scene = scatter(repeat(2randn(n), outer=K), repeat(2randn(n),outer=K), color = fill(:white, n*K), + ms1 = 0.02 + ms2 = 0.02 + global scene = scatter(repeat(R1*(2rand(n) .- 1), outer=K), repeat(R2*(2rand(n) .- 1),outer=K), color = fill((:white,0f0), n*K), backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = ms1, glowwidth = 0.005, glowcolor = :white, resolution=(600,600), limits = limits, @@ -154,7 +156,7 @@ function dom_handler(session, request) color[k*4*n + 4*i + 1] = 1.0; color[k*4*n + 4*i + 2] = 1.0; color[k*4*n + 4*i + 3] = 1.0; - if (Math.random() < 0.01) + if (Math.random() < 0.001) { positions[inew] = (2*Math.random()-1)*R1; positions[inew+1] = (2*Math.random()-1)*R2; @@ -190,14 +192,19 @@ function dom_handler(session, request) si = value; }""") - dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters"), "sigma", DOM.div(slider1, id="slider1"), DOM.p(nrs1), "beta", DOM.div(slider2, id="slider2"), DOM.p(nrs2), "gamma", - DOM.div(slider3, id="slider3"), DOM.p(nrs3), "s", DOM.div(slider4, id="slider4"), DOM.p(nrs4), "eps", DOM.div(slider5, id="slider5"), DOM.p(nrs5)) - + dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters"), DOM.table( + DOM.tr(DOM.td("sigma"), DOM.td("beta"), DOM.td("gamma"), DOM.td("s"), DOM.td("eps")), + DOM.tr( + DOM.td(DOM.div(slider1, id="slider1"), DOM.div(nrs1)), + DOM.td(DOM.div(slider2, id="slider2"), DOM.div(nrs2)), + DOM.td(DOM.div(slider3, id="slider3"), DOM.div(nrs3)), + DOM.td(DOM.div(slider4, id="slider4"), DOM.div(nrs4)), + DOM.td(DOM.div(slider5, id="slider5"), DOM.div(nrs5)) + ))) println("running...") dom end - app = JSServe.Application( dom_handler, get(ENV, "WEBIO_SERVER_HOST_URL", "127.0.0.1"), From ab6f0fa049d5d6a1dbe163f2a5fc36e0287c9ee1 Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Sun, 22 Mar 2020 16:19:53 +0100 Subject: [PATCH 23/29] small changes --- src/BridgeSDEInference.jl | 4 ++++ src/examples/sir.jl | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/BridgeSDEInference.jl b/src/BridgeSDEInference.jl index 8b9a3c8..c171080 100644 --- a/src/BridgeSDEInference.jl +++ b/src/BridgeSDEInference.jl @@ -5,6 +5,10 @@ using Statistics, Random, LinearAlgebra using ForwardDiff using ForwardDiff: value +#sir.jl +export SIR, SIRAux + + # fitzHughNagumo.jl export FitzhughDiffusion, FitzhughDiffusionAux, ℝ export regularToAlter, alterToRegular, regularToConjug, conjugToRegular, display diff --git a/src/examples/sir.jl b/src/examples/sir.jl index 5cefbff..926483c 100644 --- a/src/examples/sir.jl +++ b/src/examples/sir.jl @@ -21,7 +21,7 @@ constdiff(::SIR) = false clone(P::SIR, θ) = SIR(θ...) params(P::SIR) = [P.α, P.β, P.σ1, P.σ2] #domain(P::SIR) = LowerBoundedDomain((0.0, 0.0), (1,2)) -domain(P::SIR) = LowerBoundedDomain((0.0,), (1,)) +domain(P::SIR) = LowerBoundedDomain((0.0, 0.0), (1, 2)) # <--------------------------------------------- From 91e064d6da875e3f298ea7be9ccb4e6c19a6267c Mon Sep 17 00:00:00 2001 From: mschauer Date: Mon, 23 Mar 2020 16:53:24 +0100 Subject: [PATCH 24/29] Serve a background picture for the sliders --- scripts/nextjournal/domserv.jl | 34 ++++++++++++++++++++++---------- scripts/nextjournal/slider1.png | Bin 0 -> 274 bytes 2 files changed, 24 insertions(+), 10 deletions(-) create mode 100644 scripts/nextjournal/slider1.png diff --git a/scripts/nextjournal/domserv.jl b/scripts/nextjournal/domserv.jl index 0cdc37d..3af921b 100644 --- a/scripts/nextjournal/domserv.jl +++ b/scripts/nextjournal/domserv.jl @@ -8,7 +8,9 @@ using Random rebirth(α, R) = x -> (rand() > α ? x : (2rand(typeof(x)) .- 1).*R) const 𝕏 = SVector - +using FileIO +using Makie: band +using GLMakie using Hyperscript, Markdown using JSServe, Observables @@ -16,10 +18,16 @@ using JSServe: Application, Session, evaljs, linkjs, div, active_sessions, Asset using JSServe: @js_str, onjs, Button, TextField, Slider, JSString, Dependency, with_session using JSServe.DOM - function dom_handler(session, request) global three, scene + if @isdefined dom + dom != nothing && return dom + end + + WGLMakie.activate!() + + # slider and field for sigma sliders = JSServe.Slider(0.01:0.01:1) nrs = JSServe.NumberInput(0.0) @@ -38,6 +46,8 @@ function dom_handler(session, request) sqrtdt = sqrt(dt) particlecss = Asset(joinpath(@__DIR__,"particle.css")) + global sliderbg = Asset(joinpath(@__DIR__,"slider1.png")) + ms = 0.03 global scene = scatter(repeat(2randn(n), outer=K), repeat(2randn(n),outer=K), color = fill(:white, n*K), backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = ms, @@ -71,12 +81,10 @@ function dom_handler(session, request) si = value; var mesh = $(mesh2); var positions = mesh.geometry.attributes.offset.array; - var color = mesh.geometry.attributes.color.array; for ( var i = 0, l = positions.length; i < l; i += 2 ) { positions[i+1] = si*Math.sin(positions[i]); } mesh.geometry.attributes.offset.needsUpdate = true; - //mesh.geometry.attributes.color.needsUpdate = true; } setInterval( function (){ @@ -128,12 +136,11 @@ function dom_handler(session, request) onjs(session, sliders.value, js"""function (value){ updatekline(value); }""") + sliderbgurl = JSServe.url(sliderbg) + style_obs = Observable(" background-repeat: no-repeat; background-image: url($sliderbgurl);") - dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters"), DOM.div(sliders, id="slider"), + global dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters"), DOM.div(sliders, id="slider", style=style_obs), DOM.p(nrs)) -# JSServe.onload(session, dom, js""" -# iter = 1; -# """) println("running...") dom end @@ -145,7 +152,14 @@ app = JSServe.Application( parse(Int, get(ENV, "WEBIO_HTTP_PORT", "8081")), verbose = false ) -cl() = (close(app), "stopped") -#println("Done.") +function cl() + close(app) + global dom = nothing +end +println("Done.") # +if false + cl() + +end diff --git a/scripts/nextjournal/slider1.png b/scripts/nextjournal/slider1.png new file mode 100644 index 0000000000000000000000000000000000000000..cf625f33f62b0138b5f05de2417dc9d1bc6dc321 GIT binary patch literal 274 zcmV+t0qy>YP)Px#%t=H+RA>e5)?EsKAP`2;igx7wm*ASr@MTB@!@uu>Jxx+E%yA~q5W=x9?FwAl z>Fw;^3T!06KG;YAq#KZ9Ws2T^A9U`!!n{B62Z!7ymw;=Lk!pA>7?z5O^p^)t1w1&k z2t}>N*RA)g7}ZIA>&uyboKEU1glZ`WQg@*=sRn(v_ZXBPcXR5fJ^F-^=J42EuS%sl z?&b>#|E2UgX$g<*rx=U1loFBxkJXZ)5 Date: Wed, 25 Mar 2020 17:21:38 +0100 Subject: [PATCH 25/29] Current nextjournal code --- scripts/nextjournal/domserv_fhn.jl | 155 ++++++++++++++++++----------- 1 file changed, 98 insertions(+), 57 deletions(-) diff --git a/scripts/nextjournal/domserv_fhn.jl b/scripts/nextjournal/domserv_fhn.jl index 8cbddf7..cd2e203 100644 --- a/scripts/nextjournal/domserv_fhn.jl +++ b/scripts/nextjournal/domserv_fhn.jl @@ -6,7 +6,14 @@ using Colors using Random using WGLMakie: scatter, scatter! - +# fallback values if statistical analysis is not run +if !@isdefined thetas + thetas = [0.1, -0.8, 1.5, 0.0, 0.3] +end +if !@isdefined thetalims + thetalims = [(0.0,1.0), (-5.0,5.0), (0.0,10.0), (0.0,10.0), (0.0,1.0)] +end +WGLMakie.activate!() using Hyperscript, Markdown using JSServe, Observables @@ -14,83 +21,97 @@ using JSServe: Application, Session, evaljs, linkjs, div, active_sessions, Asset using JSServe: @js_str, onjs, Button, TextField, Slider, JSString, Dependency, with_session using JSServe.DOM -#Import dataframe - - function dom_handler(session, request) global three, scene - eps = 0.1; - s = 0.0; - gamma = 1.5; - beta = 0.0; - si = 0.0; + # fetch initial parameter (initial slider settings) + eps = thetas[1]; + s = thetas[2]; + gamma = thetas[3]; + beta = thetas[4]; + si = thetas[5]; + rebirth = 0.001; # how often a particle is "reborn" at random position + sl = 101 # slider sub-divisions + + # load histogram images to use as slider background + sliderbg = [JSServe.Asset(), JSServe.Asset(), JSServe.Asset(), JSServe.Asset(), JSServe.Asset()] # slider and field for sigma - slider1 = JSServe.Slider(0.01:0.01:1, si) - nrs1 = JSServe.NumberInput(0.0) - linkjs(session, slider1.value, nrs1.value) + slider5 = JSServe.Slider(range(thetalims[5]..., length=sl), si) + nrs5 = JSServe.NumberInput(si) + linkjs(session, slider5.value, nrs5.value) # slider and field for beta - slider2 = JSServe.Slider(0.01:0.01:10, beta) - nrs2 = JSServe.NumberInput(0.0) - linkjs(session, slider2.value, nrs2.value) + slider4 = JSServe.Slider(range(thetalims[4]..., length=sl), beta) + nrs4 = JSServe.NumberInput(beta) + linkjs(session, slider4.value, nrs4.value) # slider and field for gamma - slider3 = JSServe.Slider(1.51:0.01:10, gamma) - nrs3 = JSServe.NumberInput(1.5) + slider3 = JSServe.Slider(range(thetalims[3]..., length=sl), gamma) + nrs3 = JSServe.NumberInput(gamma) linkjs(session, slider3.value, nrs3.value) # slider and field for s - slider4 = JSServe.Slider(0.00:0.01:10, s) - nrs4 = JSServe.NumberInput(0.0) - linkjs(session, slider4.value, nrs4.value) + slider2 = JSServe.Slider(range(thetalims[2]..., length=sl), s) + nrs2 = JSServe.NumberInput(s) + linkjs(session, slider2.value, nrs2.value) # slider and field for eps - slider5 = JSServe.Slider(0.11:0.01:10, eps) - nrs5 = JSServe.NumberInput(0.1) - linkjs(session, slider5.value, nrs5.value) + slider1 = JSServe.Slider(range(thetalims[1]..., length=sl), eps) + nrs1 = JSServe.NumberInput(eps) + linkjs(session, slider1.value, nrs1.value) + + # slider and field for rebirth + slider6 = JSServe.Slider(0.0:0.0001:0.005, rebirth) + nrs6 = JSServe.NumberInput(rebirth) + linkjs(session, slider6.value, nrs6.value) # init - R = (1.5,3.0) + R = (1.5, 3.0) # plot area R1, R2 = R limits = FRect(-R[1], -R[2], 2R[1], 2R[2]) - n = 500 - K = 200 - dt = 0.0005 + n = 400 # no of particles + K = 150 # display K past positions of particle fading out + dt = 0.0005 # time step sqrtdt = sqrt(dt) - particlecss = Asset(joinpath(@__DIR__,"particle.css")) - ms1 = 0.02 - ms2 = 0.02 - global scene = scatter(repeat(R1*(2rand(n) .- 1), outer=K), repeat(R2*(2rand(n) .- 1),outer=K), color = fill((:white,0f0), n*K), + particlecss = JSServe.Asset() # style sheet + ms1 = 0.02 # markersize particles + ms2 = 0.02 # markersize isokline + + # plot particles, initially at random positions + global scene = WGLMakie.scatter(repeat(R1*(2rand(n) .- 1), outer=K), repeat(R2*(2rand(n) .- 1),outer=K), color = fill((:white,0f0), n*K), backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = ms1, glowwidth = 0.005, glowcolor = :white, - resolution=(600,600), limits = limits, + resolution=(500,500), limits = limits, ) + + # style plot axis = scene[Axis] axis[:grid, :linewidth] = (0.3, 0.3) axis[:grid, :linecolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.3),RGBA{Float32}(0.5, 0.7, 1.0, 0.3)) axis[:names][:textsize] = (0.0,0.0) axis[:ticks, :textcolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.5),RGBA{Float32}(0.5, 0.7, 1.0, 0.5)) - - splot = scene[end] - scatter!(scene, -R1:0.01:R1, (-R1:0.01:R1) .- (-R1:0.01:R1).^3 .+ s, color = RGBA{Float32}(0.5, 0.7, 1.0, 0.8), markersize=ms2) + + # plot isoklines + WGLMakie.scatter!(scene, -R1:0.01:R1, (-R1:0.01:R1) .- (-R1:0.01:R1).^3 .+ s, color = RGBA{Float32}(0.5, 0.7, 1.0, 0.8), markersize=ms2) kplot1 = scene[end] - scatter!(scene, -R1:0.01:R1, gamma*(-R1:0.01:R1) .+ beta , color = RGBA{Float32}(0.5, 0.7, 1.0, 0.8), markersize=ms2) + WGLMakie.scatter!(scene, -R1:0.01:R1, gamma*(-R1:0.01:R1) .+ beta , color = RGBA{Float32}(0.5, 0.7, 1.0, 0.8), markersize=ms2) kplot2 = scene[end] + # set up threejs scene three, canvas = WGLMakie.three_display(session, scene) js_scene = WGLMakie.to_jsscene(three, scene) mesh = js_scene.getObjectByName(string(objectid(splot))) mesh1 = js_scene.getObjectByName(string(objectid(kplot1))) mesh2 = js_scene.getObjectByName(string(objectid(kplot2))) + # init javascript evaljs(session, js""" - console.log("Hello"); - console.log("Hello"); - iter = 1; + iter = 1; // iteration number + + // fetch parameters eps = $(eps); s = $(s); gamma = $(gamma); @@ -98,6 +119,9 @@ function dom_handler(session, request) si = $(si); R1 = $(R1); R2 = $(R2); + rebirth = $(rebirth); + + // update functions for isoklines updateklinebeta = function (value){ beta = value; var mesh = $(mesh2); @@ -128,6 +152,8 @@ function dom_handler(session, request) mesh.geometry.attributes.offset.needsUpdate = true; //mesh.geometry.attributes.color.needsUpdate = true; } + + // move particles every x milliseconds setInterval( function (){ function randn_bm() { @@ -156,7 +182,7 @@ function dom_handler(session, request) color[k*4*n + 4*i + 1] = 1.0; color[k*4*n + 4*i + 2] = 1.0; color[k*4*n + 4*i + 3] = 1.0; - if (Math.random() < 0.001) + if (Math.random() < rebirth) { positions[inew] = (2*Math.random()-1)*R1; positions[inew+1] = (2*Math.random()-1)*R2; @@ -170,41 +196,56 @@ function dom_handler(session, request) mesh.geometry.attributes.color.needsUpdate = true; mesh.geometry.attributes.offset.needsUpdate = true; } - , 50); + , 15); // milliseconds """) + # react on slider movements + + onjs(session, slider1.value, js"""function (value){ + eps = value; + }""") onjs(session, slider2.value, js"""function (value){ - updateklinebeta(value); + updateklines(value); }""") - - onjs(session, slider3.value, js"""function (value){ + onjs(session, slider3.value, js"""function (value){ updateklinegamma(value); }""") - onjs(session, slider4.value, js"""function (value){ - updateklines(value); + updateklinebeta(value); }""") - onjs(session, slider5.value, js"""function (value){ - eps = value; - }""") - onjs(session, slider1.value, js"""function (value){ si = value; }""") + onjs(session, slider6.value, js"""function (value){ + rebirth = value; + }""") + + # set background for sliders + styles = [Observable("padding-left:0px; padding-top: 10px; height: 50px; width=150px; background-size: 115px 50px; background-repeat: no-repeat; + background-position: center center; background-image: url($(JSServe.url(sliderbg[j])));") for j in 1:5] + # arrange canvas and sliders as html elements dom = DOM.div(particlecss, DOM.p(canvas), DOM.p("Parameters"), DOM.table( - DOM.tr(DOM.td("sigma"), DOM.td("beta"), DOM.td("gamma"), DOM.td("s"), DOM.td("eps")), + DOM.tr(DOM.td("eps"), DOM.td("s"), DOM.td("gamma"), DOM.td("beta"), DOM.td("sigma")), DOM.tr( - DOM.td(DOM.div(slider1, id="slider1"), DOM.div(nrs1)), - DOM.td(DOM.div(slider2, id="slider2"), DOM.div(nrs2)), - DOM.td(DOM.div(slider3, id="slider3"), DOM.div(nrs3)), - DOM.td(DOM.div(slider4, id="slider4"), DOM.div(nrs4)), - DOM.td(DOM.div(slider5, id="slider5"), DOM.div(nrs5)) - ))) + DOM.td(DOM.div(slider1, id="slider1", style=styles[1]), DOM.div(nrs1)), + DOM.td(DOM.div(slider2, id="slider2", style=styles[2]), DOM.div(nrs2)), + DOM.td(DOM.div(slider3, id="slider3", style=styles[3]), DOM.div(nrs3)), + DOM.td(DOM.div(slider4, id="slider4", style=styles[4]), DOM.div(nrs4)), + DOM.td(DOM.div(slider5, id="slider5", style=styles[5] ), DOM.div(nrs5)) + ), + DOM.tr( + DOM.td("rebirth", DOM.div(slider6, id="slider6"), DOM.div(nrs6)), + ))) println("running...") dom end +# attach handler to current session +#JSServe.with_session() do session, request +# dom_handler(session, request) +#end + app = JSServe.Application( dom_handler, get(ENV, "WEBIO_SERVER_HOST_URL", "127.0.0.1"), From 6da3c8f1a1a0b0315dc84cabfaa893a9bc6b5d9f Mon Sep 17 00:00:00 2001 From: SebaGraz Date: Wed, 25 Mar 2020 22:59:43 +0100 Subject: [PATCH 26/29] Script to profile --- scripts/nextjournal/nextjournal_script.jl | 85 +++++++++++++++++++++++ 1 file changed, 85 insertions(+) create mode 100644 scripts/nextjournal/nextjournal_script.jl diff --git a/scripts/nextjournal/nextjournal_script.jl b/scripts/nextjournal/nextjournal_script.jl new file mode 100644 index 0000000..d291e51 --- /dev/null +++ b/scripts/nextjournal/nextjournal_script.jl @@ -0,0 +1,85 @@ +## Intro +using Bridge +using StaticArrays +using BridgeSDEInference +using Random, LinearAlgebra, Distributions +const State = SArray{Tuple{2},T,1,2} where {T}; + +param = :regular # regular parametrization of the trajectories (no transformation of the coordinate processes) +ε = 0.1 ; s =-0.8 ; γ =1.5 ; β = 0.0 ; σ =0.3 ; +P = FitzhughDiffusion(param, ε, s, γ, β, σ); + +include(joinpath(pathof(BridgeSDEInference), "..", "auxiliary", "data_simulation_fns.jl")) +# starting point under :regular parametrisation +x0 = State(-0.5, -0.6) + +## Simulation +# time grid +dt = 1/50000 +T = 20.0 +tt = 0.0:dt:T + +Random.seed!(4) +X, _ = simulate_segment(0.0, x0, P, tt); + +# subsampling +num_obs = 100 +skip = div(length(tt), num_obs) +Σ = [10^(-4)] +L = [1.0 0.0] +obs = (time = X.tt[1:skip:end], +values = [Float64(rand(MvNormal(L*x, Σ))[1]) for x in X.yy[1:skip:end]]); + +## Plotting +using Plots +gr() +Plots.plot(X.tt, first.(X.yy), label = "X") +Plots.plot!(X.tt, last.(X.yy), label = "Y") +Plots.scatter!(obs.time, obs.values, markersize=1.5, label = "observations") + +## Inference +θ_init = [ε, s, γ, β, σ].*(1 .+ randn(5).*[0.5, 0.5, 0.5, 0.0, 0.5]) +# Verify that the parameters lay in their domain +θ_init[1]<0 ? θ_init[1] = -θ_init[1] : θ_init[1] = θ_init[1] +θ_init[3]<0 ? θ_init[3] = -θ_init[3] : θ_init[3] = θ_init[3] +θ_init[5]<0 ? θ_init[5] = -θ_init[5] : θ_init[5] = θ_init[5]; +# Take the real β, as it is fixed. + +P_trgt = FitzhughDiffusion(param, θ_init...) +P_aux = [FitzhughDiffusionAux(param, θ_init..., t₀, u, T, v) for (t₀,T,u,v) + in zip(obs.time[1:end-1], obs.time[2:end], obs.values[1:end-1], obs.values[2:end])] + +# Container +model_setup = DiffusionSetup(P_trgt, P_aux, PartObs()); + +initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt(100)) +# Further setting +set_auxiliary!(model_setup; skip_for_save=1, adaptive_prop=NoAdaptation()); + +mcmc_setup = MCMCSetup( + Imputation(NoBlocking(), 0.975, Vern7()), + ParamUpdate(MetropolisHastingsUpdt(), 1, θ_init, + UniformRandomWalk(0.5, true), ImproperPosPrior(), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P_aux, 1)) + ), + ParamUpdate(MetropolisHastingsUpdt(), 2, θ_init, + UniformRandomWalk(0.5, false), ImproperPrior(), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P_aux, 2)) + ), + ParamUpdate(MetropolisHastingsUpdt(), 3, θ_init, + UniformRandomWalk(0.5, true), ImproperPosPrior(), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P_aux, 3)) + ), + ParamUpdate(MetropolisHastingsUpdt(), 5, θ_init, + UniformRandomWalk(0.5, true), ImproperPosPrior(), + UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P_aux, 5)) + )) + +schedule = MCMCSchedule(1*10^3, [[1,2,3,4,5]], + (save=3*10^2, verbose=10^2, warm_up=100, + readjust=(x->x%100==0), fuse=(x->false))); + +Random.seed!(4) +out = mcmc(mcmc_setup, schedule, model_setup); + + From ea77f7a21b5158cd9d67d365b41b74246324a7fc Mon Sep 17 00:00:00 2001 From: mschauer Date: Thu, 26 Mar 2020 12:08:10 +0100 Subject: [PATCH 27/29] Fixes --- scripts/nextjournal/nextjournal_script.jl | 38 ++++++++++++++--------- 1 file changed, 24 insertions(+), 14 deletions(-) diff --git a/scripts/nextjournal/nextjournal_script.jl b/scripts/nextjournal/nextjournal_script.jl index d291e51..95046cc 100644 --- a/scripts/nextjournal/nextjournal_script.jl +++ b/scripts/nextjournal/nextjournal_script.jl @@ -1,10 +1,11 @@ ## Intro +using Profile, Revise using Bridge using StaticArrays using BridgeSDEInference using Random, LinearAlgebra, Distributions const State = SArray{Tuple{2},T,1,2} where {T}; - +using BridgeSDEInference: EulerMaruyamaBounded param = :regular # regular parametrization of the trajectories (no transformation of the coordinate processes) ε = 0.1 ; s =-0.8 ; γ =1.5 ; β = 0.0 ; σ =0.3 ; P = FitzhughDiffusion(param, ε, s, γ, β, σ); @@ -27,22 +28,19 @@ num_obs = 100 skip = div(length(tt), num_obs) Σ = [10^(-4)] L = [1.0 0.0] -obs = (time = X.tt[1:skip:end], -values = [Float64(rand(MvNormal(L*x, Σ))[1]) for x in X.yy[1:skip:end]]); +obs = (time = X.tt[1:skip:end], values = [Float64(rand(MvNormal(L*x, Σ))[1]) for x in X.yy[1:skip:end]]); ## Plotting -using Plots +#=using Plots gr() Plots.plot(X.tt, first.(X.yy), label = "X") Plots.plot!(X.tt, last.(X.yy), label = "Y") Plots.scatter!(obs.time, obs.values, markersize=1.5, label = "observations") - +=# ## Inference -θ_init = [ε, s, γ, β, σ].*(1 .+ randn(5).*[0.5, 0.5, 0.5, 0.0, 0.5]) -# Verify that the parameters lay in their domain -θ_init[1]<0 ? θ_init[1] = -θ_init[1] : θ_init[1] = θ_init[1] -θ_init[3]<0 ? θ_init[3] = -θ_init[3] : θ_init[3] = θ_init[3] -θ_init[5]<0 ? θ_init[5] = -θ_init[5] : θ_init[5] = θ_init[5]; +θ_init = [ε, s, γ, β, σ].*(1 .+ (2*rand(5) .- 1).*[0.2, 0.2, 0.2, 0.0, 0.2]) + + # Take the real β, as it is fixed. P_trgt = FitzhughDiffusion(param, θ_init...) @@ -50,7 +48,20 @@ P_aux = [FitzhughDiffusionAux(param, θ_init..., t₀, u, T, v) for (t₀,T,u,v) in zip(obs.time[1:end-1], obs.time[2:end], obs.values[1:end-1], obs.values[2:end])] # Container -model_setup = DiffusionSetup(P_trgt, P_aux, PartObs()); +model_setup = DiffusionSetup(P_trgt, P_aux, PartObs()) + +# Observation scheme +L = @SMatrix [1. 0.] +Σ = @SMatrix [10^(-4)] +set_observations!(model_setup, [L for _ in P_aux], [Σ for _ in P_aux], + obs.values, obs.time) + +# Imputation grid +dt = 1/200 +set_imputation_grid!(model_setup, dt) +# Prior distribution on (X_0, Y_0) +set_x0_prior!(model_setup, KnownStartingPt(x0)); +#set_x0_prior!(model_setup, GsnStartingPt(x0, @SMatrix [.1 0; 0 .1]), x0); initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt(100)) # Further setting @@ -80,6 +91,5 @@ schedule = MCMCSchedule(1*10^3, [[1,2,3,4,5]], readjust=(x->x%100==0), fuse=(x->false))); Random.seed!(4) -out = mcmc(mcmc_setup, schedule, model_setup); - - +Profile.init() +@profile out = mcmc(mcmc_setup, schedule, model_setup); From dfab57f9ad58f702e309c01d5571f13cba449aed Mon Sep 17 00:00:00 2001 From: mschauer Date: Thu, 26 Mar 2020 12:42:19 +0100 Subject: [PATCH 28/29] Fix some type inference issues --- src/mcmc/mcmc.jl | 4 +++- src/mcmc/updates.jl | 8 +++++++- src/mcmc/workspace.jl | 18 +++++++++++------- 3 files changed, 21 insertions(+), 9 deletions(-) diff --git a/src/mcmc/mcmc.jl b/src/mcmc/mcmc.jl index 6775cc7..0507a2f 100644 --- a/src/mcmc/mcmc.jl +++ b/src/mcmc/mcmc.jl @@ -13,7 +13,9 @@ function mcmc(setup_mcmc::MCMCSetup, schedule::MCMCSchedule, setup::T) where T < ws, ll, θ = create_workspace(setup) ws_mcmc = create_workspace(setup_mcmc, schedule, θ) adpt = adaptation_object(setup, ws) - + mcmc_(setup_mcmc, schedule, setup, ws, ll, θ, ws_mcmc, adpt) +end +function mcmc_(setup_mcmc, schedule, setup, ws, ll, θ, ws_mcmc, adpt) aux = nothing for step in schedule step.save && save_imputed!(ws) diff --git a/src/mcmc/updates.jl b/src/mcmc/updates.jl index 6cb1079..cd69813 100644 --- a/src/mcmc/updates.jl +++ b/src/mcmc/updates.jl @@ -138,7 +138,13 @@ end Preconditioned Crank-Nicolson update with memory parameter `ρ`, previous vector `y` and new vector `yᵒ` """ -crank_nicolson!(yᵒ, y, ρ) = (yᵒ .= √(1-ρ)*yᵒ + √(ρ)*y) +function crank_nicolson!(yᵒ, y, ρ) + r1 = √(1-ρ) + r2 = √(ρ) + for i in eachindex(y) + yᵒ[i] = r1*yᵒ[i] + r2*y[i] + end +end """ path_log_likhd(::OS, XX, P, iRange, fpt; skipFPT=false diff --git a/src/mcmc/workspace.jl b/src/mcmc/workspace.jl index 85391ff..920679a 100644 --- a/src/mcmc/workspace.jl +++ b/src/mcmc/workspace.jl @@ -151,7 +151,7 @@ set!(x::SingleElem{T}, y::T) where T = (x.val = y) The main container of the `mcmc` function from `mcmc.jl` in which most data pertinent to sampling is stored """ -struct Workspace{ObsScheme,S,TX,TW,R,TP,TZ}# ,Q, where Q = eltype(result) +struct Workspace{ObsScheme,S,TX,TW,R,TP,TL,TZ}# ,Q, where Q = eltype(result) # Related to imputed path Wnr::Wiener{S} # Wiener, driving law XXᵒ::Vector{TX} # Diffusion proposal paths @@ -167,6 +167,7 @@ struct Workspace{ObsScheme,S,TX,TW,R,TP,TZ}# ,Q, where Q = eltype(result) time::Vector{Float64} # Storage with time axis # Related to the starting point x0_prior::TP + ll::TL z::SingleElem{TZ} #recompute_ODEs::Vector{Bool} # Info on whether to recompute H,Hν,c after resp. param updt @@ -197,10 +198,12 @@ struct Workspace{ObsScheme,S,TX,TW,R,TP,TZ}# ,Q, where Q = eltype(result) end y = XX[i].yy[end] end + y = x0_guess ll = ( logpdf(x0_prior, y) + path_log_likhd(ObsScheme(), XX, P, 1:m, fpt, skipFPT=true) + lobslikelihood(P[1], y) ) + TL = typeof(ll) XXᵒ, WWᵒ, Pᵒ = deepcopy(XX), deepcopy(WW), deepcopy(P) @@ -213,17 +216,16 @@ struct Workspace{ObsScheme,S,TX,TW,R,TP,TZ}# ,Q, where Q = eltype(result) _time = collect(Iterators.flatten(p.tt[1:skip:end-1] for p in P)) #check_if_recompute_ODEs(setup) - (workspace = new{ObsScheme,S,TX,TW,R,TP,TZ}(Wnr, XXᵒ, XX, WWᵒ, WW, Pᵒ, + new{ObsScheme,S,TX,TW,R,TP,TL,TZ}(Wnr, XXᵒ, XX, WWᵒ, WW, Pᵒ, P, fpt, skip, [], _time, - x0_prior, z), - ll = ll, θ = params(P[1].Target)) + x0_prior, ll, z) end function Workspace(ws::Workspace{ObsScheme,S,TX,TW,R̃,TP,TZ}, P::Vector{R}, Pᵒ::Vector{R}) where {ObsScheme,S,TX,TW,R̃,R,TP,TZ} - new{ObsScheme,S,TX,TW,R,TP,TZ}(ws.Wnr, ws.XXᵒ, ws.XX, ws.WWᵒ, ws.WW, Pᵒ, + new{ObsScheme,S,TX,TW,R,TP,TL,TZ}(ws.Wnr, ws.XXᵒ, ws.XX, ws.WWᵒ, ws.WW, Pᵒ, P, ws.fpt, ws.skip_for_save, ws.paths, - ws.time, ws.x0_prior, ws.z) + ws.time, ws.x0_prior, ws.ll, ws.z) end end @@ -274,5 +276,7 @@ function create_workspace(setup::MCMCSetup, schedule::MCMCSchedule, θ) end function create_workspace(setup::T) where {T <: ModelSetup} - Workspace(setup) + ws = Workspace(setup) + θ = params(ws.P[1].Target) + ws, ws.ll, θ end From 1c9f31c9ec8015dcdf6d6f1a71a39f4f5167bbda Mon Sep 17 00:00:00 2001 From: mschauer Date: Thu, 26 Mar 2020 13:43:12 +0100 Subject: [PATCH 29/29] Profile inner loop --- scripts/nextjournal/nextjournal_script.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/scripts/nextjournal/nextjournal_script.jl b/scripts/nextjournal/nextjournal_script.jl index 95046cc..fc7cfe5 100644 --- a/scripts/nextjournal/nextjournal_script.jl +++ b/scripts/nextjournal/nextjournal_script.jl @@ -92,4 +92,10 @@ schedule = MCMCSchedule(1*10^3, [[1,2,3,4,5]], Random.seed!(4) Profile.init() -@profile out = mcmc(mcmc_setup, schedule, model_setup); +setup_mcmc = mcmc_setup +setup = model_setup +ws, ll, θ = BridgeSDEInference.create_workspace(setup) +ws_mcmc = BridgeSDEInference.create_workspace(setup_mcmc, schedule, θ) +adpt = BridgeSDEInference.adaptation_object(setup, ws) +Profile.clear() +@profile out = BridgeSDEInference.mcmc_(setup_mcmc, schedule, setup, ws, ll, θ, ws_mcmc, adpt)