Merge branch 'refactor'

REQUIRE

@@ -4,5 +4,5 @@ Distributions
 Polynomials
 QuadGK
 ProgressMeter
-Klara
 ApproxFun
+StatsBase

src/PDSampler.jl

@@ -2,9 +2,8 @@ __precompile__(true)
 
 module PDSampler
 
-using Compat
-using Klara.ess
 using ApproxFun
+using StatsBase: autocov
 using Polynomials:
         Poly,
         roots,
@@ -18,13 +17,7 @@ using DataStructures:
 using ProgressMeter
 using Distributions: Beta
 
-
-const Int   = Int64
-const Float = Float64
-
 export
-    Int,
-    Float,
     # MODELS
     loglik,
     gradloglik,

src/geometry.jl

@@ -3,24 +3,24 @@ export
     Polygonal,
     nextboundary
 
-@compat abstract type Domain end
+abstract type Domain end
 
 """
     Unconstrained
 
 Signature of domains without boundaries.
 """
-immutable Unconstrained <: Domain end
+struct Unconstrained <: Domain end
 
 """
     Polygonal
 
 Encapsulate a domain with polygonal boundaries. The boundaries are defined with
 the normals to the facets and the intercepts of the corresponding hyperplanes.
 """
-immutable Polygonal <: Domain
-    normals::Matrix{Float}    # dim CxD where C=#constraints, D=#features
-    intercepts::Vector{Float} # dim C
+struct Polygonal <: Domain
+    normals::Matrix{Real} # dim CxD where C=#constraints, D=#features
+    intercepts::Vector{Real} # dim C
 end
 
 # -----------------------------------------------------------------------------
@@ -30,8 +30,10 @@ end
 
 Return (NaN, NaN) corresponding to the unconstrained case.
 """
-nextboundary{T<:Vector{Float}}(ud::Unconstrained, x::T, v::T
-                               )::Tuple{Float,Float} = (NaN, NaN)
+function nextboundary(ud::Unconstrained,
+                      x::Vector{<:Real}, v::Vector{<:Real})
+    return (NaN, NaN)
+end
 
 """
     nextboundary(pd::Polygonal, x, v)
@@ -40,17 +42,17 @@ Return the time of next boundary hit along the current trajectory when the
 domain is polygonal and return the normal of the corresponding boundary.
 The current point is `x` and the velocity `v`.
 """
-function nextboundary{T<:Vector{Float}}(pd::Polygonal,x::T, v::T
-                                         )::Tuple{Float,T}
+function nextboundary(pd::Polygonal,
+                      x::Vector{<:Real}, v::Vector{<:Real})
     # hitting time along trajectory (x+tv) for a boundary (normal,intercept) is
     # t = intercept/(normal dot v) - (x dot normal)/(normal dot v)
-    nsv  = pd.normals*v
-    hits = (pd.intercepts-pd.normals*x) ./ nsv
+    nsv  = pd.normals * v
+    hits = (pd.intercepts - pd.normals * x) ./ nsv
     # hard threshold times with nsv ~ 0 (near-parallel case),
     # remove negative times and times ~ 0 (for numerical stability)
     hits[ map(|, abs.(nsv) .< 1e-10, hits .< 1e-10) ] = Inf
     # get time of hit + index of corresponding boundary
     (t_hit, j) = findmin(hits)
     # return time of hit + normal vector to boundary
-    (t_hit, pd.normals[j,:])
+    (t_hit, pd.normals[j, :])
 end

src/ippsampler.jl

@@ -5,22 +5,22 @@ export
     nextevent_bps_q,
     nextevent_zz
 
-@compat abstract type IPPSamplingMethod end
-@compat abstract type Thinning <: IPPSamplingMethod end
+abstract type IPPSamplingMethod end
+abstract type Thinning <: IPPSamplingMethod end
 
 """
     LinearBound
 
 Indicating you have access to a global bound on the eigenvalues of the Hessian.
 """
-immutable LinearBound <: Thinning
-    xstar::Vector{Float}
-    gllstar::Vector{Float}  # grad log likelihood around xstar
-    b::Float                # N*L where L is Lipschitz constant
+struct LinearBound <: Thinning
+    xstar::Vector{Real}
+    gllstar::Vector{Real}  # grad log likelihood around xstar
+    b::Real                # N*L where L is Lipschitz constant
     a::Function
     function LinearBound(xstar, gllstar, b)
         new(xstar, gllstar, b,
-            (x,v) -> max(dot(-gllstar,v),0.) + norm(x-xstar)*b )
+            (x, v) -> max(dot(-gllstar, v), 0.) + norm(x-xstar) * b)
     end
 end
 
@@ -31,17 +31,16 @@ Object returned when calling a `nextevent_*` type function. The bouncing time
 `tau` as well as whether it should be accepted or rejected `dobounce` (set to
 true always for analytical sampling). `flipindex` is to support ZZ sampling.
 """
-immutable NextEvent
-    tau::Float         # candidate first arrival time
+struct NextEvent
+    tau::Real # candidate first arrival time
     dobounce::Function # do bounce (if accept reject step)
     flipindex::Int     # flipindex (if ZZ style)
     function NextEvent(tau, dobounce, flipindex)
-        new(tau,dobounce,flipindex)
+        new(tau, dobounce, flipindex)
     end
 end
-function NextEvent(tau; dobounce=(g,v)->true, flipindex=-1)::NextEvent
-    NextEvent(tau,dobounce,flipindex)
-end
+NextEvent(tau; dobounce=(g,v)->true, flipindex=-1) =
+    NextEvent(tau, dobounce, flipindex)
 
 # -----------------------------------------------------------------------------
 
@@ -52,40 +51,42 @@ Return a bouncing time and corresponding intensity in the specific case of a
 Gaussian density (for which the simulation of the first arrival time of the
 corresponding IPP can be done exactly).
 """
-function nextevent_bps{T<:Vector{Float}}(g::MvGaussian, x::T, v::T)::NextEvent
+function nextevent_bps(g::MvGaussian,
+                       x::Vector{<:Real}, v::Vector{<:Real})
     # precision(g)*x - precision(g)*mu(g) --> precmult, precmu in Gaussian.jl
-    a = dot(mvg_precmult(g,x)-mvg_precmu(g), v)
-    b = dot(mvg_precmult(g,v), v)
-    c = a/b
-    e = max(0.0,c)*c # so either 0 or c^2
-    tau = -c + sqrt( e + 2randexp()/b)
-    NextEvent(tau)
+    a = dot(mvg_precmult(g, x) - mvg_precmu(g), v)
+    b = dot(mvg_precmult(g, v), v)
+    c = a / b
+    e = max(0.0, c) * c # so either 0 or c^2
+    tau = -c + sqrt(e + 2randexp() / b)
+    return NextEvent(tau)
 end
 
-function nextevent_bps{T<:Vector{Float}}(g::PMFGaussian, x::T, w::T)::NextEvent
+function nextevent_bps(g::PMFGaussian,
+                       x::Vector{<:Real}, w::Vector{<:Real})
     # unmasking x, w
     xu, xv = x[1:g.d], x[g.d+1:end]
     wu, wv = w[1:g.d], w[g.d+1:end]
     # precomputing useful dot products
-    xuwv, xvwu = dot(xu,wv), dot(xv,wu)
-    wuwv, dxw  = dot(wu,wv), xuwv+xvwu
+    xuwv, xvwu = dot(xu, wv), dot(xv, wu)
+    wuwv, dxw = dot(wu, wv), xuwv + xvwu
     # real root
-    t0 = -0.5dxw/wuwv
+    t0 = -0.5dxw / wuwv
     # e(x) := ⟨x_u, x_v⟩ - r
-    ex = dot(xu,xv)-g.r
+    ex = dot(xu, xv) - g.r
     # (quadratic in t) e(x+tw) = e(x)+(⟨wu,xv⟩+⟨wv,xu⟩)t+⟨wu,wv⟩t^2
     p1 = Poly([ex, dxw, wuwv])
     # (linear in t) ⟨xv,wu⟩+⟨xu,wv⟩+2⟨wu,wv⟩t
     p2 = Poly([dxw, 2.0wuwv])
     # ⟨∇E(x+tw),w⟩ is a cubic in t (and the intensity)
-    p  = p1*p2
+    p  = p1 * p2
 
     rexp = randexp()
-    tau  = 0.0
+    tau = 0.0
 
     ### CASES (cf. models/pmf.jl)
 
-    Δ = t0^2-ex/wuwv # discriminant of p1
+    Δ = t0^2 - ex / wuwv # discriminant of p1
     if Δ <= 0
         # only single real root (t0)
         # two imaginary roots
@@ -96,7 +97,7 @@ function nextevent_bps{T<:Vector{Float}}(g::PMFGaussian, x::T, w::T)::NextEvent
         end
     else
         # three distinct real roots,
-        tm,tp = t0 + sqrt(Δ)*[-1.0,1.0]
+        tm, tp = t0 + sqrt(Δ) * [-1.0,1.0]
         if tp <= 0 # case: |/
             tau = pmf_caseA(rexp, p)
         elseif t0 <= 0 # case: |_/
@@ -107,28 +108,29 @@ function nextevent_bps{T<:Vector{Float}}(g::PMFGaussian, x::T, w::T)::NextEvent
             tau = pmf_caseD(rexp, p, tm, t0, tp)
         end
     end
-    NextEvent(tau)
+    return NextEvent(tau)
 end
 
 """
     nextevent_zz(g::MvGaussian, x, v)
 
 Same as nextevent but for the Zig Zag sampler.
 """
-function nextevent_zz{T<:Vector{Float}}(g::MvGaussian, x::T, v::T)::NextEvent
+function nextevent_zz(g::MvGaussian,
+                      x::Vector{<:Real}, v::Vector{<:Real})
     # # precision(g)*x - precision(g)*mu(g) --> precmult, precmu in Gaussian.jl
-    u1 = mvg_precmult(g,x)-mvg_precmu(g)
-    u2 = mvg_precmult(g,v)
+    u1 = mvg_precmult(g, x)-mvg_precmu(g)
+    u2 = mvg_precmult(g, v)
     taus = zeros(g.p)
     for i in 1:g.p
-        ai = u1[i]*v[i]
-        bi = u2[i]*v[i]
-        ci = ai./bi
-        ei = max(0.0,ci)*ci
-        taus[i] = -ci + sqrt(ei + 2.0randexp()/abs(bi))
+        ai = u1[i] * v[i]
+        bi = u2[i] * v[i]
+        ci = ai ./ bi
+        ei = max(0.0, ci) * ci
+        taus[i] = -ci + sqrt(ei + 2.0randexp() / abs(bi))
     end
     tau, flipindex = findmin(taus)
-    NextEvent(tau, flipindex=flipindex)
+    return NextEvent(tau, flipindex=flipindex)
 end
 
 """
@@ -137,25 +139,26 @@ end
 Return a bouncing time and corresponding intensity corresponding to a linear
 upperbound described in `lb`.
 """
-function nextevent_bps{T<:Vector{Float}}(lb::LinearBound, x::T,v::T)::NextEvent
+function nextevent_bps(lb::LinearBound,
+                       x::Vector{<:Real}, v::Vector{<:Real})
     a = lb.a(x, v)
     b = lb.b
-    @assert a>=0.0 && b>0.0 "<ippsampler/nextevent_bps/linearbound>"
-    tau       = -a/b + sqrt((a/b)^2 + 2randexp()/b)
-    lambdabar = a + b*tau
-    NextEvent(tau, dobounce=(g,v)->(rand()<-dot(g,v)/lambdabar) )
+    @assert a >= 0.0 && b > 0.0 "<ippsampler/nextevent_bps/linearbound>"
+    tau = -a / b + sqrt((a / b)^2 + 2randexp() / b)
+    lambdabar = a + b * tau
+    return NextEvent(tau, dobounce=(g,v)->(rand()<-dot(g, v)/lambdabar))
 end
 
-function nextevent_bps_q{T<:Vector{Float}}(gll::Function, x::T, v::T,
-                                            tref::Float; n=100)::NextEvent
+function nextevent_bps_q(gll::Function, x::Vector{<:Real}, v::Vector{<:Real},
+                         tref::Real; n=100)
 
-    chi(t) = max(0.0, dot(-gll(x+t*v),v))
+    chi(t) = max(0.0, dot(-gll(x + t * v), v))
     S      = Chebyshev(0.0..tref)
     p      = points(S, n)
     v      = chi.(p)
     f      = Fun(S, ApproxFun.transform(S,v))
     If     = cumsum(f) # integral from 0 to t with t < = tref
     tau    = ApproxFun.roots(If - randexp())
     tau    = (length(tau)>0) ? minimum(tau) : Inf
-    NextEvent(tau)
+    return NextEvent(tau)
 end

src/kernels.jl

@@ -15,9 +15,9 @@ export
 BPS specular reflection (in place) of a velocity `v` against a plane defined by
 its normal `n`.
 """
-function reflect_bps!{T<:Vector{Float}}(n::T, v::T)::T
-    v -= 2.0dot(n, v)*n/dot(n,n)
-    v
+function reflect_bps!(n::Vector{<:Real}, v::Vector{<:Real})
+    v .-= 2.0dot(n, v) * n / dot(n, n)
+    return v
 end
 
 """
@@ -26,10 +26,10 @@ end
 (Generalized) BPS specular reflection - flip normal component and resample
 orthogonal https://arxiv.org/abs/1706.04781
 """
-function reflect_gbps{T<:Vector{Float}}(n::T, v::T)::T
-    v2  = randn(length(n))
-    v2 -= dot(n, v2 + v)*n/dot(n,n)
-    v2
+function reflect_gbps(n::Vector{<:Real}, v::Vector{<:Real})
+    v2 = randn(length(n))
+    v2 .-= dot(n, v2 + v) * n / dot(n, n)
+    return v2
 end
 
 # """
@@ -38,7 +38,7 @@ end
 # (Generalized) BPS specular reflection flip normal component and resample orthogonal https://arxiv.org/abs/1706.04781
 #
 # """
-# function reflect_gbps_sphere!{T<:Vector{Float}}(n::T, v::T)::T
+# function reflect_gbps_sphere!{T<:Vector{Real}}(n::T, v::T)::T
 #     p = dot(n, v)*n/norm(n)^2
 #     v -= dot(n, v)*n/norm(n)^2
 #     v2 = rnorm(length(d))
@@ -52,39 +52,38 @@ end
 BPS specular reflection (in place) of a velocity `v` against a plane defined by
 its normal `n` and a mass matrix `mass`.
 """
-function reflect_bps!{T<:Vector{Float}}(n::T, v::T, mass::Matrix{Float})::T
-    v -= 2.0dot(n, v)*(mass*n)/dot(mass*n,n)
-    v
+function reflect_bps!(n::Vector{<:Real}, v::Vector{<:Real},
+                      mass::Matrix{<:Real})
+    v .-= 2.0dot(n, v) * (mass * n) / dot(mass * n, n)
+    return v
 end
 
-# ------------------------------------------------------------------------------
+# ----------------------------------------------------------------------------
 # ZZ KERNELS
 
 """
     reflect_zz!(mask, v)
 
 ZigZag reflection (in place) of v according to a mask on its indeces `mask`.
 """
-function reflect_zz!{T<:Vector{Float}}(mask::Vector{Int}, v::T)::T
-    v[mask] *= -1.0
-    v
-end
-reflect_zz!{T<:Vector{Float}}(mask::Int, v::T)::T = reflect_zz!([mask], v)
+reflect_zz!(mask::Vector{Int}, v::Vector{<:Real}) = (v[mask] .*= -1.0; v)
 
-# ------------------------------------------------------------------------------
+reflect_zz!(mask::Int, v::Vector{<:Real}) = reflect_zz!([mask], v)
+
+# --------------------------------------------------------------------------
 # Refreshment kernels
 
-function refresh_global!{T<:Vector{Float}}(v::T)::T
-    v = randn(length(v))
-    v
-end
-function refresh_restricted!{T<:Vector{Float}}(v::T)::T
+refresh_global!(v::Vector{<:Real}) = (v .= randn(length(v)); v)
+
+function refresh_restricted!(v::Vector{<:Real})
     v  = refresh_global!(v)
     v /= norm(v)
+    return v
 end
-function refresh_partial!{T<:Vector{Float}}(v::T, beta::Beta{Float})::T
+
+function refresh_partial!(v::Vector{<:Real}, beta::Beta{<:Real})
     # sample a vector with norm 1
-    w  = randn(length(v))
+    w = randn(length(v))
     w /= norm(w)
     # sample an angle
     angle = rand(beta) * 2 * pi
@@ -99,4 +98,4 @@ function refresh_partial!{T<:Vector{Float}}(v::T, beta::Beta{Float})::T
     # normalise
     v /= norm(v)
 end
-refresh_partial!(beta::Beta{Float}) = v->refresh_partial!(v,beta)
+refresh_partial!(beta::Beta{<:Real}) = v -> refresh_partial!(v, beta)