**BREAKING**: makes GP evaluation explicitly require RowVecs/ColVecs wrapper; calling GP with AbstractMatrix now errors

Project.toml

@@ -1,7 +1,7 @@
 name = "AbstractGPs"
 uuid = "99985d1d-32ba-4be9-9821-2ec096f28918"
 authors = ["JuliaGaussianProcesses Team"]
-version = "0.5.4"
+version = "0.6.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

docs/Project.toml

@@ -4,5 +4,5 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
-AbstractGPs = "0.4, 0.5"
+AbstractGPs = "0.4, 0.5, 0.6"
 Documenter = "0.27"

docs/src/api.md

@@ -48,6 +48,16 @@ fx = f(x, Σ)
 fx = AbstractGPs.FiniteGP(f, x, Σ)
 ```
 
+!!! note
+
+    When considering higher-dimensional input features, these are often stored
+    as a matrix `X`. Evaluating `f(X)` would be ambiguous, as it is not clear
+    whether to interpret rows or columns as feature vectors. Instead, use
+    `f(ColVecs(X))` or `f(RowVecs(X))` depending on your data layout. ColVecs
+    and RowVecs are lightweight wrappers provided by KernelFunctions.jl, and
+    allow for faster implementations.
+    For more details see [KernelFunctions.jl documentation on Input Types](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/api/#Vector-Valued-Inputs).
+
 The `FiniteGP` has two API levels.
 The [Primary Public API](@ref) should be supported by all `FiniteGP`s, while the [Secondary Public API](@ref) will only be supported by a subset.
 Use only the primary API when possible.

docs/src/concrete_features.md

@@ -1,26 +1,6 @@
 # Features
 
-## Plotting
-
-### Plots.jl
-
-We provide functions for plotting samples and predictions of Gaussian processes with [Plots.jl](https://github.com/JuliaPlots/Plots.jl). You can see some examples in the [One-dimensional regression](@ref) tutorial.
-
-```@docs
-AbstractGPs.RecipesBase.plot(::AbstractVector, ::AbstractGPs.FiniteGP)
-AbstractGPs.RecipesBase.plot(::AbstractGPs.FiniteGP)
-AbstractGPs.RecipesBase.plot(::AbstractVector, ::AbstractGPs.AbstractGP)
-sampleplot
-```
-
-### Makie.jl
-
-You can use the Julia package [AbstractGPsMakie.jl](https://github.com/JuliaGaussianProcesses/AbstractGPsMakie.jl) to plot Gaussian processes with [Makie.jl](https://github.com/JuliaPlots/Makie.jl).
-
-![posterior animation](https://juliagaussianprocesses.github.io/AbstractGPsMakie.jl/stable/posterior_animation.mp4)
-
-
-## Setup
+### Setup
 
 ```julia
 using AbstractGPs, Random
@@ -34,85 +14,108 @@ x = randn(rng, 10)
 y = randn(rng, 10)
 ```
 
-### Finite dimensional projection
+## Finite dimensional projection
 Look at the finite-dimensional projection of `f` at `x`, under zero-mean observation noise with variance `0.1`.
 ```julia
 fx = f(x, 0.1)
 ```
 
-### Sample from GP from the prior at `x` under noise.
+## Sample from GP from the prior at `x` under noise.
 ```julia
 y_sampled = rand(rng, fx)
 ```
 
-### Compute the log marginal probability of `y`.
+## Compute the log marginal probability of `y`.
 ```julia
 logpdf(fx, y)
 ```
 
-### Construct the posterior process implied by conditioning `f` at `x` on `y`.
+## Construct the posterior process implied by conditioning `f` at `x` on `y`.
 ```julia
 f_posterior = posterior(fx, y)
 ```
 
-### A posterior process follows the `AbstractGP` interface, so the same functions which work on the posterior as on the prior.
+## A posterior process follows the `AbstractGP` interface, so the same functions which work on the posterior as on the prior.
 
 ```julia
 rand(rng, f_posterior(x))
 logpdf(f_posterior(x), y)
 ```
 
-### Compute the VFE approximation to the log marginal probability of `y`.
-Here, z is a set of pseudo-points. 
+## Compute the VFE approximation to the log marginal probability of `y`.
+Here, z is a set of pseudo-points.
 ```julia
 z = randn(rng, 4)
 ```
 
-#### Evidence Lower BOund (ELBO)
-We provide a ready implentation of elbo w.r.t to the pseudo points. We can perform Variational Inference on pseudo-points by maximizing the ELBO term w.r.t pseudo-points `z` and any kernel parameters. For more information, see [examples](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl/tree/master/examples). 
+### Evidence Lower BOund (ELBO)
+We provide a ready implentation of elbo w.r.t to the pseudo points. We can perform Variational Inference on pseudo-points by maximizing the ELBO term w.r.t pseudo-points `z` and any kernel parameters. For more information, see [examples](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl/tree/master/examples).
 ```julia
 elbo(VFE(f(z)), fx, y)
 ```
 
-#### Construct the approximate posterior process implied by the VFE approximation.
+### Construct the approximate posterior process implied by the VFE approximation.
 The optimal pseudo-points obtained above can be used to create a approximate/sparse posterior. This can be used like a regular posterior in many cases.
 ```julia
 f_approx_posterior = posterior(VFE(f(z)), fx, y)
 ```
 
-#### An approximate posterior process is yet another `AbstractGP`, so you can do things with it like
+### An approximate posterior process is yet another `AbstractGP`, so you can do things with it like
 ```julia
 marginals(f_approx_posterior(x))
 ```
 
-### Sequential Conditioning 
+## Sequential Conditioning
 Sequential conditioning allows you to compute your posterior in an online fashion. We do this in an efficient manner by updating the cholesky factorisation of the covariance matrix and avoiding recomputing it from original covariance matrix.
 
 ```julia
 # Define GP prior
 f = GP(SqExponentialKernel())
 ```
 
-#### Exact Posterior
+### Exact Posterior
+
+Generate posterior with the first batch of data by conditioning the prior on them:
 ```julia
-# Generate posterior with the first batch of data on the prior f1.
 p_fx = posterior(f(x[1:3], 0.1), y[1:3])
+```
 
-# Generate posterior with the second batch of data considering posterior p_fx1 as the prior.
+Generate posterior with the second batch of data, considering the previous posterior `p_fx` as the prior:
+```julia
 p_p_fx = posterior(p_fx(x[4:10], 0.1), y[4:10])
 ```
 
-#### Approximate Posterior
-##### Adding observations in an sequential fashion
+### Approximate Posterior
+#### Adding observations in a sequential fashion
 ```julia
 Z1 = rand(rng, 4)
 Z2 = rand(rng, 3)
 Z = vcat(Z1, Z2)
 p_fx1 = posterior(VFE(f(Z)), f(x[1:7], 0.1), y[1:7])
 u_p_fx1 = update_posterior(p_fx1, f(x[8:10], 0.1), y[8:10])
 ```
-##### Adding pseudo-points in an sequential fashion
+
+#### Adding pseudo-points in a sequential fashion
 ```julia
 p_fx2 = posterior(VFE(f(Z1)), f(x, 0.1), y)
 u_p_fx2 = update_posterior(p_fx2, f(Z2))
 ```
+
+## Plotting
+
+### Plots.jl
+
+We provide functions for plotting samples and predictions of Gaussian processes with [Plots.jl](https://github.com/JuliaPlots/Plots.jl). You can see some examples in the [One-dimensional regression](@ref) tutorial.
+
+```@docs
+AbstractGPs.RecipesBase.plot(::AbstractVector, ::AbstractGPs.FiniteGP)
+AbstractGPs.RecipesBase.plot(::AbstractGPs.FiniteGP)
+AbstractGPs.RecipesBase.plot(::AbstractVector, ::AbstractGPs.AbstractGP)
+sampleplot
+```
+
+### Makie.jl
+
+You can use the Julia package [AbstractGPsMakie.jl](https://github.com/JuliaGaussianProcesses/AbstractGPsMakie.jl) to plot Gaussian processes with [Makie.jl](https://github.com/JuliaPlots/Makie.jl).
+
+![posterior animation](https://juliagaussianprocesses.github.io/AbstractGPsMakie.jl/stable/posterior_animation.mp4)

examples/intro-1d/Project.toml

@@ -13,7 +13,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 
 [compat]
-AbstractGPs = "0.4, 0.5"
+AbstractGPs = "0.4, 0.5, 0.6"
 AdvancedHMC = "0.2, 0.3"
 Distributions = "0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
 DynamicHMC = "2.2, 3.1"

examples/mauna-loa/Project.toml

@@ -9,7 +9,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
-AbstractGPs = "0.5"
+AbstractGPs = "0.5, 0.6"
 CSV = "0.9"
 DataFrames = "1"
 Literate = "2"

src/finite_gp_projection.jl

@@ -20,10 +20,19 @@ function FiniteGP(f::AbstractGP, x::AbstractVector, σ²::Real=default_σ²)
     return FiniteGP(f, x, Fill(σ², length(x)))
 end
 
-function FiniteGP(
-    f::AbstractGP, X::AbstractMatrix, σ²=default_σ²; obsdim::Int=KernelFunctions.defaultobs
-)
-    return FiniteGP(f, KernelFunctions.vec_of_vecs(X; obsdim=obsdim), σ²)
+function FiniteGP(f::AbstractGP, X::AbstractMatrix, σ²=default_σ²)
+    nrows, ncols = size(X)
+    throw(
+        ArgumentError(
+            """projecting a GP on a finite number of points only accepts vectors (e.g. a vector of vectors for multi-dimensional features), but was called with a $nrows × $ncols matrix.
+- If this is supposed to represent $ncols data points of dimension $nrows,
+  wrap it in ColVecs (e.g. `f(ColVecs(X), σ²)`).
+- If this is supposed to represent $nrows data points of dimension $ncols,
+  wrap it in RowVecs (e.g. `f(RowVecs(X), σ²)`).
+These are lightweight wrappers around the matrix, which allows the fastest performance. For more details, see https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/api/#Input-Types
+""",
+        ),
+    )
 end
 
 ## conversions

src/sparse_approximations.jl

@@ -281,7 +281,11 @@ function _compute_intermediates(fx::FiniteGP, y::AbstractVector{<:Real}, fz::Fin
 end
 
 function consistency_check(fx, y)
-    @assert length(fx) == length(y)
+    return length(fx) == length(y) || throw(
+        DimensionMismatch(
+            "length(fx) = $(length(fx)) and length(y) = $(length(y)) must match"
+        ),
+    )
 end
 
 function tr_Cf_invΣy(f::FiniteGP, Σy::Diagonal)

test/finite_gp_projection.jl

@@ -30,10 +30,6 @@ end
         f = GP(sin, SqExponentialKernel())
         fx, fx′ = FiniteGP(f, x, Σy), FiniteGP(f, x′, Σy′)
 
-        @test FiniteGP(f, Xmat; obsdim=1) == FiniteGP(f, RowVecs(Xmat))
-        @test FiniteGP(f, Xmat; obsdim=2) == FiniteGP(f, ColVecs(Xmat))
-        @test FiniteGP(f, Xmat, σ²; obsdim=1) == FiniteGP(f, RowVecs(Xmat), σ²)
-        @test FiniteGP(f, Xmat, σ²; obsdim=2) == FiniteGP(f, ColVecs(Xmat), σ²)
         @test mean(fx) == mean(f, x)
         @test cov(fx) == cov(f, x)
         @test var(fx) == diag(cov(fx))
@@ -232,6 +228,11 @@ end
             first(FiniteDifferences.grad(central_fdm(3, 1), Base.Fix1(logpdf, fx), y))
         @test Distributions.sqmahal!(r, fx, Y) ≈ Distributions.sqmahal(fx, Y)
     end
+    @testset "matrix inputs" begin
+        f = GP(SqExponentialKernel())
+        X = randn(5, 3)
+        @test_throws ArgumentError f(X)
+    end
 end
 
 @testset "Docs" begin