BridgeLandmarks.jl

Construct sequence of Noisefields for AHS model
db: domainbound (sources are places on square grid specified by
    (-db:2nfstd:db) x -db:2nfstd:db
nfstd: standard deviation of noise fields (the smaller: the more noise fields we use)
γ: if set to one, then the value of the  noise field on the positions is approximately 1 at all locations in the domain

Forward simulate landmarks process specified by P on grid t.
Returns driving motion W and landmarks process X
t: time grid
x0: starting point
P: landmarks specification

Perform mcmc or sgd for landmarks model using the LM-parametrisation
tt_:      time grid
(xobs0,xobsT): observations at times 0 and T (at time T this is a vector)
σobs: standard deviation of Gaussian noise assumed on xobs0 and xobsT
mT: vector of momenta at time T used for constructing guiding term
P: target process

sampler: either sgd (stochastic gradient descent) or mcmc (Markov Chain Monte Carlo)
obs_atzero: Boolean, if true there is an observation at time zero
fixinitmomentato0: Boolean, if true we assume at time zero we observe zero momenta
xinit: initial guess on starting state

ITER: number of iterations
subsamples: vector of indices of iterations that are to be saved

ρ: Crank-Nicolson parameter (ρ=0 is independence sampler)
δ: step size for MALA updates on initial state (first component for positions, second component for momenta)
prior_a: prior on parameter a
prior_c: prior on parameter c
prior_γ: prior on parameter γ
σ_a: parameter determining update proposal for a [update a to aᵒ as aᵒ = a * exp(σ_a * rnorm())]
σ_c: parameter determining update proposal for c [update c to cᵒ as cᵒ = c * exp(σ_c * rnorm())]
σ_γ: parameter determining update proposal for γ [update γ to γᵒ as γᵒ = γ * exp(σ_γ * rnorm())]

outdir: output directory for animation
pb:: Lmplotbounds (axis used for plotting landmarks evolution)
updatepars: logical flag for updating pars a, c, γ
makefig: logical flag for making figures
showmomenta: logical flag if momenta are also drawn in figures

Returns:
Xsave: saved iterations of all states at all times in tt_
parsave: saved iterations of all parameter updates ,
objvals: saved values of stochastic approximation to loglikelihood
perc_acc: acceptance percentages (bridgepath - inital state)


anim, Xsave, parsave, objvals, perc_acc = lm_mcmc_mv(tt_, (xobs0,xobsTvec), σobs, mT, P,
         sampler, obs_atzero,
         xinit, ITER, subsamples,
        (δ, prior_a, prior_c, prior_γ, σ_a, σ_c, σ_γ),
        outdir, pb; updatepars = true, makefig=true, showmomenta=false)

GuidRecursions defines a struct that contains all info required for computing the guiding term and
likelihood (including ptilde term) for a single shape

struct that contains target, auxiliary process for each shape, time grid, observation at time 0, observations
    at time T, number of shapes, and momenta in final state used for constructing the auxiliary processes
guidrec is a vector of GuidRecursions, which contains the results from the backward recursions and gpupdate step at time zero

Observe V0 = L0 X0 + N(0,Σ0) and VT = LT X0 + N(0,ΣT)
μT is just a vector of zeros (possibly remove later)

Compute y = HX where Hinv = LL' (Cholesky decomposition

Input are L and X, output is Y

y=HX is equivalent to LL'y=X, which can be solved by first backsolving LZ=X for z and next backsolving L'Y=Z

L is a lower triangular matrix with element of type UncMat X is a matrix with elements of type UncMat Returns a matrix with elements of type UncMat

Compute y = Hx where Hinv = LL' (Cholesky decomposition

Input are L and x, output is y

y=Hx is equivalent to LL'y=x, which can be solved by first backsolving Lz=x for z and next backsolving L'y=z

L is a lower triangular matrix with element of type UncMat x is a State or vector of points Returns a State (Todo: split up again and return vector for vector input)

Compute transpose of square matrix of Unc matrices

A = reshape([Unc(1:4), Unc(5:8), Unc(9:12), Unc(13:16)],2,2) B = copy(A) A conj!(B)

Good display of variable of type State

Compute B̃(t) * X (B̃ from auxiliary process) and write to out Both B̃(t) and X are of type UncMat

Compute B̃(t) * X (B̃ from auxiliary process) and write to out Both B̃(t) and X are of type UncMat

Compute tildeB(t) for landmarks auxiliary process

Compute tildeB(t) for landmarks auxiliary process

Evaluate drift bᵒ of guided proposal at (t,x), write into out

Returns matrix a(t) for Marsland-Shardlow model

Evaluate drift of landmarks auxiliary process in (t,x) and save to out x is a state and out as well

Evaluate drift of landmarks in (t,x) and save to out x is a state and out as well

Evaluate drift of landmarks in (t,x) and save to out x is a state and out as well

Evaluate σ(t,x) dW and write into out

Compute sigma(t,x) * dm where dm is a vector and sigma is the diffusion coefficient of landmarks write to out which is of type State

compute σ(t,x) * dm and write to out

Adapting Radford Neal's R implementation of Hamiltonian Monte Carlo with
stepsize ϵ and L steps

Evaluate tilder (appearing in guiding term of guided proposal) at (t,x), write into out

Multiply a(t,x) times a vector of points Returns a State (first multiply with sigma', via function σtmul, next left-multiply this vector with σ)

Multiply a(t,x) times a state Returns a state (first multiply with sigma', via function σtmul, next left-multiply this vector with σ)

Multiply a(t,x) times xin (which is of type state) Returns variable of type State

Solve L L'y =x using two backsolves, L should be lower triangular

Useful for storage of a samplepath of states
Ordering is as follows:
1) time
2) landmark nr
3) for each landmark: q1, q2 p1, p2

With m time-steps, n landmarks, this entails a vector of length m * n * 2 * d

Useful for storage of a samplepath of states
Ordering is as follows:
0) shape
1) time
2) landmark nr
3) for each landmark: q1, q2 p1, p2

With m time-steps, n landmarks, this entails a vector of length m * n * 2 * d

Extract parameters from GuidedProposal! Q

Solve backwards recursions in L, M, μ parametrisation on grid t
implicit: if true, Euler backwards is used for solving ODE for Lt, else Euler forwards is used

Case lowrank=true still gives an error: fixme!

Hamiltonian for deterministic part of landmarks model

Initialise (allocate memory) a struct of type GuidRecursions for a single shape
guidres = init_guidrec((t,obs_info,xobs0)

kernel in Hamiltonian

Guided proposal update for newly incoming observation at time zero.
Information on new observations at time zero is (L0, Σ0, xobs0)
Values just after time zero, (Lt0₊, Mt⁺0₊, μt0₊) are updated to time zero, the result being
written into (Lt0, Mt⁺0, μt0)

Version where B̃ and ã do not depend on try

Compute log tildeρ(0,x_0) for the k-th shape

plot initial and final shape, given by xobs0 and xobsT respectively

update parameter values in GuidedProposal! Q, i.e.
new values are written into Q.target and Q.aux is updated accordingly

Make ObsInfo object

Three cases:
1) obs_atzero=true: this refers to the case of observing one landmark configuration at times 0 and T
2) obs_atzero=false & fixinitmomentato0=false: case of observing multiple shapes at time T,
    both positions and momenta at time zero assumed unknown
3) obs_atzero=false & fixinitmomentato0=true: case of observing multiple shapes at time T,
    positions at time zero assumed unknown, momenta at time 0 are fixed to zero

Simulate guided proposal and compute loglikelihood for one shape
Solve sde inplace and return loglikelihood (thereby avoiding 'double' computations)

Simulate guided proposal and compute loglikelihood (vector version, multiple shapes)

solve sde inplace and return loglikelihood (thereby avoiding 'double' computations)

Stochastic approximation for loglikelihood.

Simulate guided proposal and compute loglikelihood for starting point x0,
guided proposals defined by Q and Wiener increments in W.
Guided proposals are written into X.
Writes vector of loglikelihoods into llout.
Returns sum of loglikelihoods

Compute backward recursion for all shapes and write into field Q.guidrec

update initial state
X:  current iterate of vector of sample paths
Xᵒ: vector of sample paths to write proposal into
W:  current vector of Wiener increments
ll: current value of loglikelihood
x, xᵒ, ∇x, ∇xᵒ: allocated vectors for initial state and its gradient
sampler: either sgd (not checked yet) or mcmc
Q::GuidedProposal!
δ: vector with MALA stepsize for initial state positions (δ[1]) and initial state momenta (δ[2])
updatekernel:  can be :mala_pos, :mala_mom, :rmmala_pos

For fixed Wiener increments and initial state, update parameters by random-walk-MH

update bridges for all shapes using Crank-Nicholsen scheme with parameter ρ (only in case the method is mcmc)
Newly accepted bridges are written into (X,W), loglikelihood on each segment is written into vector ll

update_path!(X,Xᵒ,W,Wᵒ,Wnew,ll,x, Q, ρ, acc_pcn)

Define z(q) = < ∇K̄(q - δ,τ), λ >
Required for Stratonovich -> Ito correction in AHS-model

Suppose one noise field nf
Returns diagonal matrix with noisefield for momentum at point location q (can be vector or Point)

For AHS model compute total noise field on position experienced at a point x.
Useful for plotting purposes.

Example usage:
    σq(Point(0.0, 0.0), nfs)
    σq([0.0; 0.0], nfs)

Suppose one noise field nf
Returns diagonal matrix with noisefield for position at point location q (can be vector or Point)

compute σ(t,x)' y, where y::State
the result is a vector of points that is written to out

compute σ(t,x)' y, where y::State
the result is a vector of points that is written to out

Compute sigma(t,x)' * y where y is a state and sigma is the diffusion coefficient of landmarks returns a vector of points of length P.nfs

Needed for b! in case P is auxiliary process

Needed for b! in case P is auxiliary process

gradient of kernel in hamiltonian

Define ∇z(q) = ∇ < ∇K̄(q - δ,τ), λ >
Required for Stratonovich -> Ito correction in AHS-model