From 700e5448c69658577e48a63e0c57ff038b10bee1 Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 12:03:46 +0100 Subject: [PATCH 01/20] Move lm_main --- scripts/lm_main.jl | 81 ++++++++++----------- src/lm_main.jl | 176 --------------------------------------------- 2 files changed, 39 insertions(+), 218 deletions(-) delete mode 100755 src/lm_main.jl diff --git a/scripts/lm_main.jl b/scripts/lm_main.jl index 50bccce..fc4e416 100755 --- a/scripts/lm_main.jl +++ b/scripts/lm_main.jl @@ -1,43 +1,36 @@ -#] add StaticArrays Distributions DelimitedFiles DataFrames CSV RCall SparseArrays LowRankApprox Trajectories -#] add ForwardDiff DiffResults TimerOutputs Plots RecursiveArrayTools NPZ +# using Bridge, StaticArrays, Distributions +# using Bridge:logpdfnormal +# using Test, Statistics, Random, LinearAlgebra +# using Bridge.Models +# using DelimitedFiles, DataFrames, RCall +# using CSV +# # install.packages(ggforce) +# using Base.Iterators, SparseArrays, LowRankApprox, Trajectories +# using ForwardDiff +# using Plots, PyPlot #using Makie +# using RecursiveArrayTools +# using NPZ # for reading python datafiles +# using Random using BridgeLandmarks -const LM = BridgeLandmarks -using Bridge, StaticArrays, Distributions -using Bridge: logpdfnormal -using Test, Statistics, Random, LinearAlgebra -using Bridge.Models -using DelimitedFiles, DataFrames, RCall -using CSV -# install.packages(ggforce) -using Base.Iterators, SparseArrays, LowRankApprox, Trajectories -using ForwardDiff +workdir = @__DIR__ +println(workdir) +cd(workdir) + +using RCall using Plots -using RecursiveArrayTools -using NPZ # for reading python datafiles using Random - -#workdir = @__DIR__ -#println(workdir) -#cd(workdir) +using Distributions +using NPZ # for reading python datafiles +using DataFrames +using DelimitedFiles +using CSV pyplot() - -const d = 2 -const TEST = false - -#include("nstate.jl") -#include("state.jl") -#include("models.jl") -#include("patches.jl") -include("plotlandmarks.jl") # keep, but presently unused as all is transferred to plotting in R include("generatedata.jl") -#include("lmguiding_mv.jl") -#include("update_initialstate.jl") -#include("lmguid.jl") # replacing lmguiding_mv and update_initialstate - +include("plotting.jl") Random.seed!(3) @@ -45,7 +38,7 @@ Random.seed!(3) ################################# start settings ################################# n = 6 # nr of landmarks models = [:ms, :ahs] -model = models[2] +model = models[1] println("model: ",model) ITER = 5 @@ -66,7 +59,7 @@ prior_γ = Exponential(5.0) datasets =["forwardsimulated", "shifted","shiftedextreme", "bear","heart","peach", "generatedstefan","forwardsimulated_multiple", "cardiac"] -dataset = datasets[3] +dataset = datasets[7] println("dataset: ",dataset) fixinitmomentato0 = true#false @@ -84,13 +77,15 @@ fixinitmomentato0 = true#false σ_c = 0.2 # update c to cᵒ as cᵒ = c * exp(σ_c * rnorm()) σ_γ = 0.2 # update γ to γᵒ as γᵒ = γ * exp(σ_γ * rnorm()) +η(n) = min(0.1, 10/sqrt(n)) # adaptation rate for adjusting tuning pars + #------------------------------------------------------------------ σobs = 0.01 # noise on observations #------------------------------------------------------------------ # set time grids dt = 0.01 -T = 1.0; t = 0.0:dt:T; tt_ = LM.tc(t,T) +T = 1.0; t = 0.0:dt:T; tt_ = BridgeLandmarks.tc(t,T) outdir = "./figs/" if false # to use later on, when plotting is transferred from R to Julia @@ -104,6 +99,10 @@ a = 2.0 # Hamiltonian kernel parameter (the larger, the stronger landmarks b c = 0.1 # multiplicative factor in kernel γ = 1.0 # Noise level +# a = 0.2 +# c = 0.1 +# γ = 0.002 + if model == :ms λ = 0.0; # Mean reversion par in MS-model = not the lambda of noise fields =# nfs = 0 # needs to have value for plotting purposes @@ -111,7 +110,7 @@ if model == :ms else db = 2.0 # domainbound nfstd = 1.0 # tau , width of noisefields - nfs = LM.construct_nfs(db, nfstd, γ) # 3rd argument gives average noise of positions (with superposition) + nfs = construct_nfs(db, nfstd, γ) # 3rd argument gives average noise of positions (with superposition) Ptrue = Landmarks(a, c, n, db, nfstd, nfs) end @@ -123,12 +122,12 @@ x0, xobs0, xobsT, Xf, Ptrue, pb, obs_atzero = generatedata(dataset,Ptrue,t,σob # step-size on mala steps for initial state if obs_atzero - δ = [0.0, 0.1] # in this case first comp is not used + δ = [0.0, 0.01] # in this case first comp is not used else δ = [0.0005, 0.005] end -(ainit, cinit, γinit) = (0.1, 0.1, 0.1) +(ainit, cinit, γinit) = (1.0, 1.0, 1.0) #(0.1, 0.1, 0.1) # (a,c,γ) if model == :ms P = MarslandShardlow(ainit, cinit, γinit, Ptrue.λ, Ptrue.n) elseif model == :ahs @@ -142,7 +141,6 @@ mT = zeros(PointF,P.n) # vector of momenta at time T used for constructing gui start = time() # to compute elapsed time - if obs_atzero # only update mom xobsT = [xobsT] # should be a vector xinit = State(xobs0, zeros(PointF,P.n)) @@ -163,17 +161,16 @@ elseif !obs_atzero & fixinitmomentato0 # multiple shapes only update pos, so i initstate_updatetypes = [:rmmala_pos] end -LM.plotlandmarkpositions[] = plotlandmarkpositions -anim, Xsave, parsave, objvals, perc_acc_pcn, accinfo = lm_mcmc(tt_, (xobs0,xobsT), σobs, mT, P, +anim, Xsave, parsave, objvals, acc_pcn, accinfo = lm_mcmc(tt_, (xobs0,xobsT), σobs, mT, P, sampler, obs_atzero, fixinitmomentato0, xinit, ITER, subsamples, - (ρ, δ, prior_a, prior_c, prior_γ, σ_a, σ_c, σ_γ), initstate_updatetypes, + (ρ, δ, prior_a, prior_c, prior_γ, σ_a, σ_c, σ_γ, η), initstate_updatetypes, adaptskip, outdir, pb; updatepars = true, makefig=true, showmomenta=false) elapsed = time() - start println("Elapsed time: ",round(elapsed/60;digits=2), " minutes") -println("Acceptance percentage pCN step: ",perc_acc_pcn) +println("Acceptance percentage pCN step: ",acc_pcn) include("./postprocessing.jl") diff --git a/src/lm_main.jl b/src/lm_main.jl deleted file mode 100755 index fc4e416..0000000 --- a/src/lm_main.jl +++ /dev/null @@ -1,176 +0,0 @@ -# using Bridge, StaticArrays, Distributions -# using Bridge:logpdfnormal -# using Test, Statistics, Random, LinearAlgebra -# using Bridge.Models -# using DelimitedFiles, DataFrames, RCall -# using CSV -# # install.packages(ggforce) -# using Base.Iterators, SparseArrays, LowRankApprox, Trajectories -# using ForwardDiff -# using Plots, PyPlot #using Makie -# using RecursiveArrayTools -# using NPZ # for reading python datafiles -# using Random - -using BridgeLandmarks - -workdir = @__DIR__ -println(workdir) -cd(workdir) - -using RCall -using Plots -using Random -using Distributions -using NPZ # for reading python datafiles -using DataFrames -using DelimitedFiles -using CSV - -pyplot() - -include("generatedata.jl") -include("plotting.jl") - - -Random.seed!(3) - -################################# start settings ################################# -n = 6 # nr of landmarks -models = [:ms, :ahs] -model = models[1] -println("model: ",model) - -ITER = 5 -subsamples = 0:1:ITER - -showplotσq = false # only for ahs model - -sampler =[:sgd, :mcmc][2] -println("sampler: ",sampler) - -#------------------------------------------------------------------ -# prior on θ = (a, c, γ) -prior_a = Exponential(5.0) -prior_c = Exponential(5.0) -prior_γ = Exponential(5.0) -#------------------------------------------------------------------ - -datasets =["forwardsimulated", "shifted","shiftedextreme", - "bear","heart","peach", - "generatedstefan","forwardsimulated_multiple", "cardiac"] -dataset = datasets[7] -println("dataset: ",dataset) - -fixinitmomentato0 = true#false -#------------------------------------------------------------------ -# for sgd (FIX LATER) -ϵ = 0.01 # sgd step size -ϵstep(i) = 1/(1+i)^(0.7) - -#------------------------------------------------------------------ -### MCMC tuning pars -ρ = 0.8 # pcN-step - -# proposal for θ = (a, c, γ) -σ_a = 0.2 # update a to aᵒ as aᵒ = a * exp(σ_a * rnorm()) -σ_c = 0.2 # update c to cᵒ as cᵒ = c * exp(σ_c * rnorm()) -σ_γ = 0.2 # update γ to γᵒ as γᵒ = γ * exp(σ_γ * rnorm()) - -η(n) = min(0.1, 10/sqrt(n)) # adaptation rate for adjusting tuning pars - -#------------------------------------------------------------------ -σobs = 0.01 # noise on observations - -#------------------------------------------------------------------ -# set time grids -dt = 0.01 -T = 1.0; t = 0.0:dt:T; tt_ = BridgeLandmarks.tc(t,T) - -outdir = "./figs/" -if false # to use later on, when plotting is transferred from R to Julia - outdir = "./figs/"* string(model) * "_" * string(sampler) *"_" * string(dataset) * "/" - if !isdir(outdir) mkdir(outdir) end -end -################################# end settings ################################# - -### Specify landmarks models -a = 2.0 # Hamiltonian kernel parameter (the larger, the stronger landmarks behave similarly) -c = 0.1 # multiplicative factor in kernel -γ = 1.0 # Noise level - -# a = 0.2 -# c = 0.1 -# γ = 0.002 - -if model == :ms - λ = 0.0; # Mean reversion par in MS-model = not the lambda of noise fields =# - nfs = 0 # needs to have value for plotting purposes - Ptrue = MarslandShardlow(a, c, γ, λ, n) -else - db = 2.0 # domainbound - nfstd = 1.0 # tau , width of noisefields - nfs = construct_nfs(db, nfstd, γ) # 3rd argument gives average noise of positions (with superposition) - Ptrue = Landmarks(a, c, n, db, nfstd, nfs) -end - -if (model == :ahs) & showplotσq - plotσq(db, nfs) -end - -x0, xobs0, xobsT, Xf, Ptrue, pb, obs_atzero = generatedata(dataset,Ptrue,t,σobs) - -# step-size on mala steps for initial state -if obs_atzero - δ = [0.0, 0.01] # in this case first comp is not used -else - δ = [0.0005, 0.005] -end - -(ainit, cinit, γinit) = (1.0, 1.0, 1.0) #(0.1, 0.1, 0.1) # (a,c,γ) -if model == :ms - P = MarslandShardlow(ainit, cinit, γinit, Ptrue.λ, Ptrue.n) -elseif model == :ahs - nfsinit = construct_nfs(Ptrue.db, Ptrue.nfstd, γinit) - P = Landmarks(ainit, cinit, Ptrue.n, Ptrue.db, Ptrue.nfstd, nfsinit) -end - - -mT = zeros(PointF,P.n) # vector of momenta at time T used for constructing guiding term -#mT = randn(PointF,P.n) - -start = time() # to compute elapsed time - -if obs_atzero # only update mom - xobsT = [xobsT] # should be a vector - xinit = State(xobs0, zeros(PointF,P.n)) - initstate_updatetypes = [:mala_mom] -elseif !obs_atzero & !fixinitmomentato0 # multiple shapes: update both pos and mom - θ = π/6 - rot = SMatrix{2,2}(cos(θ), sin(θ), -sin(θ), cos(θ)) - xinit = 1.3*State([rot * xobsT[1][i] for i in 1:P.n], zeros(PointF,P.n)) - #xinit = 1.3*State([rot * xobsTvec[1][i] for i in 1:n], randn(PointF,P.n)) - #xinit = State(xobsT[1], zeros(PointF, P.n)) - initstate_updatetypes = [:mala_mom, :rmmala_pos] -elseif !obs_atzero & fixinitmomentato0 # multiple shapes only update pos, so initialised momenta (equal toz zero) are maintained throughout the iterations - xinit = 1.2*State(xobsT[1], zeros(PointF,P.n)) - θ = π/6 - rot = SMatrix{2,2}(cos(θ), sin(θ), -sin(θ), cos(θ)) - xinit = 1.3*State([rot * xobsT[1][i] for i in 1:P.n], zeros(PointF,P.n)) - #xinit = State(xobsT[1], zeros(PointF, P.n)) - initstate_updatetypes = [:rmmala_pos] -end - - -anim, Xsave, parsave, objvals, acc_pcn, accinfo = lm_mcmc(tt_, (xobs0,xobsT), σobs, mT, P, - sampler, obs_atzero, fixinitmomentato0, - xinit, ITER, subsamples, - (ρ, δ, prior_a, prior_c, prior_γ, σ_a, σ_c, σ_γ, η), initstate_updatetypes, adaptskip, - outdir, pb; updatepars = true, makefig=true, showmomenta=false) - -elapsed = time() - start - -println("Elapsed time: ",round(elapsed/60;digits=2), " minutes") -println("Acceptance percentage pCN step: ",acc_pcn) - -include("./postprocessing.jl") From dad073981d858ee404a0d1949d6ade2b004b34d2 Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 12:04:15 +0100 Subject: [PATCH 02/20] Move generatedata.jl --- scripts/generatedata.jl | 2 +- src/generatedata.jl | 195 ---------------------------------------- 2 files changed, 1 insertion(+), 196 deletions(-) delete mode 100644 src/generatedata.jl diff --git a/scripts/generatedata.jl b/scripts/generatedata.jl index da6c1fe..d1f06bd 100644 --- a/scripts/generatedata.jl +++ b/scripts/generatedata.jl @@ -21,7 +21,7 @@ function generatedata(dataset,P,t,σobs) n = P.n if dataset=="forwardsimulated" q0 = [PointF(2.0cos(t), sin(t)) for t in (0:(2pi/n):2pi)[1:n]] #q0 = circshift(q0, (1,)) - p0 = [Point(1.0, -3.0) for i in 1:n]/n # #p0 = [randn(Point) for i in 1:n] + p0 = [5*Point(1.0, -3.0) for i in 1:n]/n # #p0 = [randn(Point) for i in 1:n] x0 = State(q0, p0) Wf, Xf = landmarksforward(t, x0, P) xobs0 = x0.q + σobs * randn(PointF,n) diff --git a/src/generatedata.jl b/src/generatedata.jl deleted file mode 100644 index d1f06bd..0000000 --- a/src/generatedata.jl +++ /dev/null @@ -1,195 +0,0 @@ -""" -Generate landmarks data, or read landmarks data. -- dataset: specifies type of data to be generated/loaded -- P: -- t: grid used for generating data -- σobs: standard deviation of observation noise - -Returns -- x0: initial state without noise -- xobs0, xobsT: (observed initial and final states with noise) -- Xf: forward simulated path -- P: only adjusted in case of 'real' data, like the bearskull data -- pb: plotting bounds used in animation - -Example - x0, xobs0, xobsT, Xf, P, pb = generatedata(dataset,P,t,σobs) - -""" -function generatedata(dataset,P,t,σobs) - - n = P.n - if dataset=="forwardsimulated" - q0 = [PointF(2.0cos(t), sin(t)) for t in (0:(2pi/n):2pi)[1:n]] #q0 = circshift(q0, (1,)) - p0 = [5*Point(1.0, -3.0) for i in 1:n]/n # #p0 = [randn(Point) for i in 1:n] - x0 = State(q0, p0) - Wf, Xf = landmarksforward(t, x0, P) - xobs0 = x0.q + σobs * randn(PointF,n) - xobsT = [Xf.yy[end].q[i] for i in 1:P.n ] + σobs * randn(PointF,n) - pb = Lmplotbounds(-2.0,2.0,-1.5,1.5) - obs_atzero = true - end - if dataset in ["shifted","shiftedextreme"] # first stretch, then rotate, then shift; finally add noise - q0 = [PointF(2.0cos(t), sin(t)) for t in (0:(2pi/n):2pi)[1:n]] #q0 = circshift(q0, (1,)) - p0 = [PointF(1.0, -3.0) for i in 1:n]/n # #p0 = [randn(Point) for i in 1:n] - x0 = State(q0, p0) - @time Wf, Xf = landmarksforward(t, x0, P) - xobs0 = x0.q + σobs * randn(PointF,n) - if dataset == "shifted" - θ, η = π/10, 0.2 - pb = Lmplotbounds(-3.0,3.0,-2.0,2.0) # verified - end - if dataset == "shiftedextreme" - θ, η = π/5, 0.4 - pb = Lmplotbounds(-3.0,3.0,-3.0,3.0) - end - rot = SMatrix{2,2}(cos(θ), sin(θ), -sin(θ), cos(θ)) - stretch = SMatrix{2,2}(1.0 + η, 0.0, 0.0, 1.0 - η) - shift = 0.0*PointF(0.1,-0.1) - xobsT = [rot * stretch * xobs0[i] + shift for i in 1:P.n ] + σobs * randn(PointF,n) - obs_atzero = true - end - if dataset=="bear" - bear0 = readdlm(joinpath(@__DIR__,"data-stefan", "bear1.csv"), ',') - bearT = readdlm(joinpath(@__DIR__,"data-stefan", "bear2.csv"), ',') - nb = size(bear0)[1] - avePoint = Point(414.0, 290.0) # average of (x,y)-coords for bear0 to center figure at origin - xobs0 = [Point(bear0[i,1], bear0[i,2]) - avePoint for i in 1:nb]/200. - xobsT = [Point(bearT[i,1], bearT[i,2]) - avePoint for i in 1:nb]/200. - # need to redefine P, because of n - if model == :ms - P = MarslandShardlow(P.a, P.c, P.γ, P.λ, nb) - else - P = Landmarks(P.a, P.c, nb, P.db, P.nfstd, P.nfs) - end - x0 = State(xobs0, rand(PointF,P.n)) - Wf, Xf = landmarksforward(t, x0, P) - pb = Lmplotbounds(-2.0,2.0,-1.5,1.5) - obs_atzero = true - end - if dataset=="heart" - q0 = [PointF(2.0cos(t), 2.0sin(t)) for t in (0:(2pi/n):2pi)[1:n]] #q0 = circshift(q0, (1,)) - p0 = [PointF(1.0, -3.0) for i in 1:n]/n # #p0 = [randn(Point) for i in 1:n] - x0 = State(q0, p0) - @time Wf, Xf = landmarksforward(t, x0, P) - xobs0 = x0.q + σobs * randn(PointF,n) - heart_xcoord(s) = 0.2*(13cos(s)-5cos(2s)-2cos(3s)-cos(4s)) - heart_ycoord(s) = 0.2*16(sin(s)^3) - qT = [PointF(heart_xcoord(t), heart_ycoord(t)) for t in (0:(2pi/n):2pi)[1:n]] #q0 = circshift(q0, (1,)) - xobsT = qT + σobs * randn(PointF,n) - pb = Lmplotbounds(-2.0,2.0,-1.5,1.5) - obs_atzero = true - end - if dataset=="peach" - q0 = [PointF(2.0cos(t), 2.0sin(t)) for t in (0:(2pi/n):2pi)[1:n]] #q0 = circshift(q0, (1,)) - p0 = [PointF(1.0, -3.0) for i in 1:n]/n # #p0 = [randn(Point) for i in 1:n] - x0 = State(q0, p0) - @time Wf, Xf = landmarksforward(t, x0, P) - xobs0 = x0.q + σobs * randn(PointF,n) - peach_xcoord(s) = (2.0 + sin(s)^3) * cos(s) - peach_ycoord(s) = (2.0 + sin(s)^3) * sin(s) - qT = [PointF(peach_xcoord(t), peach_ycoord(t)) for t in (0:(2pi/n):2pi)[1:n]] #q0 = circshift(q0, (1,)) - xobsT = qT + σobs * randn(PointF,n) - pb = Lmplotbounds(-2.0,2.0,-1.5,1.5) - obs_atzero = true - end - if dataset=="generatedstefan" - testshapes = npzread(joinpath(@__DIR__,"data-stefan", "match.npy.npz")) - xobs0vec = get(testshapes,"q0",0) - xobsTvec = get(testshapes,"v",0) - p0vec = get(testshapes,"p",0) - nb = div(length(xobs0vec),2) - - subs = 1:6:nb#1:5:nb - xobs0 = [PointF(xobs0vec[2i-1],xobs0vec[2i]) for i in subs] - xobsT = [PointF(xobsTvec[2i-1],xobsTvec[2i]) for i in subs] - p0 = [PointF(p0vec[2i-1],p0vec[2i]) for i in subs] - nb = length(subs) - - # need to redefine P, because of n - if model == :ms - P = MarslandShardlow(P.a,P.c, P.γ, P.λ, nb) - elseif model== :ahs - P = Landmarks(P.a,P.c, nb, P.db, P.nfstd, P.nfs) - end - x0 = State(xobs0, p0) - Wf, Xf = landmarksforward(t, x0, P) -# plotlandmarkpositions(Xf,P,xobs0,xobsT;db=4) - pb = Lmplotbounds(-2.0,2.0,-1.5,1.5) # verified - obs_atzero = true - end - - if dataset=="forwardsimulated_multiple" - nshapes = 3 - q0 = [PointF(2.0cos(t), sin(t)) for t in (0:(2pi/n):2pi)[1:n]] #q0 = circshift(q0, (1,)) - p0 = [Point(0.1, -0.4) for i in 1:n]/n # #p0 = [randn(Point) for i in 1:n] - #p0 = [Point(4.1, -0.4) for i in 1:n]/n # #p0 = [randn(Point) for i in 1:n] - x0 = State(q0, p0) - xobsT = Vector{PointF}[] - for k in 1:nshapes - #a = P.a * exp(0.15*randn()) - a = 5*P.a * exp(0.15*randn()) - c = P.c * exp(0.15*randn()) - γ = 0.6*getγ(P) * exp(0.01*randn()) - if model == :ms - P = MarslandShardlow(a, c, γ, P.λ, P.n) - elseif model == :ahs - nfs = construct_nfs(P.db, P.nfstd, γ) - P = Landmarks(a, c, P.n, P.db, P.nfstd, nfs) - end - Wf, Xf = landmarksforward(t, x0, P) - - push!(xobsT, [Xf.yy[end].q[i] for i in 1:P.n ] + σobs * randn(PointF,n)) - end - xobs0 = [] - pb = Lmplotbounds(-2.0,2.0,-1.5,1.5) - obs_atzero = false - end - if dataset=="cardiac" - # see https://arxiv.org/pdf/1805.06038.pdf - # 14 cardiac images, 66 landmarks per image - # The epicardial and endocardial contours where annotated with 33 landmarks along each outline resulting in 66 landmarks per image - # Each of the 14 rows is a shape with 66 landmarks, x coordinates in - # cardiac[:,:,0], y coordinates in cardiac[:,:,1], i.e. in python, all - # shapes are plotted with - # for i in range(N_samples): - # plt.plot(qs[i,:,0],qs[i,:,1]) - - cardiac = npzread(joinpath(@__DIR__,"data-stefan","cardiac.npy")) # heart data (left ventricles, the one we used in https://arxiv.org/abs/1705.10943 - cardiacx = cardiac[:,:,1] # x-coordinates of landmarks - cardiacy = cardiac[:,:,2] # y-coordinates of landmarks - nshapes = 14 - landmarksset = 1:5:66 - n = length(landmarksset) - xobsT = fill(zeros(PointF,66),nshapes) - for i in 1:nshapes # for each image - xobsT[i] = [PointF(cardiacx[i,j], cardiacy[i,j]) for j in landmarksset ] - end - x0 = State(zeros(PointF,n), zeros(PointF,n)) - pb = Lmplotbounds(-2.0,2.0,-1.5,1.5) # need to check - obs_atzero = false - xobs0 = [] - Xf = 0 - if model == :ms - P = MarslandShardlow(P.a,P.c, P.γ, P.λ, n) - elseif model== :ahs - P = Landmarks(P.a,P.c, n, P.db, P.nfstd, P.nfs) - end - - xobsTdf = DataFrame(x=extractcomp(xobsT[1],1), y =extractcomp(xobsT[1],2)) - @rput xobsTdf - R""" - library(tidyverse) - library(ggplot2) - xobsTdf %>% ggplot(aes(x=x,y=y)) + geom_point() + geom_path() - """ - pb = Lmplotbounds(-0.25,0.25,-0.25,0.25) - end - x0, xobs0, xobsT, Xf, P, pb, obs_atzero -end - -if false - cc = npzread(joinpath(@__DIR__,"data-stefan","cc.npy")) # corpus callosum data - - # there are 14 cardiac images, in https://arxiv.org/pdf/1705.10943.pdf 17 landmarks are chosen -end From 73a6e208eadd6f7c76ef1e8561fa72214003f24d Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 12:04:34 +0100 Subject: [PATCH 03/20] include update_initialstate --- src/BridgeLandmarks.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BridgeLandmarks.jl b/src/BridgeLandmarks.jl index 8e87b78..b31b368 100644 --- a/src/BridgeLandmarks.jl +++ b/src/BridgeLandmarks.jl @@ -39,6 +39,6 @@ include("patches.jl") #include("plotlandmarks.jl") # keep, but presently unused as all is transferred to plotting in R include("plotting.jl") include("lmguid.jl") # replacing lmguiding_mv and update_initialstate - +include("update_initialstate.jl") end # module From 139693c35449f70d2897f1d72ff184ab7ad8dc30 Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 12:04:46 +0100 Subject: [PATCH 04/20] Trash --- src/.Rhistory | 512 -------------------------------------------------- 1 file changed, 512 deletions(-) delete mode 100644 src/.Rhistory diff --git a/src/.Rhistory b/src/.Rhistory deleted file mode 100644 index 8958125..0000000 --- a/src/.Rhistory +++ /dev/null @@ -1,512 +0,0 @@ -mse.mle -mse.ebayes -mse.bayes -mse.posteriormean -dd <- data.frame(i=rep(1:n,4),value= -c(mle-theta0,ebayes-theta0,bayes-theta0,posteriormean-theta0), -type=rep(c('Diff Mle','Diff Emp. Bayes','Diff Bayes','Hier. Bayes'),each=n)) -ggplot(data=dd, aes(x=i,y=value,colour=type)) + geom_point()+theme_minimal()+facet_wrap(~type,nrow=1)+ -theme(legend.position="none")+ geom_hline(yintercept=0,size=1) -lambda -hist(lambda,breaks="FD") -x -1/mean(x) -alpha = 0.1*10 -beta = 0.1*10 -IT <- 10000 -lambda <- rep(0,IT) -theta <- matrix(0,IT,n) # each row contains an MCMC iteration -theta[1,] <- x # initialise at mle -for (it in 2:IT) -{ lambda[it] = rgamma(1,shape=2*n+alpha, rate= beta+sum(theta[it-1,])) -theta[it,] = x + rexp(1,rate=lambda[it]) -} -posteriormean <- as.tibble(theta) %>% filter(row_number() > IT/2) %>% summarise_all(mean) -posteriormean <- as.numeric(posteriormean[1,]) -mse.posteriormean <- sum((posteriormean-theta0)^2) -mse.mle -mse.ebayes -mse.bayes -mse.posteriormean -dd <- data.frame(i=rep(1:n,4),value= -c(mle-theta0,ebayes-theta0,bayes-theta0,posteriormean-theta0), -type=rep(c('Diff Mle','Diff Emp. Bayes','Diff Bayes','Hier. Bayes'),each=n)) -ggplot(data=dd, aes(x=i,y=value,colour=type)) + geom_point()+theme_minimal()+facet_wrap(~type,nrow=1)+ -theme(legend.position="none")+ geom_hline(yintercept=0,size=1) -n <- 10 -theta0 <- runif(n,0,1000) -# generate data -x<- rep(0,n) -for (i in 1:n) x[i] <- runif(1,0,theta0[i]) -# compute estiamtes -mle <- x -ebayes <- x+mean(x) -lambda_hyp <- 1 # hyperpar in prior -bayes <- x+1/lambda_hyp -mse.mle <- sum((mle-theta0)^2) -mse.ebayes <- sum((ebayes-theta0)^2) -mse.bayes <- sum((bayes-theta0)^2) -mse.mle -mse.ebayes -mse.bayes -# visualise results -dd <- data.frame(i=rep(1:n,3),value=c(mle-theta0,ebayes-theta0,bayes-theta0),type=rep(c('Diff Mle','Diff Emp. Bayes','Diff Bayes'),each=n)) -ggplot(data=dd, aes(x=i,y=value,colour=type)) + geom_point()+theme_minimal()+facet_wrap(~type)+ -theme(legend.position="bottom")+ geom_hline(yintercept=0,size=1.1) -# fully bayesian with Gibbs sampling -# specify hyperpars -alpha = 0.1*10 -beta = 0.1*10 -IT <- 10000 -lambda <- rep(0,IT) -theta <- matrix(0,IT,n) # each row contains an MCMC iteration -theta[1,] <- x # initialise at mle -for (it in 2:IT) -{ lambda[it] = rgamma(1,shape=2*n+alpha, rate= beta+sum(theta[it-1,])) -theta[it,] = x + rexp(1,rate=lambda[it]) -} -posteriormean <- as.tibble(theta) %>% filter(row_number() > IT/2) %>% summarise_all(mean) -posteriormean <- as.numeric(posteriormean[1,]) -mse.posteriormean <- sum((posteriormean-theta0)^2) -mse.mle -mse.ebayes -mse.bayes -mse.posteriormean -dd <- data.frame(i=rep(1:n,4),value= -c(mle-theta0,ebayes-theta0,bayes-theta0,posteriormean-theta0), -type=rep(c('Diff Mle','Diff Emp. Bayes','Diff Bayes','Hier. Bayes'),each=n)) -ggplot(data=dd, aes(x=i,y=value,colour=type)) + geom_point()+theme_minimal()+facet_wrap(~type,nrow=1)+ -theme(legend.position="none")+ geom_hline(yintercept=0,size=1) -1/mean(x) -hist(lambda) -hist(lambda) -setwd("~/.julia/dev/BridgeLandmarks/src") -setwd("~/.julia/dev/BridgeLandmarks/src") -library(tidyverse) -library(stringr) -library(gridExtra) -library(ggforce) -library(GGally) -theme_set(theme_light()) -######## read observations -obs0df <- read_delim("obs0.csv", ";", escape_double = FALSE, trim_ws = TRUE) -obsTdf <- read_delim("obsT.csv", ";", escape_double = FALSE, trim_ws = TRUE) %>% spread(key=pos,value=value) %>% -mutate(shape=as.factor(shape)) -v0 <- bind_rows(obs0df, obs0df) -vT <- bind_rows(obsTdf, obsTdf) -vT1 <- vT %>% dplyr::filter(shape=="1") -setwd("~/.julia/dev/BridgeLandmarks/src") -library(tidyverse) -library(stringr) -library(gridExtra) -library(ggforce) -library(GGally) -theme_set(theme_light()) -######## read observations -obs0df <- read_delim("obs0.csv", ";", escape_double = FALSE, trim_ws = TRUE) -obsTdf <- read_delim("obsT.csv", ";", escape_double = FALSE, trim_ws = TRUE) %>% spread(key=pos,value=value) %>% -mutate(shape=as.factor(shape)) -v0 <- bind_rows(obs0df, obs0df) -vT <- bind_rows(obsTdf, obsTdf) -vT1 <- vT %>% dplyr::filter(shape=="1") -####### read noisefields -nfsdf <- read_delim("noisefields.csv",";", escape_double = FALSE, trim_ws = TRUE) -######## read and analyse bridge iterates -d <- read_table2("iterates.csv") %>% -mutate(landmarkid=as.factor(landmarkid)) %>% -gather(key="iteratenr",value="value",contains("iter")) -d <- d %>% mutate(i = rep(1:(nrow(d)/4),each=4)) # add column to identify a particular combination (pos1,pos2,mom1,mom2) of one landmark at a specific time of a specific shape -d <- d %>% spread(key=pqtype, value=value) %>% dplyr::select(-i) %>% # can drop identifier column -mutate(iteratenr=str_replace(iteratenr,"iter",""),iterate=as.numeric(str_replace(iteratenr,"iter",""))) -d -# select all data for shape 1 -d1 <- d %>% dplyr::filter(shapes==1) -dsub <- d1 %>% dplyr::filter(iterate %in% c(0,5,10,25,600,700)) %>% -mutate(iteratenr = fct_relevel(iteratenr, c("0","5","10","25"))) -# v0 <- obsdf[1:n,] -dlabel0 <- obs0df; dlabel0$landmarkid <- unique(d$landmarkid) -# v0 <- rbind(v0,v0[1,]) -# -# vT = obsdf[seq(n+1,2*n),] -dlabelT <- obsTdf; dlabelT$landmarkid <- unique(d$landmarkid) -# vT <- rbind(vT,vT[1,]) -####### read parameter updates -parsdf <- read_delim("parameters.csv", ";", escape_double = FALSE, trim_ws = TRUE) -#------------------ figures -------------------------------------------------------------- -# plots shapes and noisefields -shapes <- ggplot() + -geom_path(data=v0, aes(x=pos1,y=pos2), colour='black')+ -geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange') + -geom_point(data=nfsdf, aes(x=locx, y=locy), color="Grey")+ -geom_circle(aes(x0 = locx, y0 = locy, r = nfstd), data = nfsdf,color="Grey",linetype="dashed")+ -theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + -geom_label(data=dlabel0, aes(x=pos1,y=pos2,label=landmarkid))+ -geom_label(data=dlabelT, aes(x=pos1,y=pos2,label=landmarkid),colour="orange") -pdf("shapes-noisefields.pdf",width=6,height=4) -show(shapes) -dev.off() -# plot paths of landmarks positions and bridges over various iterations -p4 <- dsub %>% ggplot(aes(x=pos1,y=pos2)) + -geom_path(aes(group=interaction(landmarkid,iteratenr),colour=time)) + facet_wrap(~iteratenr) + -geom_point(data=v0, aes(x=pos1,y=pos2), colour='black')+geom_point(data=vT, aes(x=pos1,y=pos2), colour='orange')+ -geom_path(data=v0, aes(x=pos1,y=pos2), colour='black')+geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange') + -theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + coord_fixed() -pdf("bridges-faceted.pdf",width=6,height=4) -show(p4) -dev.off() -# plot overlaid landmark bridges -p1 <- d1 %>% dplyr::filter(iterate %in% seq(0,max(d$iterate),by=5)) %>% ggplot() + -geom_path(aes(pos1,y=pos2,group=interaction(landmarkid,iteratenr),colour=iterate)) + -scale_colour_gradient(low="orange",high="darkblue")+ -geom_point(data=v0, aes(x=pos1,y=pos2), colour='black')+geom_point(data=vT, aes(x=pos1,y=pos2), colour='orange')+ -geom_path(data=v0, aes(x=pos1,y=pos2), colour='black',size=1.1)+geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange',size=1.1) + -theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + -geom_label(data=dlabel0, aes(x=pos1,y=pos2,label=landmarkid,hjust="outward",vjust="outward"))+ -geom_point(data=nfsdf, aes(x=locx, y=locy), color="Grey")+ -geom_circle(aes(x0 = locx, y0 = locy, r = nfstd), data = nfsdf,color="Grey",linetype="dashed")+ -coord_fixed() -pdf("bridges-overlaid.pdf",width=6,height=4) -show(p1) -dev.off() -# plot parameter updates -ppar1 <- parsdf %>% mutate(cdivgamma2=c/gamma^2) %>% gather(key=par, value=value, a, c, gamma,cdivgamma2) -ppar1$par <- factor(ppar1$par, levels=c('a', 'c', 'gamma','cdivgamma2'), labels=c("a","c",expression(gamma),expression(c/gamma^2))) -tracepars <- ppar1 %>% ggplot(aes(x=iterate, y=value)) + geom_path() + facet_wrap(~par, scales="free_y",labeller = label_parsed) + -xlab("iterate") + ylab("") + theme(strip.text.x = element_text(size = 12)) -pdf("trace-pars.pdf",width=6,height=4) -show(tracepars) -dev.off() -# pairwise scatter plots for parameter updates -ppar2 <- parsdf %>% ggplot(aes(x=a,y=c,colour=iterate)) + geom_point() + theme(legend.position = 'none') +scale_colour_gradient(low="orange",high="darkblue") -ppar3 <- parsdf %>% ggplot(aes(x=a,y=gamma,colour=iterate)) + geom_point() + theme(legend.position = 'none') +scale_colour_gradient(low="orange",high="darkblue") -ppar4 <- parsdf %>% ggplot(aes(x=c,y=gamma,colour=iterate)) + geom_point()+ theme(legend.position = 'none') +scale_colour_gradient(low="orange",high="darkblue") -pdf("scatter-pars.pdf",width=6,height=2) -grid.arrange(ppar2,ppar3,ppar4,ncol=3) -dev.off() -# plot paths of landmarks momenta -pmom <- d %>% dplyr::filter(time==0) %>% ggplot(aes(x=mom1,y=mom2,colour=iterate)) + geom_point() + -# geom_path(aes(group=interaction(landmarkid,iteratenr),colour=iterate)) + -facet_wrap(~landmarkid) +scale_colour_gradient(low="orange",high="darkblue")+theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + -geom_hline(yintercept=0, linetype="dashed")+geom_vline(xintercept=0, linetype="dashed") -pdf("momenta-faceted.pdf",width=6,height=4) -show(pmom) -dev.off() -# plot initial positions for all shapes (only interesting in case initial shape is unobserved) -dtime0 <- d %>% dplyr::filter(time==0)# ,iterate>200) -# add factor for determining which phase of sampling -dtime0 <- dtime0 %>% mutate(phase = -case_when(iterate < quantile(iterate,1/3) ~ "initial", -between(iterate,quantile(iterate,1/3),quantile(iterate,2/3)) ~ "middle", -iterate >= quantile(iterate,1/3) ~ "end") ) %>% # reorder factor levels -mutate(phase = fct_relevel(phase, "initial", "middle")) -dtimeT <- d %>% dplyr::filter(time==1) -initshapes0 <- ggplot() + -geom_point(data=vT, aes(x=pos1,y=pos2), colour='grey')+ -geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='grey', linetype="dashed",size=0.4) + -geom_point(data=v0, aes(x=pos1,y=pos2), colour='red')+ -geom_path(data=v0, aes(x=pos1,y=pos2), colour='red',size=0.8,alpha=0.8) + -geom_path(data=bind_rows(dtime0,dtime0),aes(x=pos1,y=pos2,colour=iterate)) + -facet_wrap(~phase)+ coord_fixed() -initshapes0 -pdf("initial-shapes.pdf",width=6,height=2) -initshapes0 -dev.off() -setwd("~/.julia/dev/BridgeLandmarks/src") -library(tidyverse) -library(stringr) -library(gridExtra) -library(ggforce) -library(GGally) -theme_set(theme_light()) -######## read observations -obs0df <- read_delim("obs0.csv", ";", escape_double = FALSE, trim_ws = TRUE) -obsTdf <- read_delim("obsT.csv", ";", escape_double = FALSE, trim_ws = TRUE) %>% spread(key=pos,value=value) %>% -mutate(shape=as.factor(shape)) -v0 <- bind_rows(obs0df, obs0df) -vT <- bind_rows(obsTdf, obsTdf) -vT1 <- vT %>% dplyr::filter(shape=="1") -####### read noisefields -nfsdf <- read_delim("noisefields.csv",";", escape_double = FALSE, trim_ws = TRUE) -######## read and analyse bridge iterates -d <- read_table2("iterates.csv") %>% -mutate(landmarkid=as.factor(landmarkid)) %>% -gather(key="iteratenr",value="value",contains("iter")) -d <- d %>% mutate(i = rep(1:(nrow(d)/4),each=4)) # add column to identify a particular combination (pos1,pos2,mom1,mom2) of one landmark at a specific time of a specific shape -d <- d %>% spread(key=pqtype, value=value) %>% dplyr::select(-i) %>% # can drop identifier column -mutate(iteratenr=str_replace(iteratenr,"iter",""),iterate=as.numeric(str_replace(iteratenr,"iter",""))) -d -# select all data for shape 1 -d1 <- d %>% dplyr::filter(shapes==1) -dsub <- d1 %>% dplyr::filter(iterate %in% c(0,5,10,25,600,700)) %>% -mutate(iteratenr = fct_relevel(iteratenr, c("0","5","10","25"))) -# v0 <- obsdf[1:n,] -dlabel0 <- obs0df; dlabel0$landmarkid <- unique(d$landmarkid) -# v0 <- rbind(v0,v0[1,]) -# -# vT = obsdf[seq(n+1,2*n),] -dlabelT <- obsTdf; dlabelT$landmarkid <- unique(d$landmarkid) -# vT <- rbind(vT,vT[1,]) -####### read parameter updates -parsdf <- read_delim("parameters.csv", ";", escape_double = FALSE, trim_ws = TRUE) -#------------------ figures -------------------------------------------------------------- -# plots shapes and noisefields -shapes <- ggplot() + -geom_path(data=v0, aes(x=pos1,y=pos2), colour='black')+ -geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange') + -geom_point(data=nfsdf, aes(x=locx, y=locy), color="Grey")+ -geom_circle(aes(x0 = locx, y0 = locy, r = nfstd), data = nfsdf,color="Grey",linetype="dashed")+ -theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + -geom_label(data=dlabel0, aes(x=pos1,y=pos2,label=landmarkid))+ -geom_label(data=dlabelT, aes(x=pos1,y=pos2,label=landmarkid),colour="orange") -setwd("~/.julia/dev/BridgeLandmarks/src/figs") -library(tidyverse) -library(stringr) -library(gridExtra) -library(ggforce) -library(GGally) -theme_set(theme_light()) -######## read observations -obs0df <- read_delim("obs0.csv", ";", escape_double = FALSE, trim_ws = TRUE) -obsTdf <- read_delim("obsT.csv", ";", escape_double = FALSE, trim_ws = TRUE) %>% spread(key=pos,value=value) %>% -mutate(shape=as.factor(shape)) -v0 <- bind_rows(obs0df, obs0df) -vT <- bind_rows(obsTdf, obsTdf) -vT1 <- vT %>% dplyr::filter(shape=="1") -####### read noisefields -nfsdf <- read_delim("noisefields.csv",";", escape_double = FALSE, trim_ws = TRUE) -######## read and analyse bridge iterates -d <- read_table2("iterates.csv") %>% -mutate(landmarkid=as.factor(landmarkid)) %>% -gather(key="iteratenr",value="value",contains("iter")) -d <- d %>% mutate(i = rep(1:(nrow(d)/4),each=4)) # add column to identify a particular combination (pos1,pos2,mom1,mom2) of one landmark at a specific time of a specific shape -d <- d %>% spread(key=pqtype, value=value) %>% dplyr::select(-i) %>% # can drop identifier column -mutate(iteratenr=str_replace(iteratenr,"iter",""),iterate=as.numeric(str_replace(iteratenr,"iter",""))) -d -# select all data for shape 1 -d1 <- d %>% dplyr::filter(shapes==1) -dsub <- d1 %>% dplyr::filter(iterate %in% c(0,5,10,25,600,700)) %>% -mutate(iteratenr = fct_relevel(iteratenr, c("0","5","10","25"))) -# v0 <- obsdf[1:n,] -dlabel0 <- obs0df; dlabel0$landmarkid <- unique(d$landmarkid) -# v0 <- rbind(v0,v0[1,]) -# -# vT = obsdf[seq(n+1,2*n),] -dlabelT <- obsTdf; dlabelT$landmarkid <- unique(d$landmarkid) -# vT <- rbind(vT,vT[1,]) -####### read parameter updates -parsdf <- read_delim("parameters.csv", ";", escape_double = FALSE, trim_ws = TRUE) -#------------------ figures -------------------------------------------------------------- -# plots shapes and noisefields -shapes <- ggplot() + -geom_path(data=v0, aes(x=pos1,y=pos2), colour='black')+ -geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange') + -geom_point(data=nfsdf, aes(x=locx, y=locy), color="Grey")+ -geom_circle(aes(x0 = locx, y0 = locy, r = nfstd), data = nfsdf,color="Grey",linetype="dashed")+ -theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + -geom_label(data=dlabel0, aes(x=pos1,y=pos2,label=landmarkid))+ -geom_label(data=dlabelT, aes(x=pos1,y=pos2,label=landmarkid),colour="orange") -pdf("shapes-noisefields.pdf",width=6,height=4) -show(shapes) -dev.off() -# plot paths of landmarks positions and bridges over various iterations -p4 <- dsub %>% ggplot(aes(x=pos1,y=pos2)) + -geom_path(aes(group=interaction(landmarkid,iteratenr),colour=time)) + facet_wrap(~iteratenr) + -geom_point(data=v0, aes(x=pos1,y=pos2), colour='black')+geom_point(data=vT, aes(x=pos1,y=pos2), colour='orange')+ -geom_path(data=v0, aes(x=pos1,y=pos2), colour='black')+geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange') + -theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + coord_fixed() -pdf("bridges-faceted.pdf",width=6,height=4) -show(p4) -dev.off() -# plot overlaid landmark bridges -p1 <- d1 %>% dplyr::filter(iterate %in% seq(0,max(d$iterate),by=5)) %>% ggplot() + -geom_path(aes(pos1,y=pos2,group=interaction(landmarkid,iteratenr),colour=iterate)) + -scale_colour_gradient(low="orange",high="darkblue")+ -geom_point(data=v0, aes(x=pos1,y=pos2), colour='black')+geom_point(data=vT, aes(x=pos1,y=pos2), colour='orange')+ -geom_path(data=v0, aes(x=pos1,y=pos2), colour='black',size=1.1)+geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange',size=1.1) + -theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + -geom_label(data=dlabel0, aes(x=pos1,y=pos2,label=landmarkid,hjust="outward",vjust="outward"))+ -geom_point(data=nfsdf, aes(x=locx, y=locy), color="Grey")+ -geom_circle(aes(x0 = locx, y0 = locy, r = nfstd), data = nfsdf,color="Grey",linetype="dashed")+ -coord_fixed() -pdf("bridges-overlaid.pdf",width=6,height=4) -show(p1) -dev.off() -# plot parameter updates -ppar1 <- parsdf %>% mutate(cdivgamma2=c/gamma^2) %>% gather(key=par, value=value, a, c, gamma,cdivgamma2) -ppar1$par <- factor(ppar1$par, levels=c('a', 'c', 'gamma','cdivgamma2'), labels=c("a","c",expression(gamma),expression(c/gamma^2))) -tracepars <- ppar1 %>% ggplot(aes(x=iterate, y=value)) + geom_path() + facet_wrap(~par, scales="free_y",labeller = label_parsed) + -xlab("iterate") + ylab("") + theme(strip.text.x = element_text(size = 12)) -pdf("trace-pars.pdf",width=6,height=4) -show(tracepars) -dev.off() -# pairwise scatter plots for parameter updates -ppar2 <- parsdf %>% ggplot(aes(x=a,y=c,colour=iterate)) + geom_point() + theme(legend.position = 'none') +scale_colour_gradient(low="orange",high="darkblue") -ppar3 <- parsdf %>% ggplot(aes(x=a,y=gamma,colour=iterate)) + geom_point() + theme(legend.position = 'none') +scale_colour_gradient(low="orange",high="darkblue") -ppar4 <- parsdf %>% ggplot(aes(x=c,y=gamma,colour=iterate)) + geom_point()+ theme(legend.position = 'none') +scale_colour_gradient(low="orange",high="darkblue") -pdf("scatter-pars.pdf",width=6,height=2) -grid.arrange(ppar2,ppar3,ppar4,ncol=3) -dev.off() -# plot paths of landmarks momenta -pmom <- d %>% dplyr::filter(time==0) %>% ggplot(aes(x=mom1,y=mom2,colour=iterate)) + geom_point() + -# geom_path(aes(group=interaction(landmarkid,iteratenr),colour=iterate)) + -facet_wrap(~landmarkid) +scale_colour_gradient(low="orange",high="darkblue")+theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + -geom_hline(yintercept=0, linetype="dashed")+geom_vline(xintercept=0, linetype="dashed") -pdf("momenta-faceted.pdf",width=6,height=4) -show(pmom) -dev.off() -# plot initial positions for all shapes (only interesting in case initial shape is unobserved) -dtime0 <- d %>% dplyr::filter(time==0)# ,iterate>200) -# add factor for determining which phase of sampling -dtime0 <- dtime0 %>% mutate(phase = -case_when(iterate < quantile(iterate,1/3) ~ "initial", -between(iterate,quantile(iterate,1/3),quantile(iterate,2/3)) ~ "middle", -iterate >= quantile(iterate,1/3) ~ "end") ) %>% # reorder factor levels -mutate(phase = fct_relevel(phase, "initial", "middle")) -dtimeT <- d %>% dplyr::filter(time==1) -initshapes0 <- ggplot() + -geom_point(data=vT, aes(x=pos1,y=pos2), colour='grey')+ -geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='grey', linetype="dashed",size=0.4) + -geom_point(data=v0, aes(x=pos1,y=pos2), colour='red')+ -geom_path(data=v0, aes(x=pos1,y=pos2), colour='red',size=0.8,alpha=0.8) + -geom_path(data=bind_rows(dtime0,dtime0),aes(x=pos1,y=pos2,colour=iterate)) + -facet_wrap(~phase)+ coord_fixed() -initshapes0 -pdf("initial-shapes.pdf",width=6,height=2) -initshapes0 -dev.off() -# plot paths of landmarks positions and bridges over various iterations -p4 <- dsub %>% ggplot(aes(x=pos1,y=pos2)) + -geom_path(aes(group=interaction(landmarkid,iteratenr),colour=time)) + facet_wrap(~iteratenr) + -geom_point(data=v0, aes(x=pos1,y=pos2), colour='black')+geom_point(data=vT, aes(x=pos1,y=pos2), colour='orange')+ -geom_path(data=v0, aes(x=pos1,y=pos2), colour='black')+geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange') + -theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + coord_fixed() -pdf("bridges-faceted.pdf",width=6,height=4) -show(p4) -dev.off() -# plot overlaid landmark bridges -p1 <- d1 %>% dplyr::filter(iterate %in% seq(0,max(d$iterate),by=5)) %>% ggplot() + -geom_path(aes(pos1,y=pos2,group=interaction(landmarkid,iteratenr),colour=iterate)) + -scale_colour_gradient(low="orange",high="darkblue")+ -geom_point(data=v0, aes(x=pos1,y=pos2), colour='black')+geom_point(data=vT, aes(x=pos1,y=pos2), colour='orange')+ -geom_path(data=v0, aes(x=pos1,y=pos2), colour='black',size=1.1)+geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange',size=1.1) + -theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + -geom_label(data=dlabel0, aes(x=pos1,y=pos2,label=landmarkid,hjust="outward",vjust="outward"))+ -geom_point(data=nfsdf, aes(x=locx, y=locy), color="Grey")+ -geom_circle(aes(x0 = locx, y0 = locy, r = nfstd), data = nfsdf,color="Grey",linetype="dashed")+ -coord_fixed() -pdf("bridges-overlaid.pdf",width=6,height=4) -show(p1) -dev.off() -# plot parameter updates -ppar1 <- parsdf %>% mutate(cdivgamma2=c/gamma^2) %>% gather(key=par, value=value, a, c, gamma,cdivgamma2) -ppar1$par <- factor(ppar1$par, levels=c('a', 'c', 'gamma','cdivgamma2'), labels=c("a","c",expression(gamma),expression(c/gamma^2))) -tracepars <- ppar1 %>% ggplot(aes(x=iterate, y=value)) + geom_path() + facet_wrap(~par, scales="free_y",labeller = label_parsed) + -xlab("iterate") + ylab("") + theme(strip.text.x = element_text(size = 12)) -pdf("trace-pars.pdf",width=6,height=4) -show(tracepars) -dev.off() -# pairwise scatter plots for parameter updates -ppar2 <- parsdf %>% ggplot(aes(x=a,y=c,colour=iterate)) + geom_point() + theme(legend.position = 'none') +scale_colour_gradient(low="orange",high="darkblue") -ppar3 <- parsdf %>% ggplot(aes(x=a,y=gamma,colour=iterate)) + geom_point() + theme(legend.position = 'none') +scale_colour_gradient(low="orange",high="darkblue") -ppar4 <- parsdf %>% ggplot(aes(x=c,y=gamma,colour=iterate)) + geom_point()+ theme(legend.position = 'none') +scale_colour_gradient(low="orange",high="darkblue") -pdf("scatter-pars.pdf",width=6,height=2) -grid.arrange(ppar2,ppar3,ppar4,ncol=3) -dev.off() -# plot paths of landmarks momenta -pmom <- d %>% dplyr::filter(time==0) %>% ggplot(aes(x=mom1,y=mom2,colour=iterate)) + geom_point() + -# geom_path(aes(group=interaction(landmarkid,iteratenr),colour=iterate)) + -facet_wrap(~landmarkid) +scale_colour_gradient(low="orange",high="darkblue")+theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + -geom_hline(yintercept=0, linetype="dashed")+geom_vline(xintercept=0, linetype="dashed") -pdf("momenta-faceted.pdf",width=6,height=4) -show(pmom) -dev.off() -# plot initial positions for all shapes (only interesting in case initial shape is unobserved) -dtime0 <- d %>% dplyr::filter(time==0)# ,iterate>200) -# add factor for determining which phase of sampling -dtime0 <- dtime0 %>% mutate(phase = -case_when(iterate < quantile(iterate,1/3) ~ "initial", -between(iterate,quantile(iterate,1/3),quantile(iterate,2/3)) ~ "middle", -iterate >= quantile(iterate,1/3) ~ "end") ) %>% # reorder factor levels -mutate(phase = fct_relevel(phase, "initial", "middle")) -dtimeT <- d %>% dplyr::filter(time==1) -initshapes0 <- ggplot() + -geom_point(data=vT, aes(x=pos1,y=pos2), colour='grey')+ -geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='grey', linetype="dashed",size=0.4) + -geom_point(data=v0, aes(x=pos1,y=pos2), colour='red')+ -geom_path(data=v0, aes(x=pos1,y=pos2), colour='red',size=0.8,alpha=0.8) + -geom_path(data=bind_rows(dtime0,dtime0),aes(x=pos1,y=pos2,colour=iterate)) + -facet_wrap(~phase)+ coord_fixed() -initshapes0 -pdf("initial-shapes.pdf",width=6,height=2) -initshapes0 -dev.off() -rstudioapi::getSourceEditorContext()$path -dirname(rstudioapi::getSourceEditorContext()$path) -setwd(dirname(rstudioapi::getSourceEditorContext()$path)) -n <- 100 -theta0 <- runif(n,0,1000) -x<- rep(0,n) -for (i in 1:n) x[i] <- runif(1,0,theta0[i]) -mle <- x -ebayes <- x+mean(x) -lambda_hyp <- 1 # hyperpar in prior -bayes <- x+1/lambda_hyp -mse.mle <- sum((mle-theta0)^2) -mse.ebayes <- sum((ebayes-theta0)^2) -mse.bayes <- sum((bayes-theta0)^2) -mse.mle -mse.ebayes -mse.bayes -dd <- data.frame(i=rep(1:n,3),value=c(mle-theta0,ebayes-theta0,bayes-theta0),type=rep(c('Diff Mle','Diff Emp. Bayes','Diff Bayes'),each=n)) -ggplot(data=dd, aes(x=i,y=value,colour=type)) + geom_point()+theme_minimal()+facet_wrap(~type)+ -theme(legend.position="bottom")+ geom_hline(yintercept=0,size=1.1) -# fully bayesian with Gibbs sampling -# specify hyperpars -alpha = 0.1*10 -beta = 0.1*10 -IT <- 10000 -lambda <- rep(0,IT) -theta <- matrix(0,IT,n) # each row contains an MCMC iteration -theta[1,] <- x # initialise at mle -for (it in 2:IT) -{ lambda[it] = rgamma(1,shape=2*n+alpha, rate= beta+sum(theta[it-1,])) -theta[it,] = x + rexp(1,rate=lambda[it]) -} -posteriormean <- as.tibble(theta) %>% filter(row_number() > IT/2) %>% summarise_all(mean) -posteriormean <- as.numeric(posteriormean[1,]) -mse.posteriormean <- sum((posteriormean-theta0)^2) -mse.mle -mse.ebayes -mse.bayes -mse.posteriormean -dd <- data.frame(i=rep(1:n,4),value= -c(mle-theta0,ebayes-theta0,bayes-theta0,posteriormean-theta0), -type=rep(c('Diff Mle','Diff Emp. Bayes','Diff Bayes','Hier. Bayes'),each=n)) -ggplot(data=dd, aes(x=i,y=value,colour=type)) + geom_point()+theme_minimal()+facet_wrap(~type,nrow=1)+ -theme(legend.position="none")+ geom_hline(yintercept=0,size=1) -hist(lambda,breaks="FD") -1/mean(x) -# specify hyperpars -alpha = 0.1 -beta = 2 -IT <- 10000 -lambda <- rep(0,IT) -theta <- matrix(0,IT,n) # each row contains an MCMC iteration -theta[1,] <- x # initialise at mle -for (it in 2:IT) -{ lambda[it] = rgamma(1,shape=2*n+alpha, rate= beta+sum(theta[it-1,])) -theta[it,] = x + rexp(1,rate=lambda[it]) -} -posteriormean <- as.tibble(theta) %>% filter(row_number() > IT/2) %>% summarise_all(mean) -posteriormean <- as.numeric(posteriormean[1,]) -mse.posteriormean <- sum((posteriormean-theta0)^2) -mse.mle -mse.ebayes -mse.bayes -mse.posteriormean -dd <- data.frame(i=rep(1:n,4),value= -c(mle-theta0,ebayes-theta0,bayes-theta0,posteriormean-theta0), -type=rep(c('Diff Mle','Diff Emp. Bayes','Diff Bayes','Hier. Bayes'),each=n)) -ggplot(data=dd, aes(x=i,y=value,colour=type)) + geom_point()+theme_minimal()+facet_wrap(~type,nrow=1)+ -theme(legend.position="none")+ geom_hline(yintercept=0,size=1) -hist(lambda,breaks="FD") From 1074707c710895ad0652bc822993bdfe31f6c21a Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 12:05:03 +0100 Subject: [PATCH 05/20] Update compat --- Project.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 35807b9..328be5c 100644 --- a/Project.toml +++ b/Project.toml @@ -22,7 +22,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] Bridge = "0.10" -DataFrames = "0.20" +DataFrames = "0.20" Distributions = "0.22" ForwardDiff = "0.10" GaussianDistributions = "0.3" @@ -30,8 +30,8 @@ JLD = "0.9" PDMats = "0.9" Plots = "0.28" RecursiveArrayTools = "2.0" -StaticArrays = "0.12" -TimerOutputs = "0.5" +StaticArrays = "0.11, 0.12" +TimerOutputs = "0.5" julia = "1" [extras] From fe3cb118aec35126c12edb61f14b64a75f40771f Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 12:18:59 +0100 Subject: [PATCH 06/20] Clean warnings --- src/BridgeLandmarks.jl | 2 +- src/models.jl | 2 +- src/nstate.jl | 5 +++++ src/state.jl | 6 ------ 4 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/BridgeLandmarks.jl b/src/BridgeLandmarks.jl index b31b368..f99d419 100644 --- a/src/BridgeLandmarks.jl +++ b/src/BridgeLandmarks.jl @@ -39,6 +39,6 @@ include("patches.jl") #include("plotlandmarks.jl") # keep, but presently unused as all is transferred to plotting in R include("plotting.jl") include("lmguid.jl") # replacing lmguiding_mv and update_initialstate -include("update_initialstate.jl") +include("updatematching.jl") end # module diff --git a/src/models.jl b/src/models.jl index d1b58ae..b496aac 100755 --- a/src/models.jl +++ b/src/models.jl @@ -1,5 +1,5 @@ #### Landmarks specification -import Bridge: _b, _b!, B!, σ!, b!, σ, b +import Bridge: _b, _b!, B!, σ!, b!, σ, b, auxiliary struct MarslandShardlow{T} <: ContinuousTimeProcess{State{PointF}} a::T # kernel std parameter diff --git a/src/nstate.jl b/src/nstate.jl index 826498e..3c02f0e 100644 --- a/src/nstate.jl +++ b/src/nstate.jl @@ -95,6 +95,11 @@ import Base: *, +, /, - import LinearAlgebra: dot # import Bridge: outer, inner + +function outer(x::NState, y::NState) + [outer(x[i],y[j]) for i in eachindex(x), j in eachindex(y)] +end + *(c::Number, x::NState) = NState(c*x.x) *(x::NState, c::Number) = NState(x.x*c) +(x::NState, y::NState) = NState(x.x + y.x) diff --git a/src/state.jl b/src/state.jl index affba1f..0d716c0 100644 --- a/src/state.jl +++ b/src/state.jl @@ -22,12 +22,6 @@ function deepmat2unc(A::Matrix) # d is the dimension of the square subblocks end -function outer(x::State, y::State) - [outer(x[i],y[j]) for i in eachindex(x), j in eachindex(y)] -end - -norm(x::State) = norm(vec(x)) - """ Good display of variable of type State From 5fb481d43e954bb6a13c028e7a4290d40c171b2b Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 12:19:54 +0100 Subject: [PATCH 07/20] mv makefigs.R --- {src => scripts}/makefigs.R | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename {src => scripts}/makefigs.R (100%) diff --git a/src/makefigs.R b/scripts/makefigs.R similarity index 100% rename from src/makefigs.R rename to scripts/makefigs.R From bb4c3f3fbcc208b7b244777885d8e031a01af482 Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 12:21:15 +0100 Subject: [PATCH 08/20] Mv data-stefan --- {src/data-stefan => data-stefan}/cardiac.npy | Bin {src/data-stefan => data-stefan}/cc.npy | Bin {src/data-stefan => data-stefan}/explanation.pdf | Bin {src/data-stefan => data-stefan}/match.npy.npz | Bin 4 files changed, 0 insertions(+), 0 deletions(-) rename {src/data-stefan => data-stefan}/cardiac.npy (100%) rename {src/data-stefan => data-stefan}/cc.npy (100%) rename {src/data-stefan => data-stefan}/explanation.pdf (100%) rename {src/data-stefan => data-stefan}/match.npy.npz (100%) diff --git a/src/data-stefan/cardiac.npy b/data-stefan/cardiac.npy similarity index 100% rename from src/data-stefan/cardiac.npy rename to data-stefan/cardiac.npy diff --git a/src/data-stefan/cc.npy b/data-stefan/cc.npy similarity index 100% rename from src/data-stefan/cc.npy rename to data-stefan/cc.npy diff --git a/src/data-stefan/explanation.pdf b/data-stefan/explanation.pdf similarity index 100% rename from src/data-stefan/explanation.pdf rename to data-stefan/explanation.pdf diff --git a/src/data-stefan/match.npy.npz b/data-stefan/match.npy.npz similarity index 100% rename from src/data-stefan/match.npy.npz rename to data-stefan/match.npy.npz From 6f71f4159f9918e2bbaf3ba3e144ae0b8995fcd4 Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 12:24:00 +0100 Subject: [PATCH 09/20] Move experiements --- .../exp-cardiac/an_cardiac.jl | 0 .../exp-cardiac/cardiac.npy | Bin .../exp-cardiac/data_cardiac.jld | Bin .../exp-cardiac/gen_cardiac.jl | 0 .../exp-cardiac/makefigs_cardiac.R | 0 .../exp-dgs/an_exp-dgs.jl | 0 .../exp-dgs/data_exp-dgs.jld | Bin .../exp-dgs/gen_exp-dgs.jl | 0 .../exp-dgs/makefigs_exp-dgs.R | 0 .../exp-dgs/match.npy.npz | Bin .../exp0/animate_introexample.mp4 | Bin {src/experiments => experiments}/exp0/exp0.jl | 0 .../exp0/landmarks_evolution.pdf | Bin .../exp0/makefigs_exp0.R | 0 .../exp1-1D/an_exp1-1D.jl | 0 .../exp1-1D/data_exp1-1D.jld | Bin .../exp1-1D/gen_exp1-1D.jl | 0 .../exp1-1D/makefigs_exp1-1D.R | 0 .../exp1-ckl/an_exp1-ckl.jl | 0 .../exp1-ckl/data_exp1-ckl.jld | Bin .../exp1-ckl/data_exp1ckl.jld | Bin .../exp1-ckl/gen_exp1-ckl.jl | 0 .../exp1-ckl/makefigs_exp1-ckl.R | 0 .../exp1-ckl/ms-parsestimated/acceptance.pdf | Bin .../exp1-ckl/ms-parsestimated/bridges-overlaid.pdf | Bin .../exp1-ckl/ms-parsestimated/noisefields.csv | 0 .../exp1-translated/an_exp1-translated.jl | 0 .../exp1-translated/data_exp-translated.jld | Bin .../exp1-translated/data_exp1-translated.jld | Bin .../exp1-translated/gen_exp1-translated.jl | 0 .../exp1-translated/makefigs_exp1-translated.R | 0 .../exp1-try/an_exp1-try.jl | 0 .../exp1-try/data_exp1-try.jld | Bin .../exp1-try/forwardhalfway.jld | Bin .../exp1-try/gen_exp1-try.jl | 0 .../exp1-try/makefigs_exp1-try.R | 0 {src/experiments => experiments}/exp1/an_exp1.jl | 0 {src/experiments => experiments}/exp1/data_exp1.jld | Bin {src/experiments => experiments}/exp1/gen_exp1.jl | 0 .../exp1/makefigs_exp1.R | 0 {src/experiments => experiments}/exp1/ms/accdf.csv | 0 .../exp1/ms/acceptance.pdf | Bin .../exp1/ms/bridges-faceted.pdf | Bin .../exp1/ms/bridges-overlaid.pdf | Bin {src/experiments => experiments}/exp1/ms/info.txt | 0 .../exp1/ms/iterates.csv | 0 .../exp1/ms/momenta-faceted.pdf | Bin .../exp1/ms/noisefields.csv | 0 {src/experiments => experiments}/exp1/ms/obs0.csv | 0 .../exp1/ms/scatter-pars.pdf | Bin .../exp1/ms/trace-pars.pdf | Bin {src/experiments => experiments}/exp2/accdf.csv | 0 .../experiments => experiments}/exp2/acceptance.pdf | Bin {src/experiments => experiments}/exp2/an_exp2.jl | 0 .../exp2/bridges-faceted.pdf | Bin .../exp2/bridges-overlaid.pdf | Bin {src/experiments => experiments}/exp2/data_exp2.jld | Bin {src/experiments => experiments}/exp2/gen_exp2.jl | 0 {src/experiments => experiments}/exp2/info.txt | 0 {src/experiments => experiments}/exp2/iterates.csv | 0 .../exp2/makefigs_exp2.R | 0 .../exp2/momenta-faceted.pdf | Bin .../prior/prior_initpos.jl | 0 63 files changed, 0 insertions(+), 0 deletions(-) rename {src/experiments => experiments}/exp-cardiac/an_cardiac.jl (100%) rename {src/experiments => experiments}/exp-cardiac/cardiac.npy (100%) rename {src/experiments => experiments}/exp-cardiac/data_cardiac.jld (100%) rename {src/experiments => experiments}/exp-cardiac/gen_cardiac.jl (100%) rename {src/experiments => experiments}/exp-cardiac/makefigs_cardiac.R (100%) rename {src/experiments => experiments}/exp-dgs/an_exp-dgs.jl (100%) rename {src/experiments => experiments}/exp-dgs/data_exp-dgs.jld (100%) rename {src/experiments => experiments}/exp-dgs/gen_exp-dgs.jl (100%) rename {src/experiments => experiments}/exp-dgs/makefigs_exp-dgs.R (100%) rename {src/experiments => experiments}/exp-dgs/match.npy.npz (100%) rename {src/experiments => experiments}/exp0/animate_introexample.mp4 (100%) rename {src/experiments => experiments}/exp0/exp0.jl (100%) rename {src/experiments => experiments}/exp0/landmarks_evolution.pdf (100%) rename {src/experiments => experiments}/exp0/makefigs_exp0.R (100%) rename {src/experiments => experiments}/exp1-1D/an_exp1-1D.jl (100%) rename {src/experiments => experiments}/exp1-1D/data_exp1-1D.jld (100%) rename {src/experiments => experiments}/exp1-1D/gen_exp1-1D.jl (100%) rename {src/experiments => experiments}/exp1-1D/makefigs_exp1-1D.R (100%) rename {src/experiments => experiments}/exp1-ckl/an_exp1-ckl.jl (100%) rename {src/experiments => experiments}/exp1-ckl/data_exp1-ckl.jld (100%) rename {src/experiments => experiments}/exp1-ckl/data_exp1ckl.jld (100%) rename {src/experiments => experiments}/exp1-ckl/gen_exp1-ckl.jl (100%) rename {src/experiments => experiments}/exp1-ckl/makefigs_exp1-ckl.R (100%) rename {src/experiments => experiments}/exp1-ckl/ms-parsestimated/acceptance.pdf (100%) rename {src/experiments => experiments}/exp1-ckl/ms-parsestimated/bridges-overlaid.pdf (100%) rename {src/experiments => experiments}/exp1-ckl/ms-parsestimated/noisefields.csv (100%) rename {src/experiments => experiments}/exp1-translated/an_exp1-translated.jl (100%) rename {src/experiments => experiments}/exp1-translated/data_exp-translated.jld (100%) rename {src/experiments => experiments}/exp1-translated/data_exp1-translated.jld (100%) rename {src/experiments => experiments}/exp1-translated/gen_exp1-translated.jl (100%) rename {src/experiments => experiments}/exp1-translated/makefigs_exp1-translated.R (100%) rename {src/experiments => experiments}/exp1-try/an_exp1-try.jl (100%) rename {src/experiments => experiments}/exp1-try/data_exp1-try.jld (100%) rename {src/experiments => experiments}/exp1-try/forwardhalfway.jld (100%) rename {src/experiments => experiments}/exp1-try/gen_exp1-try.jl (100%) rename {src/experiments => experiments}/exp1-try/makefigs_exp1-try.R (100%) rename {src/experiments => experiments}/exp1/an_exp1.jl (100%) rename {src/experiments => experiments}/exp1/data_exp1.jld (100%) rename {src/experiments => experiments}/exp1/gen_exp1.jl (100%) rename {src/experiments => experiments}/exp1/makefigs_exp1.R (100%) rename {src/experiments => experiments}/exp1/ms/accdf.csv (100%) rename {src/experiments => experiments}/exp1/ms/acceptance.pdf (100%) rename {src/experiments => experiments}/exp1/ms/bridges-faceted.pdf (100%) rename {src/experiments => experiments}/exp1/ms/bridges-overlaid.pdf (100%) rename {src/experiments => experiments}/exp1/ms/info.txt (100%) rename {src/experiments => experiments}/exp1/ms/iterates.csv (100%) rename {src/experiments => experiments}/exp1/ms/momenta-faceted.pdf (100%) rename {src/experiments => experiments}/exp1/ms/noisefields.csv (100%) rename {src/experiments => experiments}/exp1/ms/obs0.csv (100%) rename {src/experiments => experiments}/exp1/ms/scatter-pars.pdf (100%) rename {src/experiments => experiments}/exp1/ms/trace-pars.pdf (100%) rename {src/experiments => experiments}/exp2/accdf.csv (100%) rename {src/experiments => experiments}/exp2/acceptance.pdf (100%) rename {src/experiments => experiments}/exp2/an_exp2.jl (100%) rename {src/experiments => experiments}/exp2/bridges-faceted.pdf (100%) rename {src/experiments => experiments}/exp2/bridges-overlaid.pdf (100%) rename {src/experiments => experiments}/exp2/data_exp2.jld (100%) rename {src/experiments => experiments}/exp2/gen_exp2.jl (100%) rename {src/experiments => experiments}/exp2/info.txt (100%) rename {src/experiments => experiments}/exp2/iterates.csv (100%) rename {src/experiments => experiments}/exp2/makefigs_exp2.R (100%) rename {src/experiments => experiments}/exp2/momenta-faceted.pdf (100%) rename {src/experiments => experiments}/prior/prior_initpos.jl (100%) diff --git a/src/experiments/exp-cardiac/an_cardiac.jl b/experiments/exp-cardiac/an_cardiac.jl similarity index 100% rename from src/experiments/exp-cardiac/an_cardiac.jl rename to experiments/exp-cardiac/an_cardiac.jl diff --git a/src/experiments/exp-cardiac/cardiac.npy b/experiments/exp-cardiac/cardiac.npy similarity index 100% rename from src/experiments/exp-cardiac/cardiac.npy rename to experiments/exp-cardiac/cardiac.npy diff --git a/src/experiments/exp-cardiac/data_cardiac.jld b/experiments/exp-cardiac/data_cardiac.jld similarity index 100% rename from src/experiments/exp-cardiac/data_cardiac.jld rename to experiments/exp-cardiac/data_cardiac.jld diff --git a/src/experiments/exp-cardiac/gen_cardiac.jl b/experiments/exp-cardiac/gen_cardiac.jl similarity index 100% rename from src/experiments/exp-cardiac/gen_cardiac.jl rename to experiments/exp-cardiac/gen_cardiac.jl diff --git a/src/experiments/exp-cardiac/makefigs_cardiac.R b/experiments/exp-cardiac/makefigs_cardiac.R similarity index 100% rename from src/experiments/exp-cardiac/makefigs_cardiac.R rename to experiments/exp-cardiac/makefigs_cardiac.R diff --git a/src/experiments/exp-dgs/an_exp-dgs.jl b/experiments/exp-dgs/an_exp-dgs.jl similarity index 100% rename from src/experiments/exp-dgs/an_exp-dgs.jl rename to experiments/exp-dgs/an_exp-dgs.jl diff --git a/src/experiments/exp-dgs/data_exp-dgs.jld b/experiments/exp-dgs/data_exp-dgs.jld similarity index 100% rename from src/experiments/exp-dgs/data_exp-dgs.jld rename to experiments/exp-dgs/data_exp-dgs.jld diff --git a/src/experiments/exp-dgs/gen_exp-dgs.jl b/experiments/exp-dgs/gen_exp-dgs.jl similarity index 100% rename from src/experiments/exp-dgs/gen_exp-dgs.jl rename to experiments/exp-dgs/gen_exp-dgs.jl diff --git a/src/experiments/exp-dgs/makefigs_exp-dgs.R b/experiments/exp-dgs/makefigs_exp-dgs.R similarity index 100% rename from src/experiments/exp-dgs/makefigs_exp-dgs.R rename to experiments/exp-dgs/makefigs_exp-dgs.R diff --git a/src/experiments/exp-dgs/match.npy.npz b/experiments/exp-dgs/match.npy.npz similarity index 100% rename from src/experiments/exp-dgs/match.npy.npz rename to experiments/exp-dgs/match.npy.npz diff --git a/src/experiments/exp0/animate_introexample.mp4 b/experiments/exp0/animate_introexample.mp4 similarity index 100% rename from src/experiments/exp0/animate_introexample.mp4 rename to experiments/exp0/animate_introexample.mp4 diff --git a/src/experiments/exp0/exp0.jl b/experiments/exp0/exp0.jl similarity index 100% rename from src/experiments/exp0/exp0.jl rename to experiments/exp0/exp0.jl diff --git a/src/experiments/exp0/landmarks_evolution.pdf b/experiments/exp0/landmarks_evolution.pdf similarity index 100% rename from src/experiments/exp0/landmarks_evolution.pdf rename to experiments/exp0/landmarks_evolution.pdf diff --git a/src/experiments/exp0/makefigs_exp0.R b/experiments/exp0/makefigs_exp0.R similarity index 100% rename from src/experiments/exp0/makefigs_exp0.R rename to experiments/exp0/makefigs_exp0.R diff --git a/src/experiments/exp1-1D/an_exp1-1D.jl b/experiments/exp1-1D/an_exp1-1D.jl similarity index 100% rename from src/experiments/exp1-1D/an_exp1-1D.jl rename to experiments/exp1-1D/an_exp1-1D.jl diff --git a/src/experiments/exp1-1D/data_exp1-1D.jld b/experiments/exp1-1D/data_exp1-1D.jld similarity index 100% rename from src/experiments/exp1-1D/data_exp1-1D.jld rename to experiments/exp1-1D/data_exp1-1D.jld diff --git a/src/experiments/exp1-1D/gen_exp1-1D.jl b/experiments/exp1-1D/gen_exp1-1D.jl similarity index 100% rename from src/experiments/exp1-1D/gen_exp1-1D.jl rename to experiments/exp1-1D/gen_exp1-1D.jl diff --git a/src/experiments/exp1-1D/makefigs_exp1-1D.R b/experiments/exp1-1D/makefigs_exp1-1D.R similarity index 100% rename from src/experiments/exp1-1D/makefigs_exp1-1D.R rename to experiments/exp1-1D/makefigs_exp1-1D.R diff --git a/src/experiments/exp1-ckl/an_exp1-ckl.jl b/experiments/exp1-ckl/an_exp1-ckl.jl similarity index 100% rename from src/experiments/exp1-ckl/an_exp1-ckl.jl rename to experiments/exp1-ckl/an_exp1-ckl.jl diff --git a/src/experiments/exp1-ckl/data_exp1-ckl.jld b/experiments/exp1-ckl/data_exp1-ckl.jld similarity index 100% rename from src/experiments/exp1-ckl/data_exp1-ckl.jld rename to experiments/exp1-ckl/data_exp1-ckl.jld diff --git a/src/experiments/exp1-ckl/data_exp1ckl.jld b/experiments/exp1-ckl/data_exp1ckl.jld similarity index 100% rename from src/experiments/exp1-ckl/data_exp1ckl.jld rename to experiments/exp1-ckl/data_exp1ckl.jld diff --git a/src/experiments/exp1-ckl/gen_exp1-ckl.jl b/experiments/exp1-ckl/gen_exp1-ckl.jl similarity index 100% rename from src/experiments/exp1-ckl/gen_exp1-ckl.jl rename to experiments/exp1-ckl/gen_exp1-ckl.jl diff --git a/src/experiments/exp1-ckl/makefigs_exp1-ckl.R b/experiments/exp1-ckl/makefigs_exp1-ckl.R similarity index 100% rename from src/experiments/exp1-ckl/makefigs_exp1-ckl.R rename to experiments/exp1-ckl/makefigs_exp1-ckl.R diff --git a/src/experiments/exp1-ckl/ms-parsestimated/acceptance.pdf b/experiments/exp1-ckl/ms-parsestimated/acceptance.pdf similarity index 100% rename from src/experiments/exp1-ckl/ms-parsestimated/acceptance.pdf rename to experiments/exp1-ckl/ms-parsestimated/acceptance.pdf diff --git a/src/experiments/exp1-ckl/ms-parsestimated/bridges-overlaid.pdf b/experiments/exp1-ckl/ms-parsestimated/bridges-overlaid.pdf similarity index 100% rename from src/experiments/exp1-ckl/ms-parsestimated/bridges-overlaid.pdf rename to experiments/exp1-ckl/ms-parsestimated/bridges-overlaid.pdf diff --git a/src/experiments/exp1-ckl/ms-parsestimated/noisefields.csv b/experiments/exp1-ckl/ms-parsestimated/noisefields.csv similarity index 100% rename from src/experiments/exp1-ckl/ms-parsestimated/noisefields.csv rename to experiments/exp1-ckl/ms-parsestimated/noisefields.csv diff --git a/src/experiments/exp1-translated/an_exp1-translated.jl b/experiments/exp1-translated/an_exp1-translated.jl similarity index 100% rename from src/experiments/exp1-translated/an_exp1-translated.jl rename to experiments/exp1-translated/an_exp1-translated.jl diff --git a/src/experiments/exp1-translated/data_exp-translated.jld b/experiments/exp1-translated/data_exp-translated.jld similarity index 100% rename from src/experiments/exp1-translated/data_exp-translated.jld rename to experiments/exp1-translated/data_exp-translated.jld diff --git a/src/experiments/exp1-translated/data_exp1-translated.jld b/experiments/exp1-translated/data_exp1-translated.jld similarity index 100% rename from src/experiments/exp1-translated/data_exp1-translated.jld rename to experiments/exp1-translated/data_exp1-translated.jld diff --git a/src/experiments/exp1-translated/gen_exp1-translated.jl b/experiments/exp1-translated/gen_exp1-translated.jl similarity index 100% rename from src/experiments/exp1-translated/gen_exp1-translated.jl rename to experiments/exp1-translated/gen_exp1-translated.jl diff --git a/src/experiments/exp1-translated/makefigs_exp1-translated.R b/experiments/exp1-translated/makefigs_exp1-translated.R similarity index 100% rename from src/experiments/exp1-translated/makefigs_exp1-translated.R rename to experiments/exp1-translated/makefigs_exp1-translated.R diff --git a/src/experiments/exp1-try/an_exp1-try.jl b/experiments/exp1-try/an_exp1-try.jl similarity index 100% rename from src/experiments/exp1-try/an_exp1-try.jl rename to experiments/exp1-try/an_exp1-try.jl diff --git a/src/experiments/exp1-try/data_exp1-try.jld b/experiments/exp1-try/data_exp1-try.jld similarity index 100% rename from src/experiments/exp1-try/data_exp1-try.jld rename to experiments/exp1-try/data_exp1-try.jld diff --git a/src/experiments/exp1-try/forwardhalfway.jld b/experiments/exp1-try/forwardhalfway.jld similarity index 100% rename from src/experiments/exp1-try/forwardhalfway.jld rename to experiments/exp1-try/forwardhalfway.jld diff --git a/src/experiments/exp1-try/gen_exp1-try.jl b/experiments/exp1-try/gen_exp1-try.jl similarity index 100% rename from src/experiments/exp1-try/gen_exp1-try.jl rename to experiments/exp1-try/gen_exp1-try.jl diff --git a/src/experiments/exp1-try/makefigs_exp1-try.R b/experiments/exp1-try/makefigs_exp1-try.R similarity index 100% rename from src/experiments/exp1-try/makefigs_exp1-try.R rename to experiments/exp1-try/makefigs_exp1-try.R diff --git a/src/experiments/exp1/an_exp1.jl b/experiments/exp1/an_exp1.jl similarity index 100% rename from src/experiments/exp1/an_exp1.jl rename to experiments/exp1/an_exp1.jl diff --git a/src/experiments/exp1/data_exp1.jld b/experiments/exp1/data_exp1.jld similarity index 100% rename from src/experiments/exp1/data_exp1.jld rename to experiments/exp1/data_exp1.jld diff --git a/src/experiments/exp1/gen_exp1.jl b/experiments/exp1/gen_exp1.jl similarity index 100% rename from src/experiments/exp1/gen_exp1.jl rename to experiments/exp1/gen_exp1.jl diff --git a/src/experiments/exp1/makefigs_exp1.R b/experiments/exp1/makefigs_exp1.R similarity index 100% rename from src/experiments/exp1/makefigs_exp1.R rename to experiments/exp1/makefigs_exp1.R diff --git a/src/experiments/exp1/ms/accdf.csv b/experiments/exp1/ms/accdf.csv similarity index 100% rename from src/experiments/exp1/ms/accdf.csv rename to experiments/exp1/ms/accdf.csv diff --git a/src/experiments/exp1/ms/acceptance.pdf b/experiments/exp1/ms/acceptance.pdf similarity index 100% rename from src/experiments/exp1/ms/acceptance.pdf rename to experiments/exp1/ms/acceptance.pdf diff --git a/src/experiments/exp1/ms/bridges-faceted.pdf b/experiments/exp1/ms/bridges-faceted.pdf similarity index 100% rename from src/experiments/exp1/ms/bridges-faceted.pdf rename to experiments/exp1/ms/bridges-faceted.pdf diff --git a/src/experiments/exp1/ms/bridges-overlaid.pdf b/experiments/exp1/ms/bridges-overlaid.pdf similarity index 100% rename from src/experiments/exp1/ms/bridges-overlaid.pdf rename to experiments/exp1/ms/bridges-overlaid.pdf diff --git a/src/experiments/exp1/ms/info.txt b/experiments/exp1/ms/info.txt similarity index 100% rename from src/experiments/exp1/ms/info.txt rename to experiments/exp1/ms/info.txt diff --git a/src/experiments/exp1/ms/iterates.csv b/experiments/exp1/ms/iterates.csv similarity index 100% rename from src/experiments/exp1/ms/iterates.csv rename to experiments/exp1/ms/iterates.csv diff --git a/src/experiments/exp1/ms/momenta-faceted.pdf b/experiments/exp1/ms/momenta-faceted.pdf similarity index 100% rename from src/experiments/exp1/ms/momenta-faceted.pdf rename to experiments/exp1/ms/momenta-faceted.pdf diff --git a/src/experiments/exp1/ms/noisefields.csv b/experiments/exp1/ms/noisefields.csv similarity index 100% rename from src/experiments/exp1/ms/noisefields.csv rename to experiments/exp1/ms/noisefields.csv diff --git a/src/experiments/exp1/ms/obs0.csv b/experiments/exp1/ms/obs0.csv similarity index 100% rename from src/experiments/exp1/ms/obs0.csv rename to experiments/exp1/ms/obs0.csv diff --git a/src/experiments/exp1/ms/scatter-pars.pdf b/experiments/exp1/ms/scatter-pars.pdf similarity index 100% rename from src/experiments/exp1/ms/scatter-pars.pdf rename to experiments/exp1/ms/scatter-pars.pdf diff --git a/src/experiments/exp1/ms/trace-pars.pdf b/experiments/exp1/ms/trace-pars.pdf similarity index 100% rename from src/experiments/exp1/ms/trace-pars.pdf rename to experiments/exp1/ms/trace-pars.pdf diff --git a/src/experiments/exp2/accdf.csv b/experiments/exp2/accdf.csv similarity index 100% rename from src/experiments/exp2/accdf.csv rename to experiments/exp2/accdf.csv diff --git a/src/experiments/exp2/acceptance.pdf b/experiments/exp2/acceptance.pdf similarity index 100% rename from src/experiments/exp2/acceptance.pdf rename to experiments/exp2/acceptance.pdf diff --git a/src/experiments/exp2/an_exp2.jl b/experiments/exp2/an_exp2.jl similarity index 100% rename from src/experiments/exp2/an_exp2.jl rename to experiments/exp2/an_exp2.jl diff --git a/src/experiments/exp2/bridges-faceted.pdf b/experiments/exp2/bridges-faceted.pdf similarity index 100% rename from src/experiments/exp2/bridges-faceted.pdf rename to experiments/exp2/bridges-faceted.pdf diff --git a/src/experiments/exp2/bridges-overlaid.pdf b/experiments/exp2/bridges-overlaid.pdf similarity index 100% rename from src/experiments/exp2/bridges-overlaid.pdf rename to experiments/exp2/bridges-overlaid.pdf diff --git a/src/experiments/exp2/data_exp2.jld b/experiments/exp2/data_exp2.jld similarity index 100% rename from src/experiments/exp2/data_exp2.jld rename to experiments/exp2/data_exp2.jld diff --git a/src/experiments/exp2/gen_exp2.jl b/experiments/exp2/gen_exp2.jl similarity index 100% rename from src/experiments/exp2/gen_exp2.jl rename to experiments/exp2/gen_exp2.jl diff --git a/src/experiments/exp2/info.txt b/experiments/exp2/info.txt similarity index 100% rename from src/experiments/exp2/info.txt rename to experiments/exp2/info.txt diff --git a/src/experiments/exp2/iterates.csv b/experiments/exp2/iterates.csv similarity index 100% rename from src/experiments/exp2/iterates.csv rename to experiments/exp2/iterates.csv diff --git a/src/experiments/exp2/makefigs_exp2.R b/experiments/exp2/makefigs_exp2.R similarity index 100% rename from src/experiments/exp2/makefigs_exp2.R rename to experiments/exp2/makefigs_exp2.R diff --git a/src/experiments/exp2/momenta-faceted.pdf b/experiments/exp2/momenta-faceted.pdf similarity index 100% rename from src/experiments/exp2/momenta-faceted.pdf rename to experiments/exp2/momenta-faceted.pdf diff --git a/src/experiments/prior/prior_initpos.jl b/experiments/prior/prior_initpos.jl similarity index 100% rename from src/experiments/prior/prior_initpos.jl rename to experiments/prior/prior_initpos.jl From a64313c98d7278c167fdac21892d57de4376eaee Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 12:24:47 +0100 Subject: [PATCH 10/20] Move postprocessing --- scripts/postprocessing.jl | 193 +++++++++++++++++++++++--------------- src/postprocessing.jl | 135 -------------------------- 2 files changed, 115 insertions(+), 213 deletions(-) delete mode 100644 src/postprocessing.jl diff --git a/scripts/postprocessing.jl b/scripts/postprocessing.jl index e070419..6392eef 100644 --- a/scripts/postprocessing.jl +++ b/scripts/postprocessing.jl @@ -1,98 +1,135 @@ -## postprocessing +""" + Write bridge iterates to file + Xsave: contains values of bridges + tt_: grid on which bridges are simulated + n: nr of bridges + nshapes: nr of shapes + subsamples: indices of iterates that are kept and saved in Xsave + outdir: output directory +""" +function write_mcmc_iterates(Xsave, tt_, n, nshapes, subsamples, outdir) + nshapes = length(xobsT) + iterates = reshape(vcat(Xsave...),2*d*length(tt_)*n*nshapes, length(subsamples)) # each column contains samplepath of an iteration -fn = "_" * string(model) * "_" * string(sampler) *"_" * string(dataset) -gif(anim, outdir*fn*".gif", fps = 50) -mp4(anim, outdir*fn*".mp4", fps = 50) + # total number of rows is nshapes * length(tt_) * n * (2d) -# make a fig for acceptance probs in parameter and initial state updating -accdf = DataFrame(kernel = map(x->x.kernel, accinfo), acc = map(x->x.acc, accinfo), iter = 1:length(accinfo)) -@rput accdf -@rput outdir -R""" -library(tidyverse) -library(ggplot2) -theme_set(theme_light()) -p <- accdf %>% mutate(kernel=as.character(kernel)) %>% - ggplot(aes(x=iter, y=acc)) + geom_point() + - facet_wrap(~kernel) -ggsave(paste0(outdir,"acceptance.pdf"),p) -""" + # Ordering in each column is as follows: + # 0) shape + # 1) time + # 2) landmark nr + # 3) for each landmark: q1, q2 p1, p2 + #if !(d==2) error("pqtype in write_mcmc_iterates only implemented in case d=2") end + shapes = repeat(1:nshapes, inner=length(tt_)*2d*n) + #times = repeat(tt_,inner=2d*n*nshapes) # original buggy version when nshapes=1 + times = repeat(tt_, inner=2*d*n, outer=nshapes) + landmarkid = repeat(1:n, inner=2d, outer=length(tt_)*nshapes) + if d==1 + pqtype = repeat(["pos1", "mom1"], outer=length(tt_)*n*nshapes) + elseif d==2 + pqtype = repeat(["pos1", "pos2", "mom1", "mom2"], outer=length(tt_)*n*nshapes) + elseif d==3 + pqtype = repeat(["pos1", "pos2", "pos3", "mom1", "mom2", "mom3"], outer=length(tt_)*n*nshapes) + end + out = hcat(times,pqtype,landmarkid,shapes,iterates) + headline = "time " * "pqtype " * "landmarkid " * "shapes " * prod(map(x -> "iter"*string(x)*" ",subsamples)) + headline = chop(headline,tail=1) * "\n" + fn = outdir*"iterates.csv" + f = open(fn,"w") + write(f, headline) + writedlm(f,out) + close(f) +end -######### write mcmc iterates of bridges to csv file -nshapes = length(xobsT) -iterates = reshape(vcat(Xsave...),2*d*length(tt_)*P.n*nshapes, length(subsamples)) # each column contains samplepath of an iteration -# Ordering in each column is as follows: -# 0) shape -# 1) time -# 2) landmark nr -# 3) for each landmark: q1, q2 p1, p2 -pqtype = repeat(["pos1", "pos2", "mom1", "mom2"], length(tt_)*P.n*nshapes) -times = repeat(tt_,inner=2d*P.n*nshapes) -landmarkid = repeat(1:P.n, inner=2d, outer=length(tt_)*nshapes) -shapes = repeat(1:nshapes, inner=length(tt_)*2d*P.n) -out = hcat(times,pqtype,landmarkid,shapes,iterates) -headline = "time " * "pqtype " * "landmarkid " * "shapes " * prod(map(x -> "iter"*string(x)*" ",subsamples)) -headline = chop(headline,tail=1) * "\n" +""" + Write info to txt file +""" +function write_info(model,ITER, n, tt_, updatescheme, Σobs, tp, ρ, δ, perc_acc_pcn, elapsed, outdir) + f = open(outdir*"info.txt","w") + write(f, "Model: ", string(model),"\n") + write(f, "Number of iterations: ",string(ITER),"\n") + write(f, "Number of landmarks: ",string(n),"\n") + write(f, "Length time grid: ", string(length(tt_)),"\n") + write(f, "Endpoint: ",string(tt_[end]),"\n") + write(f, "updatescheme: ", string(updatescheme),"\n") + write(f, "Noise Sigma: ",string(Σobs),"\n") + write(f, "tuningpars_mcmc: ", string(tp),"\n") + write(f, "Final value of rho (Crank-Nicholsen parameter: ",string(ρ),"\n") + write(f, "Final value of MALA parameter (delta): ",string(δ),"\n") -fn = outdir*"iterates.csv" -f = open(fn,"w") -write(f, headline) -writedlm(f,out) -close(f) + write(f, "skip in evaluation of loglikelihood: ",string(sk),"\n") + write(f, "Average acceptance percentage pCN update steps: ",string(perc_acc_pcn),"\n\n") + write(f, "Elapsed time: ", string(elapsed),"\n") + close(f) +end -########### write info to txt file -fn = outdir* "info_" * string(model) * "_" * string(sampler) *"_" * string(dataset)*".txt" -f = open(fn,"w") -write(f, "Dataset: ", string(dataset),"\n") -write(f, "Sampler: ", string(sampler), "\n") +function write_acc(accinfo,accpcn,nshapes,outdir) + # extract number of distinct update steps in accinfo + nunique = length(unique(map(x->x[1], accinfo))) + niterates = div(length(accinfo),nunique) + accdf = DataFrame(kernel = map(x->Symbol(x.kernel), accinfo), acc = map(x->x.acc, accinfo), iter = repeat(1:niterates, inner= nunique)) + accpcndf = DataFrame(kernel = fill(Symbol("pCN"),length(accpcn)), acc=accpcn, iter=repeat(1:niterates,inner=nshapes)) + append!(accdf, accpcndf) + CSV.write(outdir*"accdf.csv", accdf; delim=";") +end -write(f, "Number of iterations: ",string(ITER),"\n") -write(f, "Number of landmarks: ",string(P.n),"\n") -write(f, "Length time grid: ", string(length(tt_)),"\n") -write(f, "Mesh width: ",string(dt),"\n") -write(f, "Noise Sigma: ",string(σobs),"\n") -write(f, "rho (Crank-Nicholsen parameter: ",string(ρ),"\n") -write(f, "MALA parameter (delta): ",string(δ),"\n") -write(f, "skip in evaluation of loglikelihood: ",string(LM.sk),"\n") -write(f, "Average acceptance percentage pCN update steps: ",string(perc_acc_pcn),"\n\n") -#write(f, "Backward type parametrisation in terms of nu and H? ",string(νHparam),"\n") -close(f) -######## write observations to file -# if obs_atzero -# obsdf = DataFrame(x=vcat( extractcomp(xobs0,1), extractcomp(xobsT,1)), -# y= vcat( extractcomp(xobs0,2), extractcomp(xobsT,2)), -# time=repeat(["0","T"], inner=P.n)) -# CSV.write(outdir*"observations.csv", obsdf; delim=";") -# else +""" + Write observations to file +""" +function write_observations(xobs0, xobsT, n, nshapes, x0,outdir) valueT = vcat(map(x->deepvec(x), xobsT)...) # merge all observations at time T in one vector - posT = repeat(["pos1","pos2"], P.n*nshapes) - shT = repeat(1:nshapes, inner=d*P.n) - obsTdf = DataFrame(pos=posT,shape=shT, value=valueT,landmark=repeat(1:P.n,inner=d,outer=nshapes)) + if d==1 + posT = repeat(["pos1"], n*nshapes) + elseif d==2 + posT = repeat(["pos1","pos2"], n*nshapes) + elseif d==3 + posT = repeat(["pos1","pos2","pos3"], n*nshapes) + end + shT = repeat(1:nshapes, inner=d*n) + obsTdf = DataFrame(pos=posT,shape=shT, value=valueT,landmark=repeat(1:n,inner=d,outer=nshapes)) q0 = map(x->vec(x),x0.q) p0 = map(x->vec(x),x0.p) - obs0df = DataFrame(pos1=extractcomp(q0,1), pos2=extractcomp(q0,2), mom1=extractcomp(p0,1) , mom2=extractcomp(p0,2),landmark=1:P.n) - + if d==1 + obs0df = DataFrame(pos1=extractcomp(q0,1), mom1=extractcomp(p0,1),landmark=1:n) + elseif d==2 + obs0df = DataFrame(pos1=extractcomp(q0,1), pos2=extractcomp(q0,2), mom1=extractcomp(p0,1) , mom2=extractcomp(p0,2),landmark=1:n) + elseif d==3 + obs0df = DataFrame(pos1=extractcomp(q0,1), pos2=extractcomp(q0,2), pos3=extractcomp(q0,3), + mom1=extractcomp(p0,1) , mom2=extractcomp(p0,2), mom3=extractcomp(p0,3),landmark=1:n) + end CSV.write(outdir*"obs0.csv", obs0df; delim=";") CSV.write(outdir*"obsT.csv", obsTdf; delim=";") -# end +end -######## write parameter iterates to file -parsdf = DataFrame(a=extractcomp(parsave,1),c=extractcomp(parsave,2), +""" + Write parameter iterates to file +""" +function write_params(parsave,subsamples,outdir) + parsdf = DataFrame(a=extractcomp(parsave,1),c=extractcomp(parsave,2), gamma=extractcomp(parsave,3), iterate=subsamples) -CSV.write(outdir*"parameters.csv", parsdf; delim=";") + CSV.write(outdir*"parameters.csv", parsdf; delim=";") +end -####### write noisefields to file -if isa(P,Landmarks) - nfsloc = [P.nfs[j].δ for j in eachindex(P.nfs)] - nfsdf = DataFrame(locx = extractcomp(nfsloc,1), - locy = extractcomp(nfsloc,2), - nfstd=fill(P.nfstd,length(P.nfs))) -elseif isa(P,MarslandShardlow) - nfsdf =DataFrame(locx=Int64[], locy=Int64[], nfstd=Int64[]) +""" + Write noisefields to file (empty in case of MS-model) +""" +function write_noisefields(P,outdir) + if isa(P,Landmarks) + nfsloc = [P.nfs[j].δ for j in eachindex(P.nfs)] + if d==2 + nfsdf = DataFrame(locx = extractcomp(nfsloc,1), + locy = extractcomp(nfsloc,2), + nfstd=fill(P.nfstd,length(P.nfs))) + elseif d==1 + nfsdf = DataFrame(locx = extractcomp(nfsloc,1), + nfstd=fill(P.nfstd,length(P.nfs))) + end + elseif isa(P,MarslandShardlow) + nfsdf =DataFrame(locx=Int64[], locy=Int64[], nfstd=Int64[]) + end + CSV.write(outdir*"noisefields.csv", nfsdf; delim=";") end -CSV.write(outdir*"noisefields.csv", nfsdf; delim=";") diff --git a/src/postprocessing.jl b/src/postprocessing.jl deleted file mode 100644 index 6392eef..0000000 --- a/src/postprocessing.jl +++ /dev/null @@ -1,135 +0,0 @@ -""" - Write bridge iterates to file - Xsave: contains values of bridges - tt_: grid on which bridges are simulated - n: nr of bridges - nshapes: nr of shapes - subsamples: indices of iterates that are kept and saved in Xsave - outdir: output directory -""" -function write_mcmc_iterates(Xsave, tt_, n, nshapes, subsamples, outdir) - nshapes = length(xobsT) - iterates = reshape(vcat(Xsave...),2*d*length(tt_)*n*nshapes, length(subsamples)) # each column contains samplepath of an iteration - - # total number of rows is nshapes * length(tt_) * n * (2d) - - # Ordering in each column is as follows: - # 0) shape - # 1) time - # 2) landmark nr - # 3) for each landmark: q1, q2 p1, p2 - #if !(d==2) error("pqtype in write_mcmc_iterates only implemented in case d=2") end - shapes = repeat(1:nshapes, inner=length(tt_)*2d*n) - #times = repeat(tt_,inner=2d*n*nshapes) # original buggy version when nshapes=1 - times = repeat(tt_, inner=2*d*n, outer=nshapes) - landmarkid = repeat(1:n, inner=2d, outer=length(tt_)*nshapes) - if d==1 - pqtype = repeat(["pos1", "mom1"], outer=length(tt_)*n*nshapes) - elseif d==2 - pqtype = repeat(["pos1", "pos2", "mom1", "mom2"], outer=length(tt_)*n*nshapes) - elseif d==3 - pqtype = repeat(["pos1", "pos2", "pos3", "mom1", "mom2", "mom3"], outer=length(tt_)*n*nshapes) - end - out = hcat(times,pqtype,landmarkid,shapes,iterates) - headline = "time " * "pqtype " * "landmarkid " * "shapes " * prod(map(x -> "iter"*string(x)*" ",subsamples)) - headline = chop(headline,tail=1) * "\n" - - fn = outdir*"iterates.csv" - f = open(fn,"w") - write(f, headline) - writedlm(f,out) - close(f) -end - - -""" - Write info to txt file -""" -function write_info(model,ITER, n, tt_, updatescheme, Σobs, tp, ρ, δ, perc_acc_pcn, elapsed, outdir) - f = open(outdir*"info.txt","w") - write(f, "Model: ", string(model),"\n") - write(f, "Number of iterations: ",string(ITER),"\n") - write(f, "Number of landmarks: ",string(n),"\n") - write(f, "Length time grid: ", string(length(tt_)),"\n") - write(f, "Endpoint: ",string(tt_[end]),"\n") - write(f, "updatescheme: ", string(updatescheme),"\n") - write(f, "Noise Sigma: ",string(Σobs),"\n") - write(f, "tuningpars_mcmc: ", string(tp),"\n") - write(f, "Final value of rho (Crank-Nicholsen parameter: ",string(ρ),"\n") - write(f, "Final value of MALA parameter (delta): ",string(δ),"\n") - - write(f, "skip in evaluation of loglikelihood: ",string(sk),"\n") - write(f, "Average acceptance percentage pCN update steps: ",string(perc_acc_pcn),"\n\n") - write(f, "Elapsed time: ", string(elapsed),"\n") - close(f) -end - -function write_acc(accinfo,accpcn,nshapes,outdir) - # extract number of distinct update steps in accinfo - nunique = length(unique(map(x->x[1], accinfo))) - niterates = div(length(accinfo),nunique) - accdf = DataFrame(kernel = map(x->Symbol(x.kernel), accinfo), acc = map(x->x.acc, accinfo), iter = repeat(1:niterates, inner= nunique)) - accpcndf = DataFrame(kernel = fill(Symbol("pCN"),length(accpcn)), acc=accpcn, iter=repeat(1:niterates,inner=nshapes)) - append!(accdf, accpcndf) - CSV.write(outdir*"accdf.csv", accdf; delim=";") -end - - -""" - Write observations to file -""" -function write_observations(xobs0, xobsT, n, nshapes, x0,outdir) - valueT = vcat(map(x->deepvec(x), xobsT)...) # merge all observations at time T in one vector - if d==1 - posT = repeat(["pos1"], n*nshapes) - elseif d==2 - posT = repeat(["pos1","pos2"], n*nshapes) - elseif d==3 - posT = repeat(["pos1","pos2","pos3"], n*nshapes) - end - shT = repeat(1:nshapes, inner=d*n) - obsTdf = DataFrame(pos=posT,shape=shT, value=valueT,landmark=repeat(1:n,inner=d,outer=nshapes)) - - q0 = map(x->vec(x),x0.q) - p0 = map(x->vec(x),x0.p) - if d==1 - obs0df = DataFrame(pos1=extractcomp(q0,1), mom1=extractcomp(p0,1),landmark=1:n) - elseif d==2 - obs0df = DataFrame(pos1=extractcomp(q0,1), pos2=extractcomp(q0,2), mom1=extractcomp(p0,1) , mom2=extractcomp(p0,2),landmark=1:n) - elseif d==3 - obs0df = DataFrame(pos1=extractcomp(q0,1), pos2=extractcomp(q0,2), pos3=extractcomp(q0,3), - mom1=extractcomp(p0,1) , mom2=extractcomp(p0,2), mom3=extractcomp(p0,3),landmark=1:n) - end - CSV.write(outdir*"obs0.csv", obs0df; delim=";") - CSV.write(outdir*"obsT.csv", obsTdf; delim=";") -end - - -""" - Write parameter iterates to file -""" -function write_params(parsave,subsamples,outdir) - parsdf = DataFrame(a=extractcomp(parsave,1),c=extractcomp(parsave,2), - gamma=extractcomp(parsave,3), iterate=subsamples) - CSV.write(outdir*"parameters.csv", parsdf; delim=";") -end - -""" - Write noisefields to file (empty in case of MS-model) -""" -function write_noisefields(P,outdir) - if isa(P,Landmarks) - nfsloc = [P.nfs[j].δ for j in eachindex(P.nfs)] - if d==2 - nfsdf = DataFrame(locx = extractcomp(nfsloc,1), - locy = extractcomp(nfsloc,2), - nfstd=fill(P.nfstd,length(P.nfs))) - elseif d==1 - nfsdf = DataFrame(locx = extractcomp(nfsloc,1), - nfstd=fill(P.nfstd,length(P.nfs))) - end - elseif isa(P,MarslandShardlow) - nfsdf =DataFrame(locx=Int64[], locy=Int64[], nfstd=Int64[]) - end - CSV.write(outdir*"noisefields.csv", nfsdf; delim=";") -end From 371d37b1b64be92966e9ba8b966aa9a4b1b83ed3 Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 12:25:42 +0100 Subject: [PATCH 11/20] Mv animation --- {src => scripts}/landmarkbridge_animation.jl | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename {src => scripts}/landmarkbridge_animation.jl (100%) diff --git a/src/landmarkbridge_animation.jl b/scripts/landmarkbridge_animation.jl similarity index 100% rename from src/landmarkbridge_animation.jl rename to scripts/landmarkbridge_animation.jl From e7e06f64fc911cc567f6e7a8031245dd0e5f2da9 Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 12:41:39 +0100 Subject: [PATCH 12/20] data-stefan has moved --- scripts/generatedata.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/scripts/generatedata.jl b/scripts/generatedata.jl index d1f06bd..509125e 100644 --- a/scripts/generatedata.jl +++ b/scripts/generatedata.jl @@ -50,8 +50,8 @@ function generatedata(dataset,P,t,σobs) obs_atzero = true end if dataset=="bear" - bear0 = readdlm(joinpath(@__DIR__,"data-stefan", "bear1.csv"), ',') - bearT = readdlm(joinpath(@__DIR__,"data-stefan", "bear2.csv"), ',') + bear0 = readdlm(joinpath(@__DIR__, "..", "data-stefan", "bear1.csv"), ',') + bearT = readdlm(joinpath(@__DIR__, "..", "data-stefan", "bear2.csv"), ',') nb = size(bear0)[1] avePoint = Point(414.0, 290.0) # average of (x,y)-coords for bear0 to center figure at origin xobs0 = [Point(bear0[i,1], bear0[i,2]) - avePoint for i in 1:nb]/200. @@ -94,7 +94,7 @@ function generatedata(dataset,P,t,σobs) obs_atzero = true end if dataset=="generatedstefan" - testshapes = npzread(joinpath(@__DIR__,"data-stefan", "match.npy.npz")) + testshapes = npzread(joinpath(@__DIR__, "..", "data-stefan", "match.npy.npz")) xobs0vec = get(testshapes,"q0",0) xobsTvec = get(testshapes,"v",0) p0vec = get(testshapes,"p",0) @@ -155,7 +155,7 @@ function generatedata(dataset,P,t,σobs) # for i in range(N_samples): # plt.plot(qs[i,:,0],qs[i,:,1]) - cardiac = npzread(joinpath(@__DIR__,"data-stefan","cardiac.npy")) # heart data (left ventricles, the one we used in https://arxiv.org/abs/1705.10943 + cardiac = npzread(joinpath(@__DIR__, "..", "data-stefan","cardiac.npy")) # heart data (left ventricles, the one we used in https://arxiv.org/abs/1705.10943 cardiacx = cardiac[:,:,1] # x-coordinates of landmarks cardiacy = cardiac[:,:,2] # y-coordinates of landmarks nshapes = 14 @@ -189,7 +189,7 @@ function generatedata(dataset,P,t,σobs) end if false - cc = npzread(joinpath(@__DIR__,"data-stefan","cc.npy")) # corpus callosum data + cc = npzread(joinpath(@__DIR__, "..", "data-stefan","cc.npy")) # corpus callosum data # there are 14 cardiac images, in https://arxiv.org/pdf/1705.10943.pdf 17 landmarks are chosen end From bf1d9e71f83b4770cb0129cfc5c4954a08f459ff Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 13:24:33 +0100 Subject: [PATCH 13/20] Fix exp0 --- experiments/exp0/exp0.jl | 18 ++++++++++-------- scripts/postprocessing.jl | 25 ++++++++++++++++++------- src/BridgeLandmarks.jl | 1 + 3 files changed, 29 insertions(+), 15 deletions(-) diff --git a/experiments/exp0/exp0.jl b/experiments/exp0/exp0.jl index 63cf0ca..9317f28 100644 --- a/experiments/exp0/exp0.jl +++ b/experiments/exp0/exp0.jl @@ -19,8 +19,9 @@ Random.seed!(3) workdir = @__DIR__ cd(workdir) -include(dirname(dirname(workdir))*"/postprocessing.jl") -outdir = workdir*("/") +include(joinpath(BL.dir(),"scripts", "postprocessing.jl")) +outdir = workdir +mkpath(joinpath(outdir, "forward")) model = [:ms, :ahs][1] T = 1.0; dt = 0.005; t = 0.0:dt:T; tt_ = tc(t,T) @@ -55,7 +56,7 @@ xobsT =[ BL.q(Xf.yy[end])] # [Xf.yy[end].q[i] for i in 1:n ] Xfsave = []# typeof(zeros(length(tt_) * P.n * 2 * d * nshapes))[] push!(Xfsave, BL.convert_samplepath(Xf)) -write_mcmc_iterates(Xfsave, tt_, n, nshapes, [1], outdir*"forward/") +write_mcmc_iterates(Xfsave, tt_, n, nshapes, [1], joinpath(outdir,"forward")) #--------- MCMC tuning pars --------------------------------------------------------- @@ -91,8 +92,9 @@ xobs0 = BL.q(x0) #priorθ = product_distribution(fill(Exponential(1.0),3)) priorθ = product_distribution([Exponential(ainit), Exponential(cinit), Exponential(γinit)]) κ = 100.0 -prior_momenta = MvNormalCanon(gramkernel(xobs0,P)/κ) -logpriormom(x0) = logpdf(prior_momenta, vcat(BL.p(x0)...))# +logpdf(prior_positions, vcat(BL.q(x0)...)) +priormom = prior_momenta = MvNormalCanon(gramkernel(xobs0,P)/κ) +#logpriormom(x0) = logpdf(prior_momenta, vcat(BL.p(x0)...))# +logpdf(prior_positions, vcat(BL.q(x0)...)) + ######################### @@ -104,7 +106,7 @@ start = time() # to compute elapsed time Xsave, parsave, objvals, accpcn, accinfo, δ, ρ, covθprop = lm_mcmc(tt_, (xobs0,xobsT), Σobs, mT, P, obs_atzero, fixinitmomentato0, ITER, subsamples, - xinit, tp, priorθ, logpriormom, updatescheme, + xinit, tp, priorθ, priormom, updatescheme, outdir) elapsed = time() - start @@ -114,9 +116,9 @@ perc_acc_pcn = mean(accpcn)*100 println("Acceptance percentage pCN step: ", round(perc_acc_pcn;digits=2)) write_mcmc_iterates(Xsave, tt_, n, nshapes, subsamples, outdir) write_info(model,ITER, n, tt_, updatescheme, Σobs, tp, ρ, δ, perc_acc_pcn, elapsed, outdir) -write_acc(accinfo,accpcn,nshapes,outdir) +write_acc(accinfo, accpcn, nshapes, outdir) v0 = Vector(xobs0) vT = [Vector(xobsT[1])] write_observations(v0, vT, n, nshapes, x0,outdir) -write_noisefields(P,outdir) +write_noisefields(P, outdir) diff --git a/scripts/postprocessing.jl b/scripts/postprocessing.jl index 6392eef..45fcb6d 100644 --- a/scripts/postprocessing.jl +++ b/scripts/postprocessing.jl @@ -8,6 +8,7 @@ outdir: output directory """ function write_mcmc_iterates(Xsave, tt_, n, nshapes, subsamples, outdir) + outdir[end] == "/" && error("provide pathname without trailing '/'") nshapes = length(xobsT) iterates = reshape(vcat(Xsave...),2*d*length(tt_)*n*nshapes, length(subsamples)) # each column contains samplepath of an iteration @@ -34,7 +35,7 @@ function write_mcmc_iterates(Xsave, tt_, n, nshapes, subsamples, outdir) headline = "time " * "pqtype " * "landmarkid " * "shapes " * prod(map(x -> "iter"*string(x)*" ",subsamples)) headline = chop(headline,tail=1) * "\n" - fn = outdir*"iterates.csv" + fn = joinpath(outdir, "iterates.csv") f = open(fn,"w") write(f, headline) writedlm(f,out) @@ -46,7 +47,9 @@ end Write info to txt file """ function write_info(model,ITER, n, tt_, updatescheme, Σobs, tp, ρ, δ, perc_acc_pcn, elapsed, outdir) - f = open(outdir*"info.txt","w") + outdir[end] == "/" && error("provide pathname without trailing '/'") + + f = open(joinpath(outdir,"info.txt"),"w") write(f, "Model: ", string(model),"\n") write(f, "Number of iterations: ",string(ITER),"\n") write(f, "Number of landmarks: ",string(n),"\n") @@ -65,13 +68,15 @@ function write_info(model,ITER, n, tt_, updatescheme, Σobs, tp, ρ, δ, perc_ac end function write_acc(accinfo,accpcn,nshapes,outdir) + outdir[end] == "/" && error("provide pathname without trailing '/'") + # extract number of distinct update steps in accinfo nunique = length(unique(map(x->x[1], accinfo))) niterates = div(length(accinfo),nunique) accdf = DataFrame(kernel = map(x->Symbol(x.kernel), accinfo), acc = map(x->x.acc, accinfo), iter = repeat(1:niterates, inner= nunique)) accpcndf = DataFrame(kernel = fill(Symbol("pCN"),length(accpcn)), acc=accpcn, iter=repeat(1:niterates,inner=nshapes)) append!(accdf, accpcndf) - CSV.write(outdir*"accdf.csv", accdf; delim=";") + CSV.write(joinpath(outdir, "accdf.csv"), accdf; delim=";") end @@ -79,6 +84,8 @@ end Write observations to file """ function write_observations(xobs0, xobsT, n, nshapes, x0,outdir) + outdir[end] == "/" && error("provide pathname without trailing '/'") + valueT = vcat(map(x->deepvec(x), xobsT)...) # merge all observations at time T in one vector if d==1 posT = repeat(["pos1"], n*nshapes) @@ -100,8 +107,8 @@ function write_observations(xobs0, xobsT, n, nshapes, x0,outdir) obs0df = DataFrame(pos1=extractcomp(q0,1), pos2=extractcomp(q0,2), pos3=extractcomp(q0,3), mom1=extractcomp(p0,1) , mom2=extractcomp(p0,2), mom3=extractcomp(p0,3),landmark=1:n) end - CSV.write(outdir*"obs0.csv", obs0df; delim=";") - CSV.write(outdir*"obsT.csv", obsTdf; delim=";") + CSV.write(joinpath(outdir, "obs0.csv"), obs0df; delim=";") + CSV.write(joinpath(outdir, "obsT.csv"), obsTdf; delim=";") end @@ -109,15 +116,19 @@ end Write parameter iterates to file """ function write_params(parsave,subsamples,outdir) + outdir[end] == "/" && error("provide pathname without trailing '/'") + parsdf = DataFrame(a=extractcomp(parsave,1),c=extractcomp(parsave,2), gamma=extractcomp(parsave,3), iterate=subsamples) - CSV.write(outdir*"parameters.csv", parsdf; delim=";") + CSV.write(joinpath(outdir, "parameters.csv"), parsdf; delim=";") end """ Write noisefields to file (empty in case of MS-model) """ function write_noisefields(P,outdir) + outdir[end] == "/" && error("provide pathname without trailing '/'") + if isa(P,Landmarks) nfsloc = [P.nfs[j].δ for j in eachindex(P.nfs)] if d==2 @@ -131,5 +142,5 @@ function write_noisefields(P,outdir) elseif isa(P,MarslandShardlow) nfsdf =DataFrame(locx=Int64[], locy=Int64[], nfstd=Int64[]) end - CSV.write(outdir*"noisefields.csv", nfsdf; delim=";") + CSV.write(joinpath(outdir, "noisefields.csv"), nfsdf; delim=";") end diff --git a/src/BridgeLandmarks.jl b/src/BridgeLandmarks.jl index f99d419..f360aa2 100644 --- a/src/BridgeLandmarks.jl +++ b/src/BridgeLandmarks.jl @@ -16,6 +16,7 @@ using JLD using TimerOutputs const to = TimerOutput() +dir() = joinpath(@__DIR__, "..") const d = 2 const sk = 1 # entries to skip for likelihood evaluation From e503269c1369f11a4cb93a3fc5bedd1e08831c96 Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 13:34:49 +0100 Subject: [PATCH 14/20] Fix an_exp1.jl --- experiments/exp1/an_exp1.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/experiments/exp1/an_exp1.jl b/experiments/exp1/an_exp1.jl index 1a2b650..2eddd80 100644 --- a/experiments/exp1/an_exp1.jl +++ b/experiments/exp1/an_exp1.jl @@ -18,8 +18,11 @@ Random.seed!(9) workdir = @__DIR__ cd(workdir) -include(dirname(dirname(workdir))*"/postprocessing.jl") -outdir = workdir*("/") +include(joinpath(BL.dir(),"scripts", "postprocessing.jl")) +outdir = workdir +mkpath(joinpath(outdir, "forward")) + + #-------- read data ---------------------------------------------------------- dat = load("data_exp1.jld") @@ -99,7 +102,7 @@ priormom = MvNormalCanon( vcat(BL.p(x0)...), gramkernel(x0.q,P)/κ) ######################### xobsT = [xobsT] xinit = State(xobs0, zeros(PointF,P.n)) -mT = zeros(PointF,n) +mT = zeros(PointF, n) start = time() # to compute elapsed time Xsave, parsave, objvals, accpcn, accinfo, δ, ρ, covθprop = @@ -116,8 +119,8 @@ println("Acceptance percentage pCN step: ", round(perc_acc_pcn;digits=2)) write_mcmc_iterates(Xsave, tt_, n, nshapes, subsamples, outdir) write_info(model,ITER, n, tt_, updatescheme, Σobs, tp, ρ, δ, perc_acc_pcn, elapsed, outdir) write_observations(xobs0, xobsT, n, nshapes, x0,outdir) -write_acc(accinfo,accpcn,nshapes,outdir) -write_params(parsave,0:ITER,outdir) -write_noisefields(P,outdir) +write_acc(accinfo, accpcn, nshapes,outdir) +write_params(parsave, 0:ITER, outdir) +write_noisefields(P, outdir) #show(to; compact = true, allocations = true, linechars = :ascii) From cd425bc916e07b755e1fd40157cbf245913d6ccf Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 14:04:45 +0100 Subject: [PATCH 15/20] use JLD2 --- experiments/exp1/an_exp1.jl | 4 ++-- experiments/exp1/gen_exp1.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/experiments/exp1/an_exp1.jl b/experiments/exp1/an_exp1.jl index 2eddd80..86dd59d 100644 --- a/experiments/exp1/an_exp1.jl +++ b/experiments/exp1/an_exp1.jl @@ -12,7 +12,7 @@ using DelimitedFiles using CSV using StaticArrays using LinearAlgebra -using JLD +using JLD2 Random.seed!(9) @@ -25,7 +25,7 @@ mkpath(joinpath(outdir, "forward")) #-------- read data ---------------------------------------------------------- -dat = load("data_exp1.jld") +dat = load("data_exp1.jld2") xobs0 = dat["xobs0"] xobsT = dat["xobsT"] n = dat["n"] diff --git a/experiments/exp1/gen_exp1.jl b/experiments/exp1/gen_exp1.jl index 7d2e2ea..6a09f9c 100644 --- a/experiments/exp1/gen_exp1.jl +++ b/experiments/exp1/gen_exp1.jl @@ -6,7 +6,7 @@ using Random using Distributions using StaticArrays using LinearAlgebra -using JLD +using JLD2 workdir = @__DIR__ println(workdir) @@ -37,4 +37,4 @@ stretch = SMatrix{2,2}(1.0 + ψ, 0.0, 0.0, 1.0 - ψ) shift = [0.5, 0.0] xobsT = [rot * stretch * xobs0[i] + shift for i in 1:n ] -save("data_exp1.jld", "xobs0",xobs0, "xobsT", xobsT, "n", n, "x0", x0, "nshapes", nshapes) +JLD2.@save "data_exp1.jld2" xobs0 xobsT n x0 nshapes From a7db6b7c71cec72b579b338ab6da855c775c1f89 Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 14:36:33 +0100 Subject: [PATCH 16/20] Fix deprecations --- Project.toml | 2 +- experiments/exp1/an_exp1.jl | 1 + src/lmguid.jl | 5 +++-- src/nstate.jl | 3 +++ 4 files changed, 8 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 328be5c..e9d9085 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] -Bridge = "0.10" +Bridge = "0.10, 0.11" DataFrames = "0.20" Distributions = "0.22" ForwardDiff = "0.10" diff --git a/experiments/exp1/an_exp1.jl b/experiments/exp1/an_exp1.jl index 86dd59d..67a75ef 100644 --- a/experiments/exp1/an_exp1.jl +++ b/experiments/exp1/an_exp1.jl @@ -13,6 +13,7 @@ using CSV using StaticArrays using LinearAlgebra using JLD2 +using FileIO Random.seed!(9) diff --git a/src/lmguid.jl b/src/lmguid.jl index 0c4db38..41f0927 100644 --- a/src/lmguid.jl +++ b/src/lmguid.jl @@ -224,12 +224,13 @@ function guidingbackwards!(::Lm, t, (Lt, Mt⁺, μt), Paux, obs_info; implicit=t oldtemp = (0.5*dt) * Bridge.outer(Lt[end] * σ̃T) if lowrank # TBA lowrank on σ̃T, and write into σ̃T + error("not implemented") end for i in length(t)-1:-1:1 dt = t[i+1]-t[i] if implicit - Lt[i] .= Lt[i+1]/(I - dt* B̃) + Lt[i] .= Lt[i+1]/lu(I - dt* B̃, Val(false)) else Lt[i] .= Lt[i+1] * (I + B̃ * dt) end @@ -455,7 +456,7 @@ function update_initialstate!(X,Xᵒ,W,ll,x,xᵒ,∇x, ∇xᵒ, if sampler ==:sgd # CHECK VALIDITY LATER - mask = deepvec(State(0*x0.q, 1 .- 0*x0.p)) + mask = deepvec(State(0*x0.q, onemask(x0.p))) # StateW = PointF # sample!(W, Wiener{Vector{StateW}}()) ForwardDiff.gradient!(∇x, u, x, cfg) # X gets overwritten but does not chang diff --git a/src/nstate.jl b/src/nstate.jl index 3c02f0e..fffd5b6 100644 --- a/src/nstate.jl +++ b/src/nstate.jl @@ -117,3 +117,6 @@ q(i::Int) = 2i - 1 p(i::Int) = 2i flipmomenta(x::NState) = NState(x.q, -x.p) + +onemask(x) = onemask.(x) +onemask(x::Number) = one(x) From c70309ca449575f77732cb444bc3c670eb6ef80b Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 15:47:40 +0100 Subject: [PATCH 17/20] Pivoting? --- src/lmguid.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lmguid.jl b/src/lmguid.jl index 41f0927..9d5e303 100644 --- a/src/lmguid.jl +++ b/src/lmguid.jl @@ -230,7 +230,7 @@ function guidingbackwards!(::Lm, t, (Lt, Mt⁺, μt), Paux, obs_info; implicit=t for i in length(t)-1:-1:1 dt = t[i+1]-t[i] if implicit - Lt[i] .= Lt[i+1]/lu(I - dt* B̃, Val(false)) + Lt[i] .= Lt[i+1]/lu(I - dt* B̃, Val(false)) # should we use pivoting? else Lt[i] .= Lt[i+1] * (I + B̃ * dt) end From b3c82ffff6c7c5e42b94eb41348900a930ebd3ed Mon Sep 17 00:00:00 2001 From: mschauer Date: Wed, 29 Jan 2020 16:03:44 +0100 Subject: [PATCH 18/20] Fixes to experiments --- experiments/exp1-1D/an_exp1-1D.jl | 12 ++++++++---- experiments/exp1-1D/gen_exp1-1D.jl | 4 ++-- experiments/exp1/data_exp1.jld2 | Bin 0 -> 11617 bytes experiments/exp2/an_exp2.jl | 9 ++++++--- experiments/exp2/data_exp2.jld2 | Bin 0 -> 10234 bytes experiments/exp2/gen_exp2.jl | 4 ++-- 6 files changed, 18 insertions(+), 11 deletions(-) create mode 100644 experiments/exp1/data_exp1.jld2 create mode 100644 experiments/exp2/data_exp2.jld2 diff --git a/experiments/exp1-1D/an_exp1-1D.jl b/experiments/exp1-1D/an_exp1-1D.jl index 577c529..140544e 100644 --- a/experiments/exp1-1D/an_exp1-1D.jl +++ b/experiments/exp1-1D/an_exp1-1D.jl @@ -12,17 +12,21 @@ using DelimitedFiles using CSV using StaticArrays using LinearAlgebra -using JLD +using JLD2 +using FileIO Random.seed!(9) + workdir = @__DIR__ cd(workdir) -include(dirname(dirname(workdir))*"/postprocessing.jl") -outdir = workdir*("/") +include(joinpath(BL.dir(),"scripts", "postprocessing.jl")) +outdir = workdir +mkpath(joinpath(outdir, "forward")) + #-------- read data ---------------------------------------------------------- -dat = load("data_exp1-1D.jld") +dat = load("data_exp1-1D.jld2") xobs0 = dat["xobs0"] xobsT = dat["xobsT"] n = dat["n"] diff --git a/experiments/exp1-1D/gen_exp1-1D.jl b/experiments/exp1-1D/gen_exp1-1D.jl index f44b57a..616d48e 100644 --- a/experiments/exp1-1D/gen_exp1-1D.jl +++ b/experiments/exp1-1D/gen_exp1-1D.jl @@ -6,7 +6,7 @@ using Random using Distributions using StaticArrays using LinearAlgebra -using JLD +using JLD2 workdir = @__DIR__ println(workdir) @@ -31,4 +31,4 @@ xobs0 = x0.q xobsT = [exp.(xobs0[i]) for i in 1:n ] #xobsT = [PointF(1.0), PointF(0.4)] -save("data_exp1-1D.jld", "xobs0",xobs0, "xobsT", xobsT, "n", n, "x0", x0, "nshapes", nshapes) +JLD2.@save "data_exp1-1D.jld2" xobs0 xobsT n x0 nshapes diff --git a/experiments/exp1/data_exp1.jld2 b/experiments/exp1/data_exp1.jld2 new file mode 100644 index 0000000000000000000000000000000000000000..78ca364fad9d59fe1bb0d1d8c1f1fbae71d46e48 GIT binary patch literal 11617 zcmeHN3sh4_8lEH&E+7>YLHjbIC`wHO1qxQ3i1@-Ntx{L3)dmch784-BcdPYTt*+GP z7P{4H7oXjQ#qRpxI9gv+ighdNqE?I81B*NrP#y{F%$>Osb8QRfboZP+&72c%?*Gl> z|Nr^tpP7$DrzPpNin&^|R*|4j(kTMQM~w{~{DNYEF4d$r8WbvJuu=sS0nvQDAZ4gh zH8g6>G-Xg&up%sEXq?`x(2N-j9OeJi1MBg8d$>LACgnH)q|}IZR0*>Y$Q;J;pz&E7_{?s07Q66fG2)a;g1A=0N_BfHWg!L97@150!s+~GLHbz zZ4@3);tt$^6L^~h*BtZWWF3%6Ne845WjydFpS<8gJEh-+zhN1!8v1?Q=ok%OgN%WT zgw-cKkTTkss#8YcFpnjpklty3Pw^-Kky?{ZIn}I7o}x7*;tDc`yD^+^Kfz!Q3!wrh zhO^SJe{A^J~k~mi42a5oZ;B(d?{ufxiGn)H`OaGZd7Wjb}_EW{X)O44DZJJ z<@Fw$WYjXv-54(D9qUai3Qu>2_hjv-nzd$qJU_ZgIaSaCPv0--;lXghzOQ=G4xeIj z1?%6$^;PXR_OkuaXQq%Aw}m=uxBr`r;)YyVz?_CkM~EU(U`>EdoVnJy$kK< zv-Vl+d7RDf247mSCu_%G5SglF|KV0}c;3JmBzv!vhWvI6UC+fWre05Bx?vKzWJ}WDxTd zG5~sGvrT4gsu>d$_)1D|rD`nFB?Inpi7^>?+DbY_R1g{=pShM07vv+nJK^uDl`;RM zDn&VCOy(-VQc_9g24o$a2R8s$tTQcTd$n>zpk_o05XJes#k zqLj$CW3#Z2RJcD(Bk2ry<^<|UZjLfFd4_eJmZcc<~NToy)|IYKDlcnDI!N6{YQ z4no>-I0XalN#agfW;^#B9>O-z38qaiTG)j^Cp4Ad7=fa}qDN)Oh#7rpcjshLG|-7* zB%W}r%uAwh%syNW_wkGWKeC5>?9IZXWsVLCHu*4(AaPI3$`98SbeIED5hn*%O_DZ-p+sTc@^zaanPWuM)A_Y#eq^^fC+r0h7>f^J&^ zS;`HN<8c)7>f%@QC-d6j=1qsR<7S_Nta8VrI(nveLmR9q*l8?z<+2FpG&f@>9k9Zi zag)pDI2DUf8B0{d)Cg}c{lL~5uTgZU8($2Gu%6|HRovmJ0iUIv7_hzZg~LD z98UD$Mwf~3@2Wp6-CbV~C*E>SJL6O#!Y@QtNe}O+gDsh!Gh^SZ6ydSs0(DL06)^Pn z_mQ_Et3~*jhs8m;p`~zJW_rlHFKR@1!KzCQXXh8e^H0xFx9zAE;dQ|kOFZ@2@M28T z>Wi&nTzjI{{o%HJq}V-3<9}U@uNWh}lfLwOq_3Fn{7JAFZ>wJ7vf$}5ROtFyb&ZQy zzxg)|FAdA6M(yRQmCkdjMEZ^N@%{M1D1X9 zBK-UbW5B`LZRpgC7xteRQY6Cv{bd2Fnb(f)7-I@c7v_ubq;LBT4XkZPHyamJmM3KY zBz~CUs3Y*{w+hyO0qtLe=QI6|Q2j)B3bPMU`$YH-<{y9R9}#|+`8SvPSA?tB{3z-C zi15*DemBth72zo?KH6w}h;X+W8eau8zC`%?$Kz9k?_u%Xm&UgUKm7Rm5a9>e`pT#4 zON2+T^|^ztPZ1u;)^`M5-y-}p%MVtXA4K?NmR~e9zliXjEI;Mb{3OD&S$@l+`Avim zXZbOk=0_3k!t!ew&95TdT1E3S_=eBVBD{*__bi&$f#<^+cH4ud)&5pKI7QN?Q&m zt7ha!hStN`hpt@&`+k6N8gle@edu+h>DNKKlgHBBWf9AIvcmT5q*B` z_?3BqP3V(L`OfdgHKIwwq<@hNXh9!*z9XT>xdzmDdvC_kKewVL{#lRK@Nx@!jjwkf zW$^8+=m6hO9qZSMj-*@|X2_^S1^l>Hw5%*D?)1Pil*sF4L3@<)8(;P)LBn~yTTyJG z2gs@Y4(DZJS2IfA^x3nCn{$!j^z_QW0T0oFK_%OtTYCmJZS7UaWj3I|%&?;=7mA^* zpjk5h#62W|AFlUpyaTU4U%bWAuNr+bBjwPQv`Sc5x;(yqN*TJP4)|M8{5|;8xeT9~)rt`|}!o$}XXi{5&>6@Z$Cfg>V5H$Io*syy|mhs(ei$T+YXV72f-O z^2(N7H{o~sJw1(nKpitJRhx2`1;rpemqI6aQww;N%2oc1V0J`*rhASF_?Ft;~4BX29sm5;~30elSIcP z$1#{*vwy3{;7gbLw0CJg6L^~r@~uwFQ#=B3BRt*sjY*DoD~YFsQ_&KnV%wUXJXqk} z?c^^BEVG?FM4*-It&X(ta8=pj4#%L1KcCl43)W1iGyItvf!$NMeRX9s`QbT8g? z=S9Xilj=C5&Ez)1*p%mKmvYRXVkb!q$Ji-Aj(J@ykS*#=8o45oQn{J)Du b>}Sq$?0PJg^FLTM*C|;x`M7AB z;M?4qScS08cao3s7<~QcTHxd|3NQf#6^1t@n8EB8=iih|)`tU2?wJE-#kf zLR~^6Oz9~%F-IA4Aiox&^b@H0>-|`M?%S3T4)|)W=hZw7JY_c!`*N>441VNxGkV~sHc#t;}oU<`pV1jY~;LtqSnF$BgC7(-wT zf&T{xkUWKn5Y8e9h|&Z}Izd1T6Zn&0U?^ia_n?UL6n~No6a(@ZK$2nA=SXuKD06sn zPv;X1_)rO)C!u5FxS`Le3_jCA*es^mm=Z?i(yPmR~*wF!KW|47=I=HN9@d}o)Ce2>WMJk z&(Y7G^%a2;lnFl&A-w%%aSL+#jWm}jV1fsR6F5X*1b2yH1R3?e6ehut@ck0O+eE1R zJ^_IQPVCU}#Tl3cw%y%jm$L)m62WHqeq_6d?SlKuD~039nrS{%`4#YdiZ%ncK$zpZ zA&R$;S1BZt<}|<9^*BePh0h=+5s#jM11*U6h{r{gPzL+$S54#x5+S)z?W-pK53YFa z^P@mU5Zaf3+{2<*8$Ta`A7!4p-8ozc zYm^O~SW`a0n(1?&g?TiB`K8zC=LNeUdtkTI)ps3Ge9SGfDWV7_$y}Lq%%%(WbmuDn z?E4l;Fk04kP3;7+r8aI-v$~L&Zz$8wtpQ>Wt8JWWG>AB@Vk9L!jC?FY z(hDpvqigfzKYR;L)6c96>oFTa=TBKT--&#T>W`fX6cKulj+F=wZd_e~%)Pn)PW7!t zH(b|WJfqZz45!-D#_MOpW;(~)a7Hs?$B8>B&FBElkVe;Evl>uSW#|)O>j4OCPu;4u zs2s@%?iES!;DY$HiF-FlB%+k4S+lsJTzI0rKIC3SKDY)p%5=W$2j!)icP2lq1-%+Q zDOqMK1fHP#PH}04L?P*`S~*E@fT5LCFuoIJX|ES>ZfHao@2j2Fw!!3`%%Ro z9!L8_+hN{Knf6FEjG1}u<+hx|^>83oFWAbF#tgS)S3lU55BtximVFCOE7H%66wIMB z0~gC>%@e9dSLiOCY10*Lm7S{aS4@QCsFUqG>fo{Ua6PwR-;NqO)dUB#1u=-_@j!Dr~ zWF_l(T;-=$c(^J)vEQN{iPoRdDK%?`ebN)-^RuYF*``Ye52r$De-<(QUPe=a6DM!8qrx%OoPb(PWB6>azUWFB)E5}(G<`#X(RQ6)5 zvrRqHpD(4OtJ#2y>`S5YAQVZlmjQF?M;7Q(1LnX zuSvc=& zm3|9O<;9tPL9GMmq~!1KO6;pp+-|vb4O2VOJJXjhHa==aa_TYfv)HfEzfvpRzmHstghxTm>%rWP5@oY!?@ zl>mE|n&IW&grA|FdF8XMcL=awr8EiARbIfI)ID2H-NM**b&o4^j#k5OChxC>=e&pQ zWj2@MpEN=&V@iCOY8{-h*}kUIm;+*6cl##HeU3aAO_%9oH-OHC;_^E6 zY8Vb*GO$~C5Cz!Q2b?>f2j=#R)tfI4p)}Pa`n#P|QK;cX*Kfh;*5bT8(*qdOjaSnYQVJARnsoa)AoJoZxk7Fm$o!!`tvD8An!bOg`u(#yxO{Tf z@kul~^Tl4rSL~wK@NA21+qdA9VWPniF&;$oSQb%Y8AYfy=VX3CNhiXpHCA2;tw%dq z5%LfH3XuG7O9bL&nh`xIedg1-l~8l`)M@N^Gm3WTr)yU9LbbMnr&W3#dc23@xy5Y+ zbUQ*8PqZmPPH*oydd;V^Z=X@NXjXAYCKpcL&}Y!ulULN7wSQIyQX`G(54ppjxB8TS zmNN(LUNx|~`Met<#f;S0Ed1}=rI^V&Ka@h%lB{Pt_xHe-+{@ovoh?IwX(+WX#Q!gFeGKi}Dd_I(7m07B;dG)+YANM*y?BZUf(lV0bxfP#SfVbLn+;{bn%8}Xg`EzY^&`84cb+#=WQo& z7cI)3b2A@kiVBhv!9B2khwK4|o+?x!c7O8O?|NXmkdWoUcReUMPi)+yOPz2k`}{qx zH$#Xq;Q7Kmv+c(t49Kz4X-goSnN)(#eBRlgzCt`Q) zI$FJlgH-gI!jv|>LC1_H*`F0{M=8ATuWhsagkwW zSuL=SO0hJ;8US0waeDAF6&(zCm$D(c4c3ON-K)oJMc0k|e?EM!1KiG+hkU0#h`Q}+ zc4Md&@|5FLx0ukFeh;ziVVPREA({2dRYy8gQ)zK}vUDCayVZ^7;e`Fe9Etxm{l@ct z@(U_vFED}3J&AP6|10k^5NQQJXCu;GA|c}UiKPN7{&dr?}dPHYkWb^T8b*;?B)U%&lH z4H~gFnIPrWHcIbgV-)R}o();9{SheQTiJ_KS%GY82|^n zJ%Y1xt4HaBuz&?amLnjzY5DDvbJ?Tx!IZwmmZyg={g0*WZz1iE(i6-BhIPQNq<^FI zSE%@&knxSuTT=L78+iDP(&tn7W|H`h(uY&er-FPwqx4ifOcDEuIx5!_9#Q>|Rj#>h z&{yq8ep(+CEVY1~_ye7H0fm1>iC>HHHkL5H+i21>h89`0wWyW9kT2jfFki@Z_-XhT za$Wv@=P^tlqahoI9L)}quHF1EUBg)i+%+M8FDGBEITQsU{E5&y5qF*_!uUzXN76jI z6nTA27(dFmfbOp~jXGz-I3 Date: Wed, 29 Jan 2020 16:03:53 +0100 Subject: [PATCH 19/20] Fix deprecations --- src/lmguid.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/lmguid.jl b/src/lmguid.jl index 9d5e303..6beb657 100644 --- a/src/lmguid.jl +++ b/src/lmguid.jl @@ -476,10 +476,10 @@ function update_initialstate!(X,Xᵒ,W,ll,x,xᵒ,∇x, ∇xᵒ, ForwardDiff.gradient!(∇x, u, x, cfg) # X gets overwritten but does not change ll_incl0 = sum(llout) if update==:mala_pos - mask = deepvec(State(1 .- 0*x0.q, 0*x0.p)) + mask = deepvec(State(onemask(x0.q), 0*x0.p)) stepsize = δ[1] elseif update==:mala_mom - mask = deepvec(State(0*x0.q, 1 .- 0*x0.p)) + mask = deepvec(State(0*x0.q, onemask(x0.p))) stepsize = δ[2] end mask_id = (mask .> 0.1) # get indices that correspond to positions or momenta @@ -494,7 +494,7 @@ function update_initialstate!(X,Xᵒ,W,ll,x,xᵒ,∇x, ∇xᵒ, elseif update == :rmmala_pos ForwardDiff.gradient!(∇x, u, x, cfg) # X gets overwritten but does not change ll_incl0 = sum(llout) - mask = deepvec(State(1 .- 0*x0.q, 0*x0.p)) + mask = deepvec(State(onemask(x0.q), 0*x0.p)) stepsize = δ[1] mask_id = (mask .> 0.1) # get indices that correspond to positions or momenta dK = gramkernel(x0.q,P) @@ -515,7 +515,7 @@ function update_initialstate!(X,Xᵒ,W,ll,x,xᵒ,∇x, ∇xᵒ, elseif update == :rmmala_mom ForwardDiff.gradient!(∇x, u, x, cfg) # X gets overwritten but does not change ll_incl0 = sum(llout) - mask = deepvec(State(0*x0.q, 1 .- 0*x0.p)) + mask = deepvec(State(0*x0.q, onemask(x0.p))) stepsize = δ[2] mask_id = (mask .> 0.1) # get indices that correspond to positions or momenta @@ -538,7 +538,7 @@ function update_initialstate!(X,Xᵒ,W,ll,x,xᵒ,∇x, ∇xᵒ, elseif update == :rmrw_mom u(x) # writes into llout ll_incl0 = sum(llout) - mask = deepvec(State(0*x0.q, 1 .- 0*x0.p)) + mask = deepvec(State(0*x0.q, onemask(x0.p))) stepsize = δ[2] mask_id = (mask .> 0.1) # get indices that correspond to positions or momenta @@ -846,10 +846,10 @@ if false ForwardDiff.gradient!(∇x, u, x,cfg) # X gets overwritten but does not change ll_incl0 = sum(llout) if update==:mala_pos - mask = deepvec(State(1 .- 0*x0.q, 0*x0.p)) + mask = deepvec(State(onemask(x0.q), 0*x0.p)) stepsize = δ[1] elseif update==:mala_mom - mask = deepvec(State(0*x0.q, 1 .- 0*x0.p)) + mask = deepvec(State(0*x0.q, onemask(x0.p))) stepsize = δ[2] end mask_id = (mask .> 0.1) # get indices that correspond to positions or momenta @@ -868,7 +868,7 @@ if false cfg = ForwardDiff.GradientConfig(slogρ!(Q, W, X,priormom,llout), x, ForwardDiff.Chunk{2*d*n}()) # 2*d*P.n is maximal ForwardDiff.gradient!(∇x, slogρ!(Q, W, X,priormom,llout),x,cfg) # X gets overwritten but does not change ll_incl0 = sum(llout) - mask = deepvec(State(1 .- 0*x0.q, 0*x0.p)) + mask = deepvec(State(onemask(x0.q), 0*x0.p)) stepsize = δ[1] mask_id = (mask .> 0.1) # get indices that correspond to positions or momenta dK = gramkernel(x0.q,P) #reshape([kernel(x0.q[i]- x0.q[j],P) * one(UncF) for i in 1:n for j in 1:n], n, n) @@ -893,7 +893,7 @@ if false cfg = ForwardDiff.GradientConfig(slogρ!(Q, W, X,priormom,llout), x, ForwardDiff.Chunk{2*d*n}()) # 2*d*P.n is maximal ForwardDiff.gradient!(∇x, slogρ!(Q, W, X,priormom,llout),x,cfg) # X gets overwritten but does not change ll_incl0 = sum(llout) - mask = deepvec(State(0*x0.q, 1 .- 0*x0.p)) + mask = deepvec(State(0*x0.q, onemask(x0.p))) stepsize = δ[2] mask_id = (mask .> 0.1) # get indices that correspond to positions or momenta From 4b262ea18aafbfd61d31b8f8a69279c24fe6025e Mon Sep 17 00:00:00 2001 From: mschauer Date: Thu, 30 Jan 2020 12:08:33 +0100 Subject: [PATCH 20/20] Move inprogress stuff --- .../{ => inprogress}/exp1-try/an_exp1-try.jl | 0 .../{ => inprogress}/exp1-try/data_exp1-try.jld | Bin .../{ => inprogress}/exp1-try/forwardhalfway.jld | Bin .../{ => inprogress}/exp1-try/gen_exp1-try.jl | 0 .../{ => inprogress}/exp1-try/makefigs_exp1-try.R | 0 experiments/{ => inprogress}/prior/prior_initpos.jl | 0 6 files changed, 0 insertions(+), 0 deletions(-) rename experiments/{ => inprogress}/exp1-try/an_exp1-try.jl (100%) rename experiments/{ => inprogress}/exp1-try/data_exp1-try.jld (100%) rename experiments/{ => inprogress}/exp1-try/forwardhalfway.jld (100%) rename experiments/{ => inprogress}/exp1-try/gen_exp1-try.jl (100%) rename experiments/{ => inprogress}/exp1-try/makefigs_exp1-try.R (100%) rename experiments/{ => inprogress}/prior/prior_initpos.jl (100%) diff --git a/experiments/exp1-try/an_exp1-try.jl b/experiments/inprogress/exp1-try/an_exp1-try.jl similarity index 100% rename from experiments/exp1-try/an_exp1-try.jl rename to experiments/inprogress/exp1-try/an_exp1-try.jl diff --git a/experiments/exp1-try/data_exp1-try.jld b/experiments/inprogress/exp1-try/data_exp1-try.jld similarity index 100% rename from experiments/exp1-try/data_exp1-try.jld rename to experiments/inprogress/exp1-try/data_exp1-try.jld diff --git a/experiments/exp1-try/forwardhalfway.jld b/experiments/inprogress/exp1-try/forwardhalfway.jld similarity index 100% rename from experiments/exp1-try/forwardhalfway.jld rename to experiments/inprogress/exp1-try/forwardhalfway.jld diff --git a/experiments/exp1-try/gen_exp1-try.jl b/experiments/inprogress/exp1-try/gen_exp1-try.jl similarity index 100% rename from experiments/exp1-try/gen_exp1-try.jl rename to experiments/inprogress/exp1-try/gen_exp1-try.jl diff --git a/experiments/exp1-try/makefigs_exp1-try.R b/experiments/inprogress/exp1-try/makefigs_exp1-try.R similarity index 100% rename from experiments/exp1-try/makefigs_exp1-try.R rename to experiments/inprogress/exp1-try/makefigs_exp1-try.R diff --git a/experiments/prior/prior_initpos.jl b/experiments/inprogress/prior/prior_initpos.jl similarity index 100% rename from experiments/prior/prior_initpos.jl rename to experiments/inprogress/prior/prior_initpos.jl