# GPLinearODEMaker.jl

GPLinearODEMaker (GLOM) is a package for finding the likelihood (and derivatives thereof) of multivariate Gaussian processes (GP) that are composed of a linear combination of a univariate GP and its derivatives.

where each X(t) is the latent GP and the qs are the time series of the outputs.

Here's an example using sine and cosines as the outputs to be modelled. The f, g!, and h! functions at the end give the likelihood, gradient, and Hessian, respectively.

```
import GPLinearODEMaker; GLOM = GPLinearODEMaker
kernel, n_kern_hyper = GLOM.include_kernel("se")
n = 100
xs = 20 .* sort(rand(n))
noise1 = 0.1 .* ones(n)
noise2 = 0.2 .* ones(n)
y1 = sin.(xs) .+ (noise1 .* randn(n))
y2 = cos.(xs) .+ (noise2 .* randn(n))
ys = collect(Iterators.flatten(zip(y1, y2)))
noise = collect(Iterators.flatten(zip(noise1, noise2)))
glo = GLOM.GLO(kernel, n_kern_hyper, 2, 2, xs, ys; noise = noise, a=[[1. 0.1];[0.1 1]])
total_hyperparameters = append!(collect(Iterators.flatten(glo.a)), [10])
workspace = GLOM.nlogL_matrix_workspace(glo, total_hyperparameters)
function f(non_zero_hyper::Vector{T} where T<:Real) = GLOM.nlogL_GLOM!(workspace, glo, non_zero_hyper) # feel free to add priors here to optimize on the posterior!
function g!(G::Vector{T}, non_zero_hyper::Vector{T}) where T<:Real
G[:] = GLOM.∇nlogL_GLOM!(workspace, glo, non_zero_hyper) # feel free to add priors here to optimize on the posterior!
end
function h!(H::Matrix{T}, non_zero_hyper::Vector{T}) where T<:Real
H[:, :] = GLOM.∇∇nlogL_GLOM!(workspace, glo, non_zero_hyper) # feel free to add priors here to optimize on the posterior!
end
```

You can use f, g!, and h! to optimize the GP hyperparameters with external packages like Optim.jl or Flux.jl

```
initial_x = GLOM.remove_zeros(total_hyperparameters)
using Optim
# @time result = optimize(f, initial_x, NelderMead()) # slow or wrong
# @time result = optimize(f, g!, initial_x, LBFGS()) # faster and usually right
@time result = optimize(f, g!, h!, initial_x, NewtonTrustRegion()) # fastest and usually right
fit_total_hyperparameters = GLOM.reconstruct_total_hyperparameters(glo, result.minimizer)
```

Once you have the best fit hyperparameters, you can easily calculate the GP conditioned on the data (i.e. the GP posterior)

```
n_samp_points = convert(Int64, max(500, round(2 * sqrt(2) * length(glo.x_obs))))
x_samp = collect(range(minimum(glo.x_obs); stop=maximum(glo.x_obs), length=n_samp_points))
n_total_samp_points = n_samp_points * glo.n_out
n_meas = length(glo.x_obs)
mean_GP, σ, mean_GP_obs, Σ = GLOM.GP_posteriors(glo, x_samp, fit_total_hyperparameters; return_mean_obs=true)
```

and use Plots to visualize the results

```
using Plots
function make_plot(output::Integer, label::String)
sample_output_indices = output:glo.n_out:n_total_samp_points
obs_output_indices = output:glo.n_out:length(ys)
p = scatter(xs, ys[obs_output_indices], yerror=noise1, label=label)
plot!(x_samp, mean_GP[sample_output_indices]; ribbon=σ[sample_output_indices], alpha=0.3, label="GP")
return p
end
plot(make_plot(1, "Sin"), make_plot(2, "Cos"), layout=(2,1), size=(960,540))
```

# Documentation

For more details and options, see the documentation

You can read about the first usage of this package in our paper

# Installation

The most current, tagged version of GPLinearODEMaker.jl can be easily installed using Julia's Pkg

`Pkg.add("GPLinearODEMaker")`

If you would like to contribute to the package, or just want to run the latest (untagged) version, you can use the following

`Pkg.develop("GPLinearODEMaker")`

# Citation

If you use `GPLinearODEMaker.jl`

in your work, please cite the BibTeX entry given in CITATION.bib

The formula images in this README created with this website