Reorganise tests to match KernelFunctions (#78)

* Reorganise tests to match KernelFunctions

* Rearrange includes and remove elbo.jl

* Bump patch version

Project.toml

@@ -1,7 +1,7 @@
 name = "ApproximateGPs"
 uuid = "298c2ebc-0411-48ad-af38-99e88101b606"
 authors = ["JuliaGaussianProcesses Team"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918"

src/ApproximateGPs.jl

@@ -19,9 +19,8 @@ export SparseVariationalApproximation
 export DefaultQuadrature, Analytic, GaussHermite, MonteCarlo
 
 include("utils.jl")
-include("sparse_variational.jl")
 include("expected_loglik.jl")
-include("elbo.jl")
+include("sparse_variational.jl")
 
 using ForwardDiff
 

src/sparse_variational.jl

@@ -106,5 +106,79 @@ end
 
 inducing_points(f::ApproxPosteriorGP{<:SparseVariationalApproximation}) = f.approx.fz.x
 
-_chol_cov(q::AbstractMvNormal) = cholesky(Symmetric(cov(q)))
-_chol_cov(q::MvNormal) = cholesky(q.Σ)
+"""
+    elbo(svgp::SparseVariationalApproximation, fx::FiniteGP, y::AbstractVector{<:Real}; num_data=length(y), quadrature=DefaultQuadrature())
+
+Compute the Evidence Lower BOund from [1] for the process `f = fx.f ==
+svgp.fz.f` where `y` are observations of `fx`, pseudo-inputs are given by `z =
+svgp.fz.x` and `q(u)` is a variational distribution over inducing points `u =
+f(z)`.
+
+`quadrature` selects which method is used to calculate the expected loglikelihood in
+the ELBO. The options are: `DefaultQuadrature()`, `Analytic()`, `GaussHermite()` and
+`MonteCarlo()`. For likelihoods with an analytic solution, `DefaultQuadrature()` uses this
+exact solution. If there is no such solution, `DefaultQuadrature()` either uses
+`GaussHermite()` or `MonteCarlo()`, depending on the likelihood.
+
+N.B. the likelihood is assumed to be Gaussian with observation noise `fx.Σy`.
+Further, `fx.Σy` must be isotropic - i.e. `fx.Σy = α * I`.
+
+[1] - Hensman, James, Alexander Matthews, and Zoubin Ghahramani. "Scalable
+variational Gaussian process classification." Artificial Intelligence and
+Statistics. PMLR, 2015.
+"""
+function AbstractGPs.elbo(
+    sva::SparseVariationalApproximation,
+    fx::FiniteGP{<:AbstractGP,<:AbstractVector,<:Diagonal{<:Real,<:Fill}},
+    y::AbstractVector{<:Real};
+    num_data=length(y),
+    quadrature=DefaultQuadrature(),
+)
+    @assert sva.fz.f === fx.f
+    return _elbo(quadrature, sva, fx, y, GaussianLikelihood(fx.Σy[1]), num_data)
+end
+
+function AbstractGPs.elbo(
+    ::SparseVariationalApproximation, ::FiniteGP, ::AbstractVector; kwargs...
+)
+    return error(
+        "The observation noise fx.Σy must be homoscedastic.\n To avoid this error, construct fx using: f = GP(kernel); fx = f(x, σ²)",
+    )
+end
+
+"""
+    elbo(svgp, ::SparseVariationalApproximation, lfx::LatentFiniteGP, y::AbstractVector; num_data=length(y), quadrature=DefaultQuadrature())
+
+Compute the ELBO for a LatentGP with a possibly non-conjugate likelihood.
+"""
+function AbstractGPs.elbo(
+    sva::SparseVariationalApproximation,
+    lfx::LatentFiniteGP,
+    y::AbstractVector;
+    num_data=length(y),
+    quadrature=DefaultQuadrature(),
+)
+    @assert sva.fz.f === lfx.fx.f
+    return _elbo(quadrature, sva, lfx.fx, y, lfx.lik, num_data)
+end
+
+# Compute the common elements of the ELBO
+function _elbo(
+    quadrature::QuadratureMethod,
+    sva::SparseVariationalApproximation,
+    fx::FiniteGP,
+    y::AbstractVector,
+    lik,
+    num_data::Integer,
+)
+    @assert sva.fz.f === fx.f
+    post = posterior(sva)
+    q_f = marginals(post(fx.x))
+    variational_exp = expected_loglik(quadrature, y, q_f, lik)
+
+    kl_term = KL(sva.q, sva.fz)
+
+    n_batch = length(y)
+    scale = num_data / n_batch
+    return sum(variational_exp) * scale - kl_term
+end

src/utils.jl

@@ -8,3 +8,6 @@ if VERSION < v"1.7"
         return Cholesky{T,typeof(U.data)}(U.data, 'U', 0)
     end
 end
+
+_chol_cov(q::AbstractMvNormal) = cholesky(Symmetric(cov(q)))
+_chol_cov(q::MvNormal) = cholesky(q.Σ)