Refactoring in (almost) state-less form (#106)

* Formatting and more

* Stateless updates

* Small clean up

* Renamed AbstractGPs.jl

* Updates on the likelihood

* Compilable version

* Removed joinpath and others

* More changes

* More changes

* more changes duh

* gp classification passing!

* Last updates

* Formatting

* Modify multiclass

* Moar changes

* Online example working

* Actually fixed online GP

* Fixing the sampling

* Last tutorial working!

* More fixes

* Passing examples for regression

* Fix hyperparameters update online

* Fixed inference issues

* Bayesian SVM test passing

* Fixed the prior means

* Lots of likelihoods passing tests!

* Fixes Gaussian tests

* More fixes

* Fixed multiclass

* Fixes softmax

* All the fixes in da house

* Make them int

* Fix ambiguity

* Add some formatting

* Better typing

Author: theogf

Project.toml

@@ -1,15 +1,13 @@
 name = "AugmentedGaussianProcesses"
 uuid = "38eea1fd-7d7d-5162-9d08-f89d0f2e271e"
 authors = ["Theo Galy-Fajou <theo.galyfajou@gmail.com>"]
-version = "0.10.5"
+version = "0.11.0"
 
 [deps]
 AbstractMCMC = "80f14c24-f653-4e6a-9b94-39d6b0f70001"
-AdvancedHMC = "0bf59076-c3b1-5ca4-86bd-e02cd72cde3d"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
-Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InducingPoints = "b4bd816d-b975-4295-ac05-5f2992945579"
 KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
@@ -28,11 +26,9 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 AbstractMCMC = "2, 3"
-AdvancedHMC = "0.2.13, 0.3"
 ChainRulesCore = "0.9, 1"
 Distributions = "0.21.5, 0.22, 0.23, 0.24, 0.25"
 FastGaussQuadrature = "0.4"
-Flux = "0.10, 0.11, 0.12"
 ForwardDiff = "0.10"
 InducingPoints = "0.2"
 KernelFunctions = "0.8, 0.9, 0.10"

README.md

@@ -72,10 +72,10 @@ A complete documentation is available [in the docs](https://theogf.github.io/Aug
 ```julia
 using AugmentedGaussianProcesses;
 using KernelFunctions
-model = SVGP(X_train, Y_train, SqExponentialKernel(), LogisticLikelihood(),AnalyticSVI(100), 64)
-train!(model, 100)
-Y_predic = predict_y(model, X_test) #For getting the label directly
-Y_predic_prob, Y_predic_prob_var = proba_y(model,X_test) #For getting the likelihood (and likelihood uncertainty) of predicting class 1
+model = SVGP(SqExponentialKernel(), LogisticLikelihood(), AnalyticSVI(100), 64)
+train!(model, X_train, Y_train, 100)
+Y_predic = predict_y(model, X_test) # For getting the label directly
+Y_predic_prob, Y_predic_prob_var = proba_y(model, X_test) # For getting the likelihood (and likelihood uncertainty) of predicting class 1
 ```
 
 Both [documentation](https://theogf.github.io/AugmentedGaussianProcesses.jl/stable/) and [examples/tutorials](https://nbviewer.jupyter.org/github/theogf/AugmentedGaussianProcesses.jl/tree/master/examples/) are available.

docs/examples/gpclassification.jl

@@ -15,7 +15,7 @@ data = CSV.read(data.body, DataFrame)
 data.Class[data.Class .== 2] .= -1
 data = Matrix(data)
 X = data[:, 1:2]
-Y = data[:, end];
+Y = Int.(data[:, end]);
 
 # ### We create a function to visualize the data
 
@@ -28,21 +28,19 @@ plot_data(X, Y; size=(500, 500))
 
 # ### Run sparse classification with increasing number of inducing points
 Ms = [4, 8, 16, 32, 64]
-models = Vector{AbstractGP}(undef, length(Ms) + 1)
+models = Vector{AbstractGPModel}(undef, length(Ms) + 1)
 kernel = SqExponentialKernel() ∘ ScaleTransform(1.0)
 for (i, num_inducing) in enumerate(Ms)
     @info "Training with $(num_inducing) points"
-    m = SVGP(
-        X,
-        Y,
+    global m = SVGP(
         kernel,
         LogisticLikelihood(),
         AnalyticVI(),
-        num_inducing;
+        inducingpoints(KmeansAlg(num_inducing), X);
         optimiser=false,
         Zoptimiser=false,
     )
-    @time train!(m, 20)
+    @time train!(m, X, Y, 20)
     models[i] = m
 end
 # ### Running the full model

docs/examples/gpevents.jl

@@ -5,7 +5,6 @@
 # ### Loading necessary packages
 using Plots
 using AugmentedGaussianProcesses
-const AGP = AugmentedGaussianProcesses
 using Distributions
 
 # ## Creating some random data
@@ -58,17 +57,6 @@ function plot_model(model, X, Y, title=nothing)
         lab="",
         linewidth=3.0,
     )
-    if model isa SVGP # Plot the inducing points as well
-        Plots.plot!(
-            p,
-            vec(model.f[1].Z),
-            zeros(dim(model.f[1]));
-            msize=2.0,
-            color="black",
-            t=:scatter,
-            lab="",
-        )
-    end
     return p
 end;
 

docs/examples/gpregression.jl

@@ -2,7 +2,6 @@
 # 
 # ### Loading necessary packages
 using AugmentedGaussianProcesses
-const AGP = AugmentedGaussianProcesses
 using Distributions
 using Plots
 
@@ -27,23 +26,21 @@ scatter(X, Y; lab="")
 
 Ms = [4, 8, 16, 32, 64];
 # Create an empty array of GPs
-models = Vector{AbstractGP}(undef, length(Ms) + 1);
+models = Vector{AbstractGPModel}(undef, length(Ms) + 1);
 # Chose a kernel
 kernel = SqExponentialKernel();#  + PeriodicKernel()
 # And Run sparse classification with an increasing number of inducing points
 for (index, num_inducing) in enumerate(Ms)
     @info "Training with $(num_inducing) points"
     m = SVGP(
-        X,
-        Y, # First arguments are the input and output
         kernel, # Kernel
         GaussianLikelihood(σ), # Likelihood used
         AnalyticVI(), # Inference usede to solve the problem
-        num_inducing; # Number of inducing points used
+        inducingpoints(KmeansAlg(num_inducing), X); # Inducing points initialized with kmeans
         optimiser=false, # Keep kernel parameters fixed
         Zoptimiser=false, # Keep inducing points locations fixed
     )
-    @time train!(m, 100) # Train the model for 100 iterations
+    @time train!(m, X, Y, 100) # Train the model for 100 iterations
     models[index] = m # Save the model in the array
 end
 
@@ -123,7 +120,7 @@ Plots.plot(
 likelihoods = [
     StudentTLikelihood(3.0), LaplaceLikelihood(3.0), HeteroscedasticLikelihood(1.0)
 ]
-ngmodels = Vector{AbstractGP}(undef, length(likelihoods) + 1)
+ngmodels = Vector{AbstractGPModel}(undef, length(likelihoods) + 1)
 for (i, l) in enumerate(likelihoods)
     @info "Training with the $(l)" # We need to use VGP
     m = VGP(

docs/examples/multiclassgp.jl

@@ -39,7 +39,6 @@ end
 # ### And a function to plot the data
 function plot_data(X, Y, σ)
     p = Plots.plot(size(300, 500); lab="", title="sigma = $σ")
-    ys = unique(Y)
     Plots.scatter!(eachcol(X)...; group=Y, msw=0.0, lab="")
     return p
 end
@@ -48,21 +47,21 @@ plot([plot_data(generate_mixture_data(σ)..., σ) for σ in σs]...)
 
 # ## Model training
 # ### Run sparse multiclass classification with different level of noise
-models = Vector{AbstractGP}(undef, length(σs))
+models = Vector{AbstractGPModel}(undef, length(σs))
 kernel = SqExponentialKernel()
 num_inducing = 50
 for (i, σ) in enumerate(σs)
     @info "Training with data with noise $σ"
+    X, y = generate_mixture_data(σ)
     m = SVGP(
-        generate_mixture_data(σ)...,
         kernel,
         LogisticSoftMaxLikelihood(n_class),
         AnalyticVI(),
-        num_inducing;
+        inducingpoints(KmeansAlg(num_inducing), X);
         optimiser=false,
         Zoptimiser=false,
     )
-    @time train!(m, 20)
+    @time train!(m, X, y, 20)
     models[i] = m
 end
 

docs/examples/onlinegp.jl

@@ -35,10 +35,15 @@ IP_alg = OIPS(0.8);
 model = OnlineSVGP(k, GaussianLikelihood(σ), AnalyticVI(), IP_alg; optimiser=false)
 anim = Animation()
 size_batch = 100
-for (i, (X_batch, y_batch)) in enumerate(eachbatch((X_train, y_train); obsdim=1, size=size_batch))
-    train!(model, X_batch, y_batch; iterations=3)
-    plot_model(model, X, X_test, X_train[1:(i * size_batch)], y_train[1:(i * size_batch)])
-    frame(anim)
+let state = nothing
+    for (i, (X_batch, y_batch)) in
+        enumerate(eachbatch((X_train, y_train); obsdim=1, size=size_batch))
+        _, state = train!(model, X_batch, y_batch, state; iterations=5)
+        plot_model(
+            model, X, X_test, X_train[1:(i * size_batch)], y_train[1:(i * size_batch)]
+        )
+        frame(anim)
+    end
 end
 gif(anim; fps=4)
 # This works just as well with any likelihood! Just try it out!

docs/examples/sampling.jl

@@ -9,58 +9,41 @@ using AugmentedGaussianProcesses
 using Distributions
 using LinearAlgebra
 
-# ### Loading the banana dataset from OpenML
+# ### Generating some random data
 kernel = SqExponentialKernel()
-x = range(0, 10, length=50)
+x = range(0, 10; length=50)
 K = kernelmatrix(kernel, x)
 f = rand(MvNormal(K + 1e-8I)) # Sample a random GP
 y = rand.(Bernoulli.(AGP.logistic.(f)))
-y_sign = sign.(y .- 0.5)
-
+y_sign = Int.(sign.(y .- 0.5))
 
 # ### We create a function to visualize the data
 
-function plot_data(x, y; size=(300,500))
-    Plots.scatter(x,
-                y,
-                alpha=0.2,
-                markerstrokewidth=0.0,
-                lab="",
-                size=size
-            )
+function plot_data(x, y; size=(300, 500))
+    return Plots.scatter(x, y; alpha=0.2, markerstrokewidth=0.0, lab="", size=size)
 end
-plot_data(x, y; size = (500, 500))
+plot_data(x, y; size=(500, 500))
 
 # ### Run the variational gaussian process approximation
 @info "Running full model"
-mfull = VGP(x, y_sign,
-            kernel,
-            LogisticLikelihood(),
-            AnalyticVI(),
-            optimiser = false
-            )
+mfull = VGP(x, y_sign, kernel, LogisticLikelihood(), AnalyticVI(); optimiser=false)
 @time train!(mfull, 5)
 
 # ### We can also create a sampling based model
 @info "Sampling from model"
-mmcmc = MCGP(x, y,
-            kernel,
-            LogisticLikelihood(),
-            GibbsSampling(),
-            optimiser = false
-            )
+mmcmc = MCGP(x, y, kernel, LogisticLikelihood(), GibbsSampling(); optimiser=false)
 m = mmcmc
 @time samples = sample(mmcmc, 1000)
 
 # ### We can now visualize the results of both models
 
 # ### We first plot the latent function f (truth, the VI estimate, the samples)
-p1 = plot(x, f, label="true f")
-plot!(x, samples, label="", color=:black, alpha=0.02, lab="")
-plot!(x, mean(mfull[1]), ribbon=sqrt.(var(mfull[1])), label="VI")
+p1 = plot(x, f; label="true f")
+plot!(x, samples; label="", color=:black, alpha=0.02, lab="")
+plot!(x, mean(mfull[1]); ribbon=sqrt.(var(mfull[1])), label="VI")
 # ### And we can also plot the predictions vs the data
-p2 = plot_data(x, y; size=(600,400))
+p2 = plot_data(x, y; size=(600, 400))
 μ_vi, σ_vi = proba_y(mfull, x)
 plot!(x, μ_vi; ribbon=σ_vi, label="VI")
 μ_mcmc, σ_mcmc = proba_y(mmcmc, x)
-plot!(x, μ_mcmc; ribbon=σ_mcmc, label="MCMC")
+plot!(x, μ_mcmc; ribbon=σ_mcmc, label="MCMC")

docs/make.jl

@@ -28,23 +28,24 @@ end
 
 # Make the docs
 
-makedocs(modules = [AugmentedGaussianProcesses],
-         format = Documenter.Writers.HTMLWriter.HTML(
-            assets = ["assets/icon.ico"],
-            analytics = "UA-129106538-2",
-            ),
-         sitename= "AugmentedGaussianProcesses",
-         authors = "Théo Galy-Fajou",
-         pages = [
-         "Home" => "index.md",
-         "Background" => "background.md",
-         "User Guide" => "userguide.md",
-         "Kernels" => "kernel.md",
-         "Examples" => joinpath.("examples", filter(x -> endswith(x, ".md"), readdir(MD_OUTPUT))),
-         "Julia GP Packages" => "comparison.md",
-         "API" => "api.md"
-         ]
-         )
+makedocs(;
+    modules=[AugmentedGaussianProcesses],
+    format=Documenter.Writers.HTMLWriter.HTML(;
+        assets=["assets/icon.ico"], analytics="UA-129106538-2"
+    ),
+    sitename="AugmentedGaussianProcesses",
+    authors="Théo Galy-Fajou",
+    pages=[
+        "Home" => "index.md",
+        "Background" => "background.md",
+        "User Guide" => "userguide.md",
+        "Kernels" => "kernel.md",
+        "Examples" =>
+            joinpath.("examples", filter(x -> endswith(x, ".md"), readdir(MD_OUTPUT))),
+        "Julia GP Packages" => "comparison.md",
+        "API" => "api.md",
+    ],
+)
 
 # Deploy the docs
 

docs/src/examples/gpclassification.md

@@ -20,7 +20,7 @@ Run sparse classification with increasing number of inducing points
 
 ```@example gpclassification
 Ms = [4, 8, 16, 32, 64]
-models = Vector{AbstractGP}(undef, length(Ms) + 1)
+models = Vector{AbstractGPModel}(undef, length(Ms) + 1)
 kernel = transform(SqExponentialKernel(), 1.0)
 for (i, num_inducing) in enumerate(Ms)
     @info "Training with $(num_inducing) points"