From fdc4031b1d162fa6ea7b33151062eb2fa875537e Mon Sep 17 00:00:00 2001 From: mschauer Date: Fri, 21 Sep 2018 23:40:16 +0200 Subject: [PATCH] Make project_partialbridge Julia 1.0 safe --- .../partialbridge_fitzhugh.jl | 21 ++-- project_partialbridge/partialbridge_nclar.jl | 119 +++++++++--------- 2 files changed, 75 insertions(+), 65 deletions(-) diff --git a/project_partialbridge/partialbridge_fitzhugh.jl b/project_partialbridge/partialbridge_fitzhugh.jl index 199d23f..bd4e0d0 100644 --- a/project_partialbridge/partialbridge_fitzhugh.jl +++ b/project_partialbridge/partialbridge_fitzhugh.jl @@ -1,5 +1,5 @@ -cd("/Users/Frank/.julia/dev/Bridge") -outdir="/Users/Frank/Sync/DOCUMENTS/onderzoek/code/diffbridges/out_fh/" +mkpath("output/out_fh") +outdir="output/out_fh/" using Bridge, StaticArrays, Distributions using Test, Statistics, Random, LinearAlgebra @@ -83,6 +83,7 @@ for k1 in (1:3) iterations = !(k1==3) ? 5*10^4 : 10*10^4 skip_it = 1000 subsamples = 0:skip_it:iterations + printiter = 100 if endpoint == "first" #v = -0.959 @@ -153,21 +154,25 @@ for k1 in (1:3) llo = llikelihood(Bridge.LeftRule(), Xo, Po) - if mod(iter,50)==0 - print(iter," ll $ll $llo, diff_ll: ",round(llo-ll,3))#, X[10], " ", Xo[10]) + if mod(iter, printiter) == 0 + print(iter," ll $ll $llo, diff_ll: ",round(llo-ll, digits=3))#, X[10], " ", Xo[10]) end - + accept = false if log(rand()) <= llo - ll X.yy .= Xo.yy W.yy .= Wo.yy ll = llo - #print("✓") + accept = true acc +=1 end - println() + if iter in subsamples push!(XX, copy(X)) end + if mod(iter, printiter) == 0 + accept && print("✓") + println() + end end @info "Done."*"\x7"^6 @@ -183,7 +188,7 @@ for k1 in (1:3) writedlm(f,iterates,',') close(f) - ave_acc_perc = 100*round(acc/iterations,2) + ave_acc_perc = 100*round(acc/iterations, digits=2) # write info to txt file fn = outdir*"info-"*endpoint*"-"*aux_choice*".txt" diff --git a/project_partialbridge/partialbridge_nclar.jl b/project_partialbridge/partialbridge_nclar.jl index 80b986d..6c36aa0 100644 --- a/project_partialbridge/partialbridge_nclar.jl +++ b/project_partialbridge/partialbridge_nclar.jl @@ -1,6 +1,7 @@ -# reminder, to type H*, do H\^+ -cd("/Users/Frank/.julia/dev/Bridge") -outdir="/Users/Frank/Sync/DOCUMENTS/onderzoek/code/diffbridges/out_nclar/" +# reminder, to type H⁺, do H\^+ + +mkpath("output/out_nclar") +outdir="output/out_nclar/" using Bridge, StaticArrays, Distributions using Test, Statistics, Random, LinearAlgebra @@ -16,7 +17,7 @@ tt = τ(T).(0.:dt:T) sk = 0 # skipped in evaluating loglikelihood -obs_scheme =["full","firstcomponent"][2] +obs_scheme =["full", "firstcomponent"][2] easy_conditioning = false # if true, then path to 1, else to 2 νHparam = false @@ -34,9 +35,9 @@ subsamples = 0:skip_it:iterations ρ = obs_scheme=="full" ? 0.85 : 0.95 if obs_scheme=="full" - L = SMatrix{3,3}(1.0I) -# Σ = SMatrix{3,3}(0.0I) - v = easy_conditioning ? ℝ{3}(1/32,1/4,1) : ℝ{3}(5/128,3/8,2) + L = SMatrix{3, 3}(1.0I) +# Σ = SMatrix{3, 3}(0.0I) + v = easy_conditioning ? ℝ{3}(1/32, 1/4, 1) : ℝ{3}(5/128, 3/8, 2) end if obs_scheme=="firstcomponent" L = @SMatrix [1. 0. 0.] @@ -45,7 +46,7 @@ if obs_scheme=="firstcomponent" end m, d = size(L) -Σ = SMatrix{m,m}(Σdiagel*I) +Σ = SMatrix{m, m}(Σdiagel*I) # specify target process struct NclarDiffusion <: ContinuousTimeProcess{ℝ{3}} @@ -54,8 +55,9 @@ struct NclarDiffusion <: ContinuousTimeProcess{ℝ{3}} σ::Float64 end -Bridge.b(t, x, P::NclarDiffusion) = ℝ{3}(x[2],x[3],-P.α * sin(P.ω * x[3])) -Bridge.σ(t, x, P::NclarDiffusion) = ℝ{3}(0.0, 0.0, P.σ) +Bridge.b(t, x, P::NclarDiffusion) = ℝ{3}(x[2], x[3], -P.α * sin(P.ω * x[3])) +Bridge.σ(t, P::NclarDiffusion) = ℝ{3}(0.0, 0.0, P.σ) +Bridge.σ(t, x, P::NclarDiffusion) = Bridge.σ(t, P) Bridge.constdiff(::NclarDiffusion) = true P = NclarDiffusion(2*3.0, 2pi, 1.0) @@ -74,16 +76,17 @@ end Random.seed!(42) Bridge.B(t, P::NclarDiffusionAux) = @SMatrix [0.0 1.0 0.0 ; 0.0 0.0 1.0 ; 0.0 0.0 0.0] -Bridge.β(t, P::NclarDiffusionAux) = ℝ{3}(0.0,0.0,0) -Bridge.σ(t, x, P::NclarDiffusionAux) = ℝ{3}(0.0,0.0, P.σ) +Bridge.β(t, P::NclarDiffusionAux) = ℝ{3}(0.0, 0.0, 0) +Bridge.σ(t, P::NclarDiffusionAux) = ℝ{3}(0.0, 0.0, P.σ) +Bridge.σ(t, x, P::NclarDiffusionAux) = Bridge.σ(t, P) Bridge.constdiff(::NclarDiffusionAux) = true -Bridge.b(t, x, P::NclarDiffusionAux) = Bridge.B(t,P) * x + Bridge.β(t,P) -Bridge.a(t, P::NclarDiffusionAux) = Bridge.σ(t,0,P) * Bridge.σ(t, 0, P)' +Bridge.b(t, x, P::NclarDiffusionAux) = Bridge.B(t, P) * x + Bridge.β(t, P) +Bridge.a(t, P::NclarDiffusionAux) = Bridge.σ(t, 0, P) * Bridge.σ(t, 0, P)' Pt = NclarDiffusionAux(P.α, P.ω, P.σ) # Solve Backward Recursion -Po = νHparam ? Bridge.PartialBridgeνH(tt, P, Pt, L, ℝ{m}(v),ϵ, Σ) : Bridge.PartialBridge(tt, P, Pt, L, ℝ{m}(v), Σ) +Po = νHparam ? Bridge.PartialBridgeνH(tt, P, Pt, L, ℝ{m}(v), ϵ, Σ) : Bridge.PartialBridge(tt, P, Pt, L, ℝ{m}(v), Σ) ####################### MH algorithm ################### W = sample(tt, Wiener()) @@ -92,7 +95,6 @@ Xo = copy(X) bridge!(Xo, x0, W, Po) bridge!(X, x0, W, Po) -ll = llikelihood(Bridge.LeftRule(), X, Po,skip=sk) # further initialisation Wo = copy(W) @@ -102,28 +104,31 @@ if 0 in subsamples push!(XX, copy(X)) end -acc = 0 - -for iter in 1:iterations - # Proposal - sample!(W2, Wiener()) - #ρ = rand(Uniform(0.95,1.0)) - Wo.yy .= ρ*W.yy + sqrt(1-ρ^2)*W2.yy - bridge!(Xo, x0, Wo, Po) - llo = llikelihood(Bridge.LeftRule(), Xo, Po,skip=sk) - print("ll $ll $llo, diff_ll: ",round(llo-ll,3)) - - if log(rand()) <= llo - ll - X.yy .= Xo.yy - W.yy .= Wo.yy - ll = llo - print("✓") - acc +=1 - end - println() - if iter in subsamples - push!(XX, copy(X)) +let + ll = llikelihood(Bridge.LeftRule(), X, Po, skip=sk) + acc = 0 + for iter in 1:iterations + # Proposal + sample!(W2, Wiener()) + #ρ = rand(Uniform(0.95, 1.0)) + Wo.yy .= ρ*W.yy + sqrt(1-ρ^2)*W2.yy + bridge!(Xo, x0, Wo, Po) + + llo = llikelihood(Bridge.LeftRule(), Xo, Po, skip=sk) + print("ll $ll $llo, diff_ll: ", round(llo - ll, digits = 3)) + + if log(rand()) <= llo - ll + X.yy .= Xo.yy + W.yy .= Wo.yy + ll = llo + print("✓") + acc += 1 + end + println() + if iter in subsamples + push!(XX, copy(X)) + end end end @@ -132,34 +137,34 @@ end # write mcmc iterates to csv file fn = outdir*"iterates-"*obs_scheme*".csv" -f = open(fn,"w") -head = "iteration, time, component, value \n" +f = open(fn, "w") +head = "iteration, time, component, value\n" write(f, head) -iterates = [Any[s, tt[j], d, XX[i].yy[j][d]] for d in 1:3, j in 1:length(X), (i,s) in enumerate(subsamples) ][:] -writecsv(f,iterates) +iterates = [Any[s, tt[j], d, XX[i].yy[j][d]] for d in 1:3, j in 1:length(X), (i, s) in enumerate(subsamples) ][:] +writedlm(f, iterates, ',') close(f) -ave_acc_perc = 100*round(acc/iterations,2) +ave_acc_perc = 100*round(acc/iterations, digits=2) # write info to txt file fn = outdir*"info-"*obs_scheme*".txt" -f = open(fn,"w") -write(f, "Choice of observation schemes: ",obs_scheme,"\n") -write(f, "Easy conditioning (means going up to 1 for the rough component instead of 2): ",string(easy_conditioning),"\n") -write(f, "Number of iterations: ",string(iterations),"\n") -write(f, "Skip every ",string(skip_it)," iterations, when saving to csv","\n\n") -write(f, "Starting point: ",string(x0),"\n") -write(f, "End time T: ", string(T),"\n") -write(f, "Endpoint v: ",string(v),"\n") -write(f, "Noise Sigma: ",string(Σ),"\n") -write(f, "L: ",string(L),"\n\n") -write(f,"Mesh width: ",string(dt),"\n") -write(f, "rho (Crank-Nicholsen parameter: ",string(ρ),"\n") -write(f, "Average acceptance percentage: ",string(ave_acc_perc),"\n\n") -write(f, "Backward type parametrisation in terms of nu and H? ",string(νHparam),"\n") +f = open(fn, "w") +write(f, "Choice of observation schemes: ", obs_scheme, "\n") +write(f, "Easy conditioning (means going up to 1 for the rough component instead of 2): ", string(easy_conditioning), "\n") +write(f, "Number of iterations: ", string(iterations), "\n") +write(f, "Skip every ", string(skip_it), " iterations, when saving to csv", "\n\n") +write(f, "Starting point: ", string(x0), "\n") +write(f, "End time T: ", string(T), "\n") +write(f, "Endpoint v: ", string(v), "\n") +write(f, "Noise Sigma: ", string(Σ), "\n") +write(f, "L: ", string(L), "\n\n") +write(f, "Mesh width: ", string(dt), "\n") +write(f, "rho (Crank-Nicholsen parameter: ", string(ρ), "\n") +write(f, "Average acceptance percentage: ", string(ave_acc_perc), "\n\n") +write(f, "Backward type parametrisation in terms of nu and H? ", string(νHparam), "\n") close(f) -println("Average acceptance percentage: ",ave_acc_perc,"\n") +println("Average acceptance percentage: ", ave_acc_perc, "\n") println(obs_scheme) println("Parametrisation of nu and H? ", νHparam)