Skip to content

Commit

Permalink
Merge pull request #15 from mschauer/refactor2
Browse files Browse the repository at this point in the history
Refactor
  • Loading branch information
mschauer authored Jan 30, 2020
2 parents 25e9533 + 4b262ea commit bad06bb
Show file tree
Hide file tree
Showing 84 changed files with 248 additions and 1,204 deletions.
8 changes: 4 additions & 4 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -21,17 +21,17 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"

[compat]
Bridge = "0.10"
DataFrames = "0.20"
Bridge = "0.10, 0.11"
DataFrames = "0.20"
Distributions = "0.22"
ForwardDiff = "0.10"
GaussianDistributions = "0.3"
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]
Expand Down
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
18 changes: 10 additions & 8 deletions src/experiments/exp0/exp0.jl → experiments/exp0/exp0.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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 ---------------------------------------------------------

Expand Down Expand Up @@ -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)...))


#########################

Expand All @@ -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

Expand All @@ -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)
File renamed without changes.
File renamed without changes.
Original file line number Diff line number Diff line change
Expand Up @@ -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"]
Expand Down
File renamed without changes.
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ using Random
using Distributions
using StaticArrays
using LinearAlgebra
using JLD
using JLD2

workdir = @__DIR__
println(workdir)
Expand All @@ -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
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
20 changes: 12 additions & 8 deletions src/experiments/exp1/an_exp1.jl → experiments/exp1/an_exp1.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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.jld")
dat = load("data_exp1.jld2")
xobs0 = dat["xobs0"]
xobsT = dat["xobsT"]
n = dat["n"]
Expand Down Expand Up @@ -99,7 +103,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 =
Expand All @@ -116,8 +120,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)
File renamed without changes.
Binary file added experiments/exp1/data_exp1.jld2
Binary file not shown.
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ using Random
using Distributions
using StaticArrays
using LinearAlgebra
using JLD
using JLD2

workdir = @__DIR__
println(workdir)
Expand Down Expand Up @@ -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
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
Original file line number Diff line number Diff line change
Expand Up @@ -10,14 +10,17 @@ 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_exp2.jld")
Expand Down
File renamed without changes.
File renamed without changes.
File renamed without changes.
Binary file added experiments/exp2/data_exp2.jld2
Binary file not shown.
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ using Random
using Distributions
using StaticArrays
using LinearAlgebra
using JLD
using JLD2

workdir = @__DIR__
println(workdir)
Expand Down Expand Up @@ -41,4 +41,4 @@ end
xobs0 = []


save("data_exp2.jld", "xobs0",xobs0, "xobsT", xobsT, "n", n, "x0", x0, "nshapes", nshapes)
JLD2.@save "data_exp2.jld2" xobs0 xobsT n x0 nshapes
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
12 changes: 6 additions & 6 deletions scripts/generatedata.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
File renamed without changes.
Loading

0 comments on commit bad06bb

Please sign in to comment.