zero allocation

ode.jl

using Bridge, Distributions, StaticArrays
using Base.Test

"""
    kernelrk(f, t, y, dt, P)

Ralston (1965) update (order 3 step of the Bogacki–Shampine 1989 method)
to solve ``y(t + dt) - y(t) = \\int_t^{t+dt} f(s, y(s)) ds``.
"""
function kernelrk(f, t, y, dt, P)
    k1 = f(t, y, P)
    k2 = f(t + 1/2*dt, y + 1/2*dt*k1, P)
    k3 = f(t + 3/4*dt, y + 3/4*dt*k2, P)
    y + dt*(2/9*k1 + 1/3*k2 + 4/9*k3)
end

@inline _dK(t, X, P) = P.B*X + X*P.B' - P.a

"""
    gpK!(K::SamplePath, P)

Precompute ``K = H^{-1}`` from ``(d/dt)K = BK + KB' + a`` for a guided proposal.
"""
function gpK!(K::SamplePath{T}, P) where {T}
    tt = K.tt
    yy = K.yy
    yy[end] = y::T = zero(T)
    for i in length(tt)-1:-1:1
        y = kernelrk(_dK, tt[i+1], y, tt[i] - tt[i+1], P)  
        yy[i] = y  
    end
    K
end

n = 1000
TT = 5
tt = 0.0:TT/n:TT
m = 1000
S = SVector{2,Float64}
M = SArray{Tuple{2,2},Float64,2,4}
B = M([-1 0.1; -0.2 -1])
mu = 0*S([0.2, 0.3])
sigma = M(2*[-0.212887  0.0687025
  0.193157  0.388997 ])
a = sigma*sigma'
  
P = LinPro(B, mu, sigma)
type LP{S,T}
   B::S
   a::T 
end
P2 = LP(B, a)

t = 0.5
T = 2.0
n2 = 150
tt = linspace(t, T, n2)
K = SamplePath(tt, zeros(M, length(tt)))
gpK!(K, P2)
@test norm(K.yy[1]*Bridge.H(t, T, P)-I) < 10/n2^3

@test (@allocated gpK!(K, P2)) == 0