BridgeLandmarks.jl

Landmarks{S,T} <: ContinuousTimeProcess{State{PointF}}

Arguments

• a: Hamiltonian kernel parameter
• c: kernel multiplicate parameter
• n:: Int64 number of landmarks
• db::Float64 square domain bound used for construction of noise fields
• nfstd: standard deviation of noisefields (assumed to be the same for all noisefields)
• nfs: vector of noisefields

LandmarksAux{S,T} <: ContinuousTimeProcess{State{PointF}}

Arguments

• a: Hamiltonian kernel parameter
• c: kernel multiplicate parameter
• xT: State at time T used for constructing the auxiliary process
• n: number of landmarks
• nfs: vector of noisefields

LandmarksAux(P::Landmarks, xT, gramT) = LandmarksAux(P.a,P.c, xT, P.n, P.nfs, gramT)

MarslandShardlow{T} <: ContinuousTimeProcess{State{PointF}}

Arguments

• a: Hamiltonian kernel parameter
• c: kernel multiplicate parameter
• γ: noise level
• λ: mean reversion parameter (heath-bath parameter in Marsland-Shardlow (2017))
• n: number of landmarks

MarslandShardlowAux{S,T} <: ContinuousTimeProcess{State{PointF}}

Arguments

• a: Hamiltonian kernel parameter
• c: kernel multiplicate parameter
• γ: noise level
• λ: mean reversion parameter (heath-bath parameter in Marsland-Shardlow (2017))
• xT: State at time T used for constructing the auxiliary process
• n: number of landmarks
• gramT: Gram-matrix at xT

MarslandShardlowAux(P::MarslandShardlow, xT, gramT) = MarslandShardlowAux(P.a,P.c, P.γ, P.λ, xT, P.n, gramT)

struct summarising settings for updates on mcmc algorithm.

Besides initial values for (c,γ), the values that are most likely to require adjustments are covθprop, δpos, δmom, stdev, db

• ρlowerbound: lower bound for ρ, this is a a number in [0,1). Wiener innovations are updated according to Wnew = ρinit * Wold + sqrt(1-ρinit^2) * Wind where Wind is an independnet Wiener process. At each iteration, ρ is sampled uniformly on [ρlowerbound, 1]

• η: Stepsize cooling function for the adaptation to ensure diminishing adaptation.

Set at n -> min(0.1, 1.0/√n)

• adaptskip: to define guided proposals, we need to plug-in a momentum vector at time T. This is unobserved, and by default initialised as the zero vector. Now every adaptskip iterations, mT is updated to the value obtained in that iteration. In addition, every adaptskip iterations, adaptive tuning of mcmc

parameters is performed.

• σobs: If covariance matrix of extrinsic noise is not specified, then it is set to σobs^2*I

• dt: mesh width using for discretisation (the mesh always gets transformed by a time-change s -> s*(2-s))

• cinit, : initial value for c. Note that we assume Hamltonain kernel x -> c exp(-|x|^2/(2a^2))

• γinit: initial value for gamma which quantifies the amount of intrinsic noise.

• aprior: prior on a

• cprior: prior on c

• γprior: prior on γ

• stdev: For the AHS-model we take noise fields centred at points that are both horizontally and vertically separaeted by a distance that is an integer multiple of 2stdev

• db: assign noisefields within the box [-db[1], db[1]] x -db[2], db[2]

• κ: in case of landmarksmatching, κ quantifies diffusiveness of the prior on the initial momenta

• δpos: stepsize for updating intial positions. Note that the value of δpos is adapted during mcmc-iterations.

• δmom: stepsize for updating intial momenta. Note that the value of δmom is adapted during mcmc-iterations.

• skip_saveITER: skip every skip_saveITER iteration in saving paths

struct summarising settings for updates on mcmc algorithm.

Besides initial values for (c,γ), the values that are most likely to require adjustments are covθprop, δpos, δmom

• ρlowerbound: lower bound for ρ, this is a a number in [0,1). Wiener innovations are updated according to Wnew = ρinit * Wold + sqrt(1-ρinit^2) * W_ind

where W_ind is an independnet Wiener process. At each iteration, ρ is sampled uniformly on [ρlowerbound, 1]

• η: Stepsize cooling function for the adaptation to ensure diminishing adaptation.

Set at n -> min(0.1, 1.0/√n)

• adaptskip: to define guided proposals, we need to plug-in a momentum vector at time T. This is unobserved, and by default initialised as the zero vector. Now every adaptskip iterations, mT is updated to the value obtained in that iteration. In addition, every adaptskip iterations, adaptive tuning of mcmc

parameters is performed.

• σobs: If covariance matrix of extrinsic noise is not specified, then it is set to σobs^2*I

• dt: mesh width using for discretisation (the mesh always gets transformed by a time-change s -> s*(2-s))

• cinit, : initial value for c. Note that we assume Hamltonain kernel x -> c exp(-|x|^2/(2a^2))

• γinit: initial value for gamma which quantifies the amount of intrinsic noise.

• aprior: prior on a

• cprior: prior on c

• γprior: prior on γ

• κ: in case of landmarksmatching, κ quantifies diffusiveness of the prior on the initial momenta

• δpos: stepsize for updating intial positions. Note that the value of δpos is adapted during mcmc-iterations.

• δmom: stepsize for updating intial momenta. Note that the value of δmom is adapted during mcmc-iterations.

• skip_saveITER: skip every skip_saveITER iteration in saving paths

construct_nfs(db::Array{Float64,1}, nfstd, γ)

Construct sequence of Noisefields for AHS model db: domainbound. Vector db[1], db[2], where sources are places on square grid specified by (-db[1]:2nfstd:db[1]) x -db[2]:2nfstd:db[2] 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

gramkernel(q, P; ϵ = 10^(-12))

Gram matrix for kernel with vector of landmarks given by q::Vector(PointF) ϵ*I is added to avoid numerical problems that destroy PSDness of the Gram matrix

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

landmarksmatching(
    xobs0::Array{PointF},xobsT::Array{PointF};
    pars = Pars_ms(),
    updatescheme = [:innov, :mala_mom],
    ITER = 100,
    outdir=@__DIR__,
    Σobs = nothing,
    ainit = nothing,
    printskip=20
)

Arguments

• xobs0: Array of PointF (coordinates of inital shape)
• xobsT: Array of PointF (coordinats of shape at time T)

Optional arguments

• pars: either pars_ms() orpars_ahs()` (this selects the model and default parameter values)
• updatescheme: array of mcmc updates
• ITER: number of iterations
• outdir: path of directory to which output is written
• Σobs: If provided, a length two array, where the first (second) element gives an Array{Unc} elements, where each element specifies the covariance of the extrinsic noise for each landmark at time 0 (final time T). If not provided, Σobs is constructed using the value of pars.σobs
• anit: Hamiltonian kernel parameter. If not provided, defaults to setting ainit = mean(norm.([xobs0[i]-xobs0[i-1] for i in 2:n]))
• printskip: skip every printskip observations in writing output to console

Example:

    # prepare some test data
    n = 7
    xobs0 = [PointF(2.0cos(t), sin(t))/4.0  for t in collect(0:(2pi/n):2pi)[2:end]]
    θ, ψ =  π/6, 0.1
    rot =  SMatrix{2,2}(cos(θ), sin(θ), -sin(θ), cos(θ))
    stretch = SMatrix{2,2}(1.0 + ψ, 0.0, 0.0, 1.0 - ψ)
    shift = [0.5, -1.0]
    xobsT = [rot * stretch * xobs0[i] + shift  for i in 1:n]

    landmarksmatching(xobs0,xobsT; ITER=10)
    landmarksmatching(xobs0,xobsT; ITER=10, pars=BL.Pars_ahs())

landmarksmatching(
    landmarks0::Matrix{Float64},landmarksT::Matrix{Float64};
    pars = Pars_ms(),
    updatescheme = [:innov, :mala_mom],
    ITER = 100,
    outdir=@__DIR__,
    Σobs = nothing,
    ainit = nothing
)

landmarksmatching, where the input does not consist of arrays of points, but both the initial shape and final shape are represented by a matrix. Each row of the matrix gives the coordinates of a landmark.

(t, obsinfo, mT, P, ITER, subsamples, xinit, pars, priorθ, priormom, updatescheme, outdir, printskip)

This is the main function call for doing either MCMC or SGD for landmark models Backward ODEs used are in terms of the LMμ-parametrisation

Arguments

• t: time grid
• obsinfo: struct of type ObsInfo
• mT: vector of momenta at time T used for constructing guiding term
• P: target process
• ITER: number of iterations
• subsamples: vector of indices of iterations that are to be saved
• xinit: initial value on state at time 0 (x0)
• pars: tuning pars for updates (see ?Pars_ms() and ?Pars_ahs() for documentation)
• priorθ: prior on θ = (a,c,γ)
• priormom: prior on initial momenta
• updatescheme: vector specifying mcmc-updates
• outdir output directory for animation
• printskip skip every printskip iterations in writing output to console

Returns:

• Xsave: saved iterations of all states at all times in tt_
• parsave: saved iterations of all parameter updates
• initendstates_save: saved iterations of initial and endstate
• accinfo: acceptance indicators of mcmc-updates
• δ: value of (δpos, δmom) af the final iteration of the algorithm (these are stepsize parameters for updating initial positions and momenta respectively)
• ρ: value of ρinit at the final iteration of the algorithm
• δa: value of δa at the final iteration of the algorithm (step size for updating Hamiltonian kernel parameter a)

show_updates()

Convenience function to show the available implemted updates.

template_estimation(
    xobsT::Array{PointF};
    pars = Pars_ms(),
    updatescheme = [:innov, :rmmala_pos, :parameter],
    ITER = 100,
    outdir=@__DIR__,
    Σobs = nothing,
    ainit = nothing
    xinitq = nothing
)

Arguments

• xobsT: an array, where each element of the array gives the coordinates of a shape (i.e. the landmarks, which

are represented as an array of PointF)

Optional arguments

• pars: either pars_ms()orpars_ahs()`` (this selects the model and default parameter values)
• updatescheme: array of mcmc updates
• ITER: number of iterations
• outdir: path of directory to which output is written
• Σobs: If provided, a length two array, where the first (second) element gives an Array{Unc} elements, where each element specifies the covariance of the extrinsic noise for each landmark at time 0 (final time T). If not provided, Σobs is constructed using the value of pars.σobs
• anit: Hamiltonian kernel parameter. If not provided, defaults to setting ainit = mean(norm.([xobs0[i]-xobs0[i-1] for i in 2:n]))
• xinitq::Array{PointF}: Initialisation for the shape. If not provided xobsT[1] will be taken for initialisation.
• printskip: skip every printskip observations in writing output to console

Example:

``` # Make test example data set n = 5 nshapes = 7 q0 = [PointF(2.0cos(t), sin(t)) for t in collect(0:(2pi/n):2pi)[2:end]] # initial shape is an ellipse x0 = State(q0, randn(PointF,n)) xobsT = Vector{PointF}[]

T = 1.0; dt = 0.01; t = 0.0:dt:T
Ptrue = MarslandShardlow(2.0, 0.1, 1.7, 0.0, n)
for k in 1:nshapes
        Wf, Xf = BridgeLandmarks.landmarksforward(t, x0, Ptrue)
        push!(xobsT, Xf.yy[end].q + σobs * randn(PointF,n))
end

template_estimation(xobsT)

`

template_estimation(
    landmarksT::Array{Matrix{Float64}};
    pars = Pars_ms(),
    updatescheme = [:innov, :mala_mom],
    ITER = 100,
    outdir=@__DIR__,
    Σobs = nothing,
    ainit = nothing
)

GuidRecursions{TL,TM⁺,TM, Tμ, TH, TLt0, TMt⁺0, Tμt0}

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

Arguments

Suppose t is the specified (fixed) time grid. Then the elements of the struct are:

• Lt: matrices L
• Mt⁺`: matrices M⁺ (inverses of M)
• M: matrices M
• μ`: vectors μ
• Ht: matrices H, where H = L' M L
• Lt0: L(0) (so obtained from L(0+) after gpupdate step incorporating observation xobs0)
• Mt⁺0: M⁺(0) (so obtained from M⁺(0+) after gpupdate step incorporating observation xobs0)
• μt0: μ(0) (so obtained μ(0+) after gpupdate step incorporating observation xobs0)

GuidedProposal{T,Ttarget,Taux,TL,Txobs0,TxobsT,Tnshapes,TmT,F} <: ContinuousTimeProcess{T}

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

Arguments

target::Ttarget                 # target diffusion P
aux::Vector{Taux}               # auxiliary diffusion for each shape (Ptilde for each shape)
tt::Vector{Float64}             # grid of time points on single segment (S,T]
xobs0::Txobs0                   # observation at time 0
xobsT::Vector{TelxobsT}         # vector of observations at time T
guidrec::Vector{TL}             # guided recursions on grid tt
nshapes::Int64                  # number of shapes
mT::Vector{TmT}                 # vector of artificial momenta used for constructing auxiliary process (one for each shape)
endpoint::F

Example (for landmarksmatching)

    n = 5
    P = MarslandShardlow(1.0, 2.0, 0.5, 0.0, n)
    xT = State(rand(PointF,n),rand(PointF,n))
    Paux = [BridgeLandmarks.auxiliary(P,xT)]
    tt = collect(0:0.05:1.0)
    xobs0 = rand(PointF,n)
    xobsT = [rand(PointF,n)]
    σobs = 0.01
    Σobs = fill([σobs^2 * one(UncF) for i in 1:n],2)
    obsinfo = BridgeLandmarks.set_obsinfo(xobs0,xobsT,Σobs, true,false)
    gr = [BridgeLandmarks.init_guidrec(tt,obsinfo)]
    nshapes = 1
    mT = [rand(PointF,n)]
    Q = BridgeLandmarks.GuidedProposal(P,Paux,tt,xobs0,xobsT,gr,nshapes,mT)

struct containing information on the observations

We assue observations V0 and VT, where

• V0 = L0 * X0 + N(0,Σ0)
• VT = LT * X0 + N(0,ΣT)

In addition, μT is a vector of zeros (for initialising the backward ODE on μ) (possibly remove later)

Base.:*(H::InverseCholesky, X)

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

Input are L and X, output is Y

y = H X 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

Base.:*(H::InverseCholesky, x::State)

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)

Base.show(io::IO, state::State)

Good display of variable of type State

Bridge.B!(t,X,out, Paux::LandmarksAux, xT)

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

Bridge.B!(t,X,out, Paux::MarslandShardlowAux)

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

Bridge.B(t, Paux::LandmarksAux, xT)

Compute tildeB(t) for landmarks auxiliary process

Bridge.B(t, Paux::MarslandShardlowAux)

Compute tildeB(t) for landmarks auxiliary process

_b!((i,t), x::State, out::State, Q::GuidedProposal,k)

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

Bridge.a!(t, x_, out, P::Union{Landmarks,LandmarksAux})

Bridge.a(t,  P::Union{MarslandShardlow, MarslandShardlowAux})

Returns matrix a(t) for Marsland-Shardlow model

auxiliary(Q::GuidedProposal,k::Int64) = Q.aux[k]

Extract auxiliary process of k-th shape.

auxiliary(P::Union{MarslandShardlow, Landmarks}, xT, gramT)

Construct auxiliary process corresponding to P and xT

Bridge.b!(t, x, out, Paux::LandmarksAux)

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

Bridge.b!(t, x, out, Paux::MarslandShardlowAux)

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

Bridge.b!(t, x, out, P::MarslandShardlow)

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

constdiff(Q::GuidedProposal)

If true, both the target and auxiliary process have constant diffusion coefficient.

target(Q::GuidedProposal) = Q.target

σ!(t, x, dw, out, Q::GuidedProposal)

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

Bridge.σ!(t, x_, dm, out, P::Union{Landmarks,LandmarksAux})

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

Bridge.σ!(t, x, dm, out, P::MarslandShardlowAux, xT)

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

Bridge.σ!(t, x, dm, out, P::MarslandShardlow)

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

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

K̄(q,τ)

K̄(q,τ) = exp(-Bridge.inner(q)/(2*τ^2)) Kernel for noisefields of AHS-model

_r!((i,t), x::State, out::State, Q::GuidedProposal,k)

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

adapt_pospar_step!(δ, n, accinfo, η, update; adaptskip = 15, targetaccept=0.5)

Adapt the step size for updates, where tuning par has to be positive. The adaptation is multiplication/divsion by exp(η(n)).

adaptpcnstep(n, accpcn, ρ, nshapes, η; adaptskip = 15, targetaccept=0.5)

Adjust pcN-parameter ρ adaptive every adaptskip steps. For that ρ is first tranformed to (-∞, ∞), then updated, then mapped back to (0,1)

adjust_to_newpars(Q::GuidedProposal,(aᵒ,cᵒ,γᵒ),obsinfo)

Provide new parameter values for GuidedProposal Q, these are written into fields target and aux. Returns a new instance of GuidedProposal, adjusted to the new set of parameters.

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

amul(t, x::State, xin::State, P::Union{MarslandShardlow, MarslandShardlowAux})

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

amul(t, x::State, xin::Vector{<:Point}, P::Union{Landmarks,LandmarksAux})

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 σ)

cholinverse!(L, x)

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

construct_gp_xobsT(Q, xobsTᵒ)

Update xobsTᵒ into auxiliary process of Q, following by recomputing the backwards ODEs

convert_samplepath(X)

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

convert_samplepath(Xvec::Vector)

Useful for storage of a samplepath of states Ordering is as follows:

1. shape
2. time
3. landmark nr
4. for each landmark: q1, q2 p1, p2

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

function copypaths!(X,Xᵒ)

Write Xᵒ into X

deepvalue(x)

If x is a vector of Float64, x is returned. If x is a vector of Dual-numbers, its Float64 part is returned.

deepvalue(x::State)

Extract Float64 part of elements in State (so in case of Dual numbers, derivative part is dropped)

extract_initial_and_endstate(iteratenr::Int64, Xpath)

Extract from a sample path its states at time zero and time T (endtime)

iteratnr:: integer that indicates iteration.nr in mcmc-algorithm Xpath:: Bridge.Samplepath

Returns a matrix with first column iteratenr repeated, second column state at time zero, third column state a time T

The second and third column as ordered as follows (in case of two dimensional landmarks): landmark 1 pos1 landmark 1 pos2 landmark 1 mom1 landmark 1 mom2

landmark 2 pos1 landmark 2 pos2 landmark 2 mom1 landmark 2 mom2

...

getpars(Q::GuidedProposal)

Extract parameters from GuidedProposal Q, that is, (a,c,γ)`

gp!(::LeftRule,  Xᵒ, x0, W, Q::GuidedProposal,k; skip = sk, ll0 = true)

Simulate guided proposal as specified in Q and compute loglikelihood for one shape, starting from x0, using Wiener increments W

Write into

Xᵒ

Returns

logliklihood.

gp!(::LeftRule,  X::Vector, q, p , W, Q::GuidedProposal; skip = sk, ll0 = true)

Simulate guided proposal and compute loglikelihood (vector version, multiple shapes) q and p make up a state (could in fact be turned into the intial state by merge_state(q,p))

gp!(::LeftRule,  X, x0, W, Q::GuidedProposal; skip = sk, ll0 = true)

Simulate guided proposal and compute loglikelihood (vector version, multiple shapes) x0 is assumed to be of type State

gp_update!(Lt0₊, Mt⁺0₊::Array{Pnt,2}, μt0₊, (L0, Σ0, xobs0), Lt0, Mt⁺0, μt0) where Pnt

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)

gram_matrix!(out, q, P)

Write into out the Gram matrix for kernel with vector of landmarks given by q::Vector(PointF)

gram_matrix(q, P)

Gram matrix for kernel with vector of landmarks given by q::Vector(PointF)

guidingbackwards!(::Lm, t, (Lt, Mt⁺, μt), Paux, obsinfo; implicit=true, lowrank=false)

Solve backwards recursions in L, M, μ parametrisation on grid t

Arguments

• t: time grid
• (Lt, Mt⁺, μt): containers to write the solutions into
• Paux: auxiliary process
• obsinfo: of type ObsInfo containing information on the observations
• implicit: if true an implicit Euler backwards scheme is used (else explicit forward)

Case lowrank=true still gives an error: fixme!

hamiltonian(x, P)

Hamiltonian for deterministic part of landmarks model

hamiltonian(x::NState, P::MarslandShardlow)

Returns Hamiltonian at x for deterministic system (no Wiener noise)

kernel(q, P::LandmarkModel)

kernel in Hamiltonian: P.c * exp(-Bridge.inner(q)/(2*P.a^2))

Version where B̃ and ã do not depend on try

lρtilde(x0, Q,k)

Compute log ρ̃(0,x_0,k), where k indexes shape

merge_state(q,p)

Merge q and p parts of state to an element of type NState.

Note that merge does the inverse operation

Example

y = State(rand(PointF,3),rand(PointF,3))
q, p = split_state(y)
yy = merge_state(q,p)
yy-y

onemask(x::Number) = one(x)

set_guidrec(Q::GuidedProposal, gr) = GuidedProposal(Q.target, Q.aux, X.tt, Q.xobs, Q.xobsT, gr,Q.nshapes, Q.endpoint)

Update guidrec field of Q to gr

set_obsinfo(xobs0,xobsT,Σobs, obs_atzero::Bool,fixinitmomentato0::Bool)

Constructor for ObsInfo.

Arguments

• n: number of landmarks
• obs_atzero: Boolean, if true, the initial configuration is observed
• fixinitmomenta0: Boolean, if true, the initial momenta are fixed to zero
• Σobs: 2-element array where Σobs0 = Σobs[1] and ΣobsT = Σobs[2] Both Σobs0 and ΣobsT are arrays of length n of type UncF that give the observation covariance matrix on each landmark
• xobs0: in case obs_atzero=true, it is provided and passed through; in other cases it is constructed such that the backward ODEs are initialised correctly.

Note that there are three cases:

• obs_atzero=true: this refers to the case of observing one landmark configuration at times 0 and T
• obs_atzero=false & fixinitmomentato0=false: case of observing multiple shapes at time T, both positions and momenta at time zero assumed unknown
• 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

slogρ!(q, p, Q, W, X, priormom,llout)

Main use of this function is to get gradient information of the loglikelihood with respect ot the initial state.

Arguments

• q and p make up the initial state (which could be obtained from x0 = merge_state(q,p))
• Q: GuidePropsoal!
• W: vector of innovations
• X: sample path
• priormom: prior on the initial momenta
• llout: vector where loglikelihoods for each shape are written into

Writes into

• The guided proposal is written into X.
• The computed loglikelihood for each shape is written into llout

Returns

• loglikelihood (summed over all shapes) + the logpriordensity of the initial momenta

slogρ_mom!(q, Q, W, X,priormom, llout) = (p) -> slogρ!(q, p , Q, W, X, priormom, llout)

slogρ_pos!(p, Q, W, X,priormom, llout) = (q) -> slogρ!(q, p , Q, W, X, priormom, llout)

split_state(x::NState)

Split a state into its q and p parts, where the latter are reinterpreted as vectors.

Note that merge does the inverse operation

Example

y = State(rand(PointF,3),rand(PointF,3))
q, p = split_state(y)
yy = merge_state(q,p)
yy-y

tame(x) = x/(I + norm(x))

Useful for tamed version of Langevin algorithm

trun(x,h) = x/max(1.0,0.5*h*norm(x))

useful function for truncated version of mala

update_cyclicmatching(X, ll,obs_info, Xᵒ, W, Q)

writes into / modifies

• X
• ll

Returns

Q, obsinfo, accept

update_guidrec!(Q, obsinfo)

Q::GuidedProposal obsinfo::ObsInfo

Computes backwards recursion for (L,M⁺,μ), including gp-update step at time 0. Next, the guidrec field or Q is updated. Returns Q.

update_initialstate!(X,Xᵒ,W,ll,x,xᵒ,∇, ∇ᵒ, Q::GuidedProposal, δ, update,priormom, (dK, inv_dK), At)

Arguments

• 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 (both actual and proposed values)
• sampler: either sgd (not checked yet) or mcmc
• Q: GuidedProposal
• δ: vector with MALA stepsize for initial state positions (δ[1]) and initial state momenta (δ[2])
• update: can be :mala_pos, :mala_mom, :rmmalapos,:rmmalamom,:rmrw_mom
• priormom: prior on initial momenta
• (dK, inv_dK): gramkernel at x0.q and its inverse

Writes into / modifies

• X: landmarks paths
• x: initial state of the landmark paths
• `∇: gradient of loglikelihood with respect to initial state x
• ll: a vector with as k-th element the loglikelihood for the k-th shape

Returns

• 0/1 indicator (reject/accept)

update_mT(Q, mTv, obsinfo)

Update State vector of auxiliary process for each shape. For the k-th shape, the momentum gets replaced with mTv[k]

Example

See documentation for GuidedProposal to contruct an example instance, say Q

```julia mTv = [zeros(PointF,5)]

BridgeLandmarks.update_mT!(Q, mTv, obsinfo)

update_pars!(X, Xᵒ,W, Q, x, ll, priorθ, covθprop, obsinfo, At)

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

Writes into / modifies

X: landmark paths under landmarks process with parameter θ ll: vector of loglikelihoods (one for each shape)

Returns

Q, accept: where accept is 1/0 according to accept/reject in the MH-step.

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

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 All Wiener increments are always updated.

Write into

• X: diffusion paths
• W: innovations
• ll: vector of loglikehoods (k-th element is for the k-th shape)

write_acc(accinfo,outdir)

Write dataframe with acceptance numbers of updatesteps to a csv file

Write info to txt file

write_initendstates(initendstates_save)

initendstatessave is an Array of matrices, each matrix corresponding to one iteration, with ordering as specified in `extractinitialandendstate`

write_mcmc_iterates(Xsave, tt_, n, nshapes, subsamples, outdir)

Write bridge iterates to file

Arguments:

• 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

write_noisefields(P,outdir)

Write noisefields to file (empty in case of MS-model)

write_observations(xobs0, xobsT, n, nshapes, outdir)

Write observations to csv file

write_params(parsave,subsamples,outdir)

Write parameter iterates to file

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

σp(q, p, nf::Noisefield) = -Diagonal(p .* nf.γ .* ∇K̄(q - nf.δ,nf.τ))

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

σq(q, nf::Noisefield) = Diagonal(nf.γ * K̄(q - nf.δ,nf.τ))

Suppose one noise field nf Returns diagonal matrix with noisefield for position 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)

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

σt!(t, x_, y::State{Pnt}, out, P::MarslandShardlow) where Pnt

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

σ̃(t,  Paux::LandmarksAux, xT)

Return matrix σ̃(t) for LandmarksAux

σ̃(t,  P::Union{MarslandShardlow, MarslandShardlowAux})

Return sparse matrix matrix σ̃(t)

∇K̄(q, qT, τ)

Needed for b! in case P is auxiliary process

∇K̄(q,τ)

Gradient of kernel for noisefields

Needed for b! in case P is auxiliary process

∇kernel(q, P::LandmarkModel)

gradient of kernel in hamiltonian

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