Merge pull request #24 from mschauer/weeksworth

Recent updates to the code

REQUIRE

@@ -2,4 +2,5 @@ julia 0.7
 Polynomials
 Distributions
 StaticArrays
-RecipesBase
+RecipesBase
+Colors

docs/src/library.md

@@ -162,5 +162,8 @@ Bridge.sizedtype
 Bridge.piecewise
 Bridge.BridgePre!
 Bridge.aeuler
+Bridge.MeanCov
+Bridge.upsample
+Bridge.viridis 
 ```
 

extra/makie.jl

@@ -5,6 +5,15 @@ import Bridge: mcsvd3, visualize_uncertainty
 convert_arguments(P::Type{<:Union{Lines,Scatter}}, X::SamplePath{<:AbstractVector}) = convert_arguments(P, X.yy)
 convert_arguments(P::Type{<:Union{Lines,Scatter}}, X::SamplePath{<:Real}) = convert_arguments(P, X.tt, X.yy)
 
+function band!(scene, x, ylower, yupper; nargs...)
+    n = length(x)
+    coordinates = [x ylower; x yupper]
+    ns = 1:n-1
+    ns2 = n+1:2n-1
+    connectivity = [ns ns .+ 1 ns2;
+                    ns2 ns2 .+ 1 ns .+ 1]
+    mesh!(scene, coordinates, connectivity; shading = false, nargs...)
+end
 
 """
     mcsvd3(mcstates) -> mean, q, sv

project/partialbridge.jl

@@ -1,5 +1,3 @@
-cd("/Users/Frank/.julia/dev/Bridge")
-
 using Bridge, StaticArrays, Distributions
 using Test, Statistics, Random, LinearAlgebra
 using Bridge.Models
@@ -104,7 +102,7 @@ for iter in 1:iterations
 end
 
 # write mcmc iterates to csv file 
-f = open("/Users/Frank/Dropbox/DiffBridges/Rcode/integrated_diff/iterates.csv","w")
+f = open("output/integrated_diff/iterates.csv","w")
 head = "iteration, time, component, value \n"
 write(f, head)
 iterates = [Any[s, tt[j], d, XX[i].yy[j][d]] for d in 1:2, j in 1:length(X), (i,s) in enumerate(subsamples) ][:]
@@ -136,9 +134,9 @@ library(tidyverse)
 library(grDevices)
 
 theme_set(theme_minimal())
-setwd("~/Dropbox/DiffBridges/Rcode/integrated_diff")
+setwd("output/integrated_diff")
 
-d <- read.csv("~/Dropbox/DiffBridges/Rcode/integrated_diff/iterates.csv")
+d <- read.csv("output/integrated_diff/iterates.csv")
 dsub <- d[d$iteration %in% c(seq(0,1000,by=200),seq(40000,50000,by=2000)),]
 
 
@@ -148,10 +146,8 @@ dsub %>% ggplot(aes(x=time,y=value,colour=iteration)) +
   facet_wrap(~component,ncol=1,scales='free_y')+
   scale_colour_gradient(low='yellow',high='Darkred')
 dev.off()    
-
-#setwd("/Users/Frank/.julia/dev/Bridge")
 """
-cd("/Users/Frank/.julia/dev/Bridge")
+
 # import variables in Julia 
 @rget d
 

src/Bridge.jl

@@ -41,6 +41,7 @@ include("chol.jl")
 
 using StaticArrays
 using Distributions
+using Colors
 
 
 
@@ -94,6 +95,7 @@ hasbi(::Any) = false
 hasai(::Any) = false
 
 include("expint.jl")
+#include("setcol.jl")
 
 include("fsa.jl")
 include("gaussian.jl")
@@ -116,6 +118,7 @@ include("guip!.jl")
 include("partialbridge.jl")
 
 include("euler.jl")
+include("sde.jl")
 
 include("sde!.jl")
 include("levy.jl")

src/euler.jl

@@ -61,6 +61,9 @@ struct EulerMaruyama! <: SDESolver
 end
 
 
+struct EulerMaruyamaWithIndex! <: SDESolver
+end
+
 
 function solve!(::StochasticHeun, Y, u, W::SamplePath, P::ProcessOrCoefficients)
 

src/mclog.jl

@@ -106,4 +106,69 @@ function mcmarginalstats(states)
         append!(Xstd, map(x->sqrt.(diag(x)), vv))
     end 
     Xmean, Xstd
-end
+end
+
+"""
+    MeanCov(itr)
+
+Iterator interface for online mean and covariance
+Iterates are triples `mean, λ, cov/λ`  
+   
+"""
+struct MeanCov{T}
+    iter::T
+end
+import Base: length, IteratorSize
+length(x::MeanCov) = length(x.iter)
+IteratorSize(::Type{MeanCov{T}}) where {T} = IteratorSize(T) isa Base.HasShape ?  Base.HasLength() : IteratorSize(T)
+
+
+@inline function iterate(mc::MeanCov) 
+    u = iterate(mc.iter)
+    if u === nothing
+        return nothing
+    end
+    x, state = u
+    n = 0
+    delta = copy(x)
+    m = delta/(n+1)
+
+    m2 = outer(delta, x - m)
+    (m, 1/n, m2), (m, m2, delta, n + 1, state)
+end
+
+
+iterate(mc::MeanCov, mcstate) = iterate_(mc, ismutable(typeof(mcstate[1])), mcstate...)    
+
+function iterate_(mc::MeanCov, ::Val{false}, m, m2, delta, n, state ) 
+    u = iterate(mc.iter, state)
+    if u === nothing
+        return nothing
+    end
+    x, state = u
+    delta = x - m
+    m  += delta/(n+1)
+    m2 += outer(delta, x - m)
+    (m, 1/n, m2), (m, m2, delta, n + 1, state)
+end
+
+
+function iterate_(mc::MeanCov, ::Val{true}, m, m2, delta, n, state ) 
+    u = iterate(mc.iter, state)
+    if u === nothing
+        return nothing
+    end
+    x, state = u
+
+    for i in eachindex(m)
+        delta[i] = x[i] - m[i]
+        m[i] += (delta[i])/(n+1)
+    end
+    for i in eachindex(m)
+        for j in eachindex(m)
+            m2[i,j] += outer(delta[i], x[j]-m[j])
+        end
+    end
+    (m, 1/n, m2), (m, m2, delta, n + 1, state)
+end
+

src/misc.jl

@@ -111,6 +111,36 @@ function piecewise(Y::SamplePath, tend = Y.tt[end])
     tt, repeat(Y.yy, inner=2)
 end
 
+"""
+    upsample(x, td, t) 
+
+If `x` is piecewise constant with jumps at `td`, return values of `x`
+at times `t`.
+"""
+upsample(x, td, t) = [x[s == td[end] ? end : searchsortedlast(td, s)] for s in t]
+
+"""
+    viridis
+
+Map `s` onto the first `maxviri` viridis colors 
+"""
+function viridis(s, alpha = 0.9f0, maxviri = 200, RGBA=RGBA) 
+    l, u = extrema(s)
+    map(x->RGBA(Float32.(Bridge._viridis[1+floor(Int,x*maxviri)])..., alpha), (s .- l)/(u .- l))
+end
+
+
+function iterateall(it)
+    u = iterate(it)
+    u === nothing && error("Empty iterator")
+    x, state = u
+    while true
+        u = iterate(it, state)
+        u === nothing && return x
+        x, state = u
+    end
+end
+
 
 """
     _viridis

src/ode!.jl

@@ -29,4 +29,20 @@ function solvebackward!(::Bridge.R3!, f!, X, xT)
         yy[.., i] = y
     end
     X
+end
+
+function solvebackward!(::Bridge.R3!, f!, X::SamplePath, xT)
+    tt = X.tt
+    n = length(tt)
+    yy = X.yy
+    y = copy(xT)
+    z = yy[n] 
+    z .= y 
+    ws = (copy(y), copy(y), copy(y), copy(y)) # y2, k1, k2, k3
+    for i in n-1:-1:1
+        kernelr3!(f!, tt[i+1], y, ws, y, tt[i] - tt[i+1])  
+        z = yy[i] 
+        z .= y
+    end
+    X
 end

src/ode.jl

@@ -73,7 +73,7 @@ function kernelbs3(t, y, dt, P, k = F(t, y, P))
     yº, k4, err
 end
 
-function solvebackward!(method, F, X, xT, P) 
+function solvebackward!(method::ODESolver, F, X, xT, P) 
     _solvebackward!(method, F, X, xT, P) 
 end
 @inline function _solvebackward!(::R3, F, X::SamplePath{T}, xT, P) where {T}

src/poisson.jl

@@ -20,7 +20,7 @@ Sampling method for `InhomogPoisson` by the 'thinning' algorithm.
 #### Examples:
 
 ```
-sample(Thinning(λmax), T, InhomogPoisson(λ))
+sample(ThinningAlg(λmax), T, InhomogPoisson(λ))
 ```
 """
 struct ThinningAlg <: PoissonSampler

src/sde!.jl

@@ -7,7 +7,7 @@ struct BridgePre! <: SDESolver
 end
 
 
-function solve!(::EulerMaruyama!, Y, u::T, W::AbstractPath, P::ProcessOrCoefficients) where {T}
+function solve!(solver::Union{EulerMaruyama!,EulerMaruyamaWithIndex!}, Y, u::T, W::AbstractPath, P::ProcessOrCoefficients) where {T}
     N = length(W)
     N != length(Y) && error("Y and W differ in length.")
 
@@ -21,6 +21,7 @@ function solve!(::EulerMaruyama!, Y, u::T, W::AbstractPath, P::ProcessOrCoeffici
     tmp1 = copy(y)
     tmp2 = copy(y)
     dw = W.yy[.., 1]
+    dw2 = copy(dw)
     for i in 1:N-1
         t¯ = tt[i]
         dt = tt[i+1] - t¯ 
@@ -30,8 +31,13 @@ function solve!(::EulerMaruyama!, Y, u::T, W::AbstractPath, P::ProcessOrCoeffici
         for k in eachindex(dw)
             @inbounds dw[k] = W.yy[k, i+1] - W.yy[k, i]
         end
-        b!(t¯, y, tmp1, P)
-        σ!(t¯, y, dw, tmp2, P)
+        if solver isa EulerMaruyamaWithIndex!
+            bi(i, y, tmp1, P)
+            σi!(i, y, dw, tmp2, P)
+        else
+            b!(t¯, y, tmp1, P)
+            σ!(t¯, y, dw, tmp2, P)
+        end
         for k in eachindex(y)
             @inbounds y[k] = y[k] + tmp1[k]*dt + tmp2[k]
         end
@@ -41,6 +47,38 @@ function solve!(::EulerMaruyama!, Y, u::T, W::AbstractPath, P::ProcessOrCoeffici
 end
 
 
+function solve!(::EulerMaruyama!, Y, u::T, W::SamplePath, P::ProcessOrCoefficients) where {T}
+    N = length(W)
+    N != length(Y) && error("Y and W differ in length.")
+
+    tt = Y.tt
+    tt[:] = W.tt
+    yy = Y.yy
+    y::T = copy(u)
+
+    size(Y.yy) != (length(y), N) && error("Starting point has wrong length.")
+    #assert(size(W.yy) == (length(y), N))
+    tmp1 = copy(y)
+    tmp2 = copy(y)
+    dw = W.yy[.., 1]
+    for i in 1:N-1
+        t¯ = tt[i]
+        dt = tt[i+1] - t¯ 
+        for k in eachindex(tmp1)
+            @inbounds yy[k, i] = y[k]
+        end
+        for k in eachindex(dw)
+            @inbounds dw[k] = W.yy[k, i+1] - W.yy[k, i]
+        end
+        b!(t¯, y, tmp1, P)
+        σ!(t¯, y, dw, tmp2, P)
+        for k in eachindex(y)
+            @inbounds y[k] = y[k] + tmp1[k]*dt + tmp2[k]
+        end
+    end
+    yy[.., N] = y
+    Y
+end
 
 function bridge!(Y, u, W::VSamplePath, P::GuidedBridge!)
     W.tt === P.tt && error("Time axis mismatch between bridge P and driving W.") # not strictly an error

src/sde.jl

@@ -0,0 +1,42 @@
+#=
+function solve!(::EulerMaruyama, Y, u, W::SamplePath, P::ProcessOrCoefficients)
+    N = length(W)
+    N != length(Y) && error("Y and W differ in length.")
+
+    tt, yy = Y.tt, Y.yy
+    tt .= W.tt
+    ww = W.yy
+ 
+    y::typeof(u) = u
+
+    for i in 1:N-1
+        dt = tt[i+1] - tt[i]
+        yy[i] = y
+        y = y + b(tt[i], y, P)*dt + _scale((ww[i+1]-ww[i]), σ(tt[i], y,  P))
+    end
+    yy[N] = y
+    Y
+end
+=#
+
+"""
+Currently only timedependent sigma, as Ito correction is necessary
+"""
+function solvebackward!(::EulerMaruyama, Y, u, W::SamplePath, P::ProcessOrCoefficients)
+    N = length(W)
+    N != length(Y) && error("Y and W differ in length.")
+
+    tt, yy = Y.tt, Y.yy
+    tt .= W.tt
+    ww = W.yy
+
+    y::typeof(u) = u
+
+    for i in N:-1:2
+        dt = tt[i-1] - tt[i]
+        yy[i] = y
+        y = y + b(tt[i], y, P)*dt + _scale((ww[i-1]-ww[i]), σ(tt[i], P))
+    end
+    yy[1] = y
+    Y
+end

src/types.jl

@@ -1,3 +1,8 @@
+ismutable(_) = Val(false)
+ismutable(::Type{<:Array}) = Val(true)
+
+
+
 import Base: getindex, setindex!, length, copy, vcat, keys, values, iterate
 import Base: zero
 
@@ -66,6 +71,7 @@ end
 SamplePath(tt, yy::Vector{T}) where {T} = SamplePath{T}(tt, yy)
 
 samplepath(tt, v) = SamplePath(tt, fill(v, length(tt))) 
+samplepath(tt, v::Vector) = SamplePath(tt, [copy(v) for t in tt]) 
 
 
 copy(X::SamplePath{T}) where {T} = SamplePath{T}(copy(X.tt), copy(X.yy))