From 52eebc189f180220286c86951596b2a16cb11026 Mon Sep 17 00:00:00 2001 From: fmeulen Date: Fri, 3 May 2019 13:32:21 +0200 Subject: [PATCH 1/3] added lowrank functionality for computing tilde-a in backward recursion --- landmarks/lmguiding.jl | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/landmarks/lmguiding.jl b/landmarks/lmguiding.jl index e63b365..0d0a293 100644 --- a/landmarks/lmguiding.jl +++ b/landmarks/lmguiding.jl @@ -26,14 +26,26 @@ end struct Lm end function guidingbackwards!(::Lm, t, (Lt, Mt⁺, μt), Paux, (Lend, Mend⁺, μend)) - Mt⁺[end], Lt[end] = Σ, L + Mt⁺[end] .= Σ + Lt[end] .= L BB = Matrix(Bridge.B(0, Paux)) # does not depend on time - aa = Matrix(Bridge.a(0, Paux)) # does not depend on time + println("computing ã and its low rank approximation:") + # various ways to compute ã (which does not depend on time); + # low rank appoximation really makes sense here +# @time aa = Matrix(Bridge.a(0, Paux)) # vanilla, no lr approx +# @time aalr = pheigfact(deepmat(Matrix(Bridge.a(0, Paux)))) # low rank approx default +# @time aalr = pheigfact(deepmat(Matrix(Bridge.a(0, Paux))),rank=400) # fix rank + @time aalr = pheigfact(deepmat(Matrix(Bridge.a(0, Paux))), rtol=1e-7) # control accuracy of lr approx + println("Rank ",size(aalr[:vectors],2), " approximation to ã") + sqrt_aalr = deepmat2unc(aalr[:vectors] * diagm(0=> sqrt.(aalr[:values]))) + β = vec(Bridge.β(0,Paux)) # does not depend on time for i in length(t)-1:-1:1 dt = t[i+1]-t[i] - Lt[i] .= Lt[i+1] * (I + BB * dt) - Mt⁺[i] .= Mt⁺[i+1] + Lt[i+1]* aa * Matrix(Lt[i+1]') * dt +# Lt[i] .= Lt[i+1] * (I + BB * dt) # explicit + Lt[i] .= Lt[i+1]/(I - dt* BB) # implicit, similar computational cost +# Mt⁺[i] .= Mt⁺[i+1] + Lt[i+1]* aa * Matrix(Lt[i+1]') * dt + Mt⁺[i] .= Mt⁺[i+1] + Bridge.outer(Lt[i+1] * sqrt_aalr) * dt μt[i] .= μt[i+1] + Lt[i+1] * β * dt end (Lt[1], Mt⁺[1], μt[1]) @@ -45,7 +57,7 @@ auxiliary(Q::GuidedProposall!) = Q.aux constdiff(Q::GuidedProposall!) = constdiff(target(Q)) && constdiff(auxiliary(Q)) -function Bridge._b!((i,t), x::State, out::State, Q::GuidedProposall!) +function _b!((i,t), x::State, out::State, Q::GuidedProposall!) Bridge.b!(t, x, out, Q.target) out .+= amul(t,x,Q.Lt[i]' * (Q.Mt[i] *(Q.xobs-Q.μt[i]-Q.Lt[i]*vec(x))),Q.target) out @@ -60,7 +72,8 @@ function _r!((i,t), x::State, out::State, Q::GuidedProposall!) end function guidingterm((i,t),x::State,Q::GuidedProposall!) - Bridge.b(t,x,Q.target) + amul(t,x,Q.Lt[i]' * (Q.Mt[i] *(Q.xobs-Q.μt[i]-Q.Lt[i]*vec(x))),Q.target) + #Bridge.b(t,x,Q.target) + + amul(t,x,Q.Lt[i]' * (Q.Mt[i] *(Q.xobs-Q.μt[i]-Q.Lt[i]*vec(x))),Q.target) end """ Returns the guiding terms a(t,x)*r̃(t,x) along the path of a guided proposal From 390c1fd90ecfce52ab967b60c8140355958907cf Mon Sep 17 00:00:00 2001 From: fmeulen Date: Fri, 3 May 2019 13:33:45 +0200 Subject: [PATCH 2/3] essentially only added 'using LowRankApprox' --- landmarks/lmpar.jl | 35 +++++++++++++++++++---------------- 1 file changed, 19 insertions(+), 16 deletions(-) diff --git a/landmarks/lmpar.jl b/landmarks/lmpar.jl index bd7a98a..98dc171 100644 --- a/landmarks/lmpar.jl +++ b/landmarks/lmpar.jl @@ -12,6 +12,7 @@ using CSV using RCall using Base.Iterators using SparseArrays +using LowRankApprox models = [:ms, :ahs] model = models[2] @@ -22,14 +23,14 @@ partialobs = true const d = 2 const itostrat = false -n = 75 # nr of landmarks +n = 100 # nr of landmarks -θ = π/3# π/6 0#π/5 # angle to rotate endpoint +θ = 2π/5# π/6 0#π/5 # angle to rotate endpoint -σobs = 10^(-3) # noise on observations +σobs = 10^(-2) # noise on observations println(model) -T = 0.3#1.0#0.5 +T = 0.75#1.0#0.5 t = 0.0:0.01:T # time grid #Random.seed!(5) @@ -43,12 +44,12 @@ include("lmguiding.jl") ### Specify landmarks models a = 3.0 # the larger, the stronger landmarks behave similarly λ = 0.0; #= not the lambda of noise fields =# γ = 8.0 -db = 3.0 # domainbound -nfstd = 1.0 # tau , widht of noisefields +db = 6.0 # domainbound +nfstd = 2.0 # tau , widht of noisefields r1 = -db:nfstd:db r2 = -db:nfstd:db nfloc = Point.(collect(product(r1, r2)))[:] -nfscales = [.05Point(1.0, 1.0) for x in nfloc] # intensity +nfscales = [.1Point(1.0, 1.0) for x in nfloc] # intensity nfs = [Noisefield(δ, λ, nfstd) for (δ, λ) in zip(nfloc, nfscales)] Pms = MarslandShardlow(a, γ, λ, n) @@ -69,7 +70,7 @@ W = SamplePath(t, [copy(w0) for s in t]) sample!(W, Wiener{Vector{StateW}}()) # specify initial landmarks configuration -q0 = [Point(2.5cos(t), sin(t)) for t in (0:(2pi/n):2pi)[1:n]] #q0 = circshift(q0, (1,)) +q0 = [Point(2.0cos(t), sin(t)) for t in (0:(2pi/n):2pi)[1:n]] #q0 = circshift(q0, (1,)) p_ = 2*Point(-0.1, 0.1) p0 = [p_ for i in 1:n] # #p0 = [randn(Point) for i in 1:n] @@ -81,17 +82,19 @@ println("Solve for forward provess:") @time solve!(EulerMaruyama!(), X, x0, W, P) #@time solve!(StratonovichHeun!(), X, x0, W, P) -tc(t,T) = t.*(2-t/T) +tc(t,T) = t.* (2 .-t/T) tt_ = tc(t,T)#tc(t,T)# 0:dtimp:(T) - +#tt_ = t # observe positions without noise v0 = q(X.yy[1]) -rot = SMatrix{2,2}(cos(θ), sin(θ), -sin(θ), cos(θ)) -vT = [rot * X.yy[end].q[i] for i in 1:P.n ] +rot = SMatrix{2,2}(cos(θ), sin(θ), -sin(θ), 2*cos(θ)) +#vT = [rot * X.yy[end].q[i] for i in 1:P.n ] +trans = SMatrix{2,2}(1.5, 1.0, 0.0, 1.0) +vT = [rot * trans * X.yy[end].q[i] for i in 1:P.n ] + -#extractcomp(v,i) = map(x->x[i], v) #################### # solve backward recursion on [0,T] if partialobs==true @@ -100,9 +103,9 @@ if partialobs==true #Σ = 10 * reshape(rand(Unc,n^2),n,n) μend = zeros(Point,P.n) xobs = vT - Pahsaux = LandmarksAux(Pahs, State(vT, zero(vT))) - Pmsaux = MarslandShardlowAux(Pms, State(vT, zero(vT))) - #Pahsaux = LandmarksAux(Pahs, State(vT, 10*rand(Point,Pahs.n))) + mT = zero(vT)#rand(Point, n)# + Pahsaux = LandmarksAux(Pahs, State(vT, mT)) + Pmsaux = MarslandShardlowAux(Pms, State(vT, mT)) else # full observation case L = [(i==j)*one(Unc) for i in 1:2n, j in 1:2n] From 2825f0d33140552d10756135554f81739822c8b9 Mon Sep 17 00:00:00 2001 From: fmeulen Date: Fri, 3 May 2019 13:34:41 +0200 Subject: [PATCH 3/3] speed up computing a(t,x) by using that it is symmetric --- landmarks/models.jl | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/landmarks/models.jl b/landmarks/models.jl index 009a4ed..33cff93 100644 --- a/landmarks/models.jl +++ b/landmarks/models.jl @@ -386,6 +386,26 @@ end Bridge.a(t, x, P::Union{MarslandShardlow, MarslandShardlowAux}) = Bridge.a(t, P) +# function Bridge.a(t, x_, P::Union{Landmarks,LandmarksAux}) +# if P isa Landmarks +# x = x_ +# else +# x = P.xT +# end +# out = zeros(Unc,2P.n,2P.n) +# for i in 1:P.n +# for k in 1:P.n +# for j in 1:length(P.nfs) +# out[2i-1,2k-1] += σq(q(x,i),P.nfs[j]) * σq(q(x, k),P.nfs[j])' +# out[2i-1,2k] += σq(q(x,i),P.nfs[j]) * σp(q(x,k),p(x,k),P.nfs[j])' +# out[2i,2k-1] += σp(q(x,i),p(x,i),P.nfs[j]) * σq(q(x,k),P.nfs[j])' +# out[2i,2k] += σp(q(x,i),p(x,i),P.nfs[j]) * σp(q(x,k),p(x,k),P.nfs[j])' +# end +# end +# end +# out +# end + function Bridge.a(t, x_, P::Union{Landmarks,LandmarksAux}) if P isa Landmarks x = x_ @@ -394,18 +414,24 @@ function Bridge.a(t, x_, P::Union{Landmarks,LandmarksAux}) end out = zeros(Unc,2P.n,2P.n) for i in 1:P.n - for k in 1:P.n + for k in i:P.n for j in 1:length(P.nfs) out[2i-1,2k-1] += σq(q(x,i),P.nfs[j]) * σq(q(x, k),P.nfs[j])' out[2i-1,2k] += σq(q(x,i),P.nfs[j]) * σp(q(x,k),p(x,k),P.nfs[j])' out[2i,2k-1] += σp(q(x,i),p(x,i),P.nfs[j]) * σq(q(x,k),P.nfs[j])' out[2i,2k] += σp(q(x,i),p(x,i),P.nfs[j]) * σp(q(x,k),p(x,k),P.nfs[j])' end + end + end + for i in 2:2P.n + for k in 1:i-1 + out[i,k] = out[k,i] end end out end + Bridge.a(t, P::LandmarksAux) = Bridge.a(t, 0, P) #Bridge.a(t, P::Union{Landmarks,LandmarksAux}) = Bridge.a(t, P.xT, P::Union{Landmarks,LandmarksAux})