Merge Newlandmarks

landmarks/lmguiding.jl

@@ -35,14 +35,26 @@ end
 struct Lm  end
 
 function guidingbackwards!(::Lm, t, (Lt, Mt⁺, μt), Paux, (Lend, Mend⁺, μend))
-    Mt⁺[end], Lt[end], μt[end] = Σend, Lend, μend
-    BB = deepmat(Bridge.B(0, Paux)) # does not depend on time
-    aa = deepmat(Bridge.a(0, Paux)) # does not depend on time
-    β = deepvec(Bridge.β(0,Paux)) # does not depend on time
+    Mt⁺[end] .= Σ
+    Lt[end] .= L
+    BB = Matrix(Bridge.B(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])
@@ -54,7 +66,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
@@ -69,7 +81,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

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]

landmarks/models.jl

@@ -410,6 +410,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_
@@ -429,6 +449,11 @@ function Bridge.a(t, x_, P::Union{Landmarks,LandmarksAux})
                 out[2i,2k-1] += a21 * a12'
                 out[2i,2k] += a21 * a22'
             end
+      end
+    end
+    for i in 2:2P.n
+        for k in 1:i-1
+            out[i,k] = out[k,i]
         end
     end
     for i in 2:2P.n
@@ -439,6 +464,7 @@ function Bridge.a(t, x_, P::Union{Landmarks,LandmarksAux})
     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})