GNU Octave GP Emulation Package Manual

This manual is for the GNU Octave GP Emulation package, version 0.0.1. The sources of this manual are derived from Oliver Heimlich’s Interval package manual.

Copyright © 2016–2017 Juan Pablo Carbajal

Permission is granted to copy, distribute and/or modify this document under a Creative Commons Attribution 4.0 International License. A copy of the license is included in Creative Commons License.

Permission is granted to copy, distribute and/or modify the source files used for this document under the terms of GNU General Public License, Version 3 or any later version published by the Free Software Foundation. A copy of the license is included in GNU General Public License.

Next: , Up: (dir)   [Contents]

Table of Contents

Next: , Previous: , Up: Top   [Contents]


Welcome to the user manual of the GP Emulation package for GNU Octave. This chapter presents background information and may safely be skipped. First-time users who want to cut right to the chase should read Getting Started, which teaches basic concepts and first steps with the package.

Next: , Previous: , Up: Top   [Contents]

1 Getting Started

The following chapter take you by the hand and gives a quick overview on the GP emulation package’s basic capabilities. More detailed information can usually be found in the functions’ documentation.

1.1 Installation

Unless the package is provided by a third-party distributor, the interval package can be installed with the pkg command from within Octave. Assuming you have got a file named gpemulation-<version>.tar.gz, where <version> is a three digits version (here I use 1.0.0), and that you have this file in the folder returned by Octave’s pwd command; the installation proceeds as follows:

pkg install gpemulation-1.0.0.tar.gz
  -| For information about changes from previous versions
  -| of the gpemulation package, run 'news gpemulation'.

During this kind of installation parts of the gpemulation package are compiled for the target system, which requires development libraries for GNU Octave (version ≥ 3.8.0) to be installed. The official Octave release for Microsoft Windows already includes these dependencies. For other systems it might be necessary to install the package “liboctave-dev”, which are provided by most GNU distributions (names may vary).

In order to use the gpemulation package during an Octave session, it must have been loaded, i. e., added to the path. In the following parts of the manual it is assumed that the package has been loaded, which can be accomplished with the pkg load gpemulation command.

That’s it. After only two commands the package is ready to be used within Octave.

Next: , Previous: , Up: Top   [Contents]

2 Quick Introduction to Gaussian Processes

In this chapter I introduce Gaussian Processes succinctly and conceptually. For a more formal and detailed introduction you could refer to Guassian Processes for Machine Learning.

A stochastic process can be thought of the generalization of a probability distribution of a finite-dimensional random variable to functions. That is, the process draws a function from a distribution of functions. If the function values at each point in the functions domain follow a Gaussian distribution, the process is a Gaussian Process (GP). The figure below shows 500 draws from a GP. The vertical lines mark points at which the histograms are calculated. These histograms are well described with a Gaussian curve (shown in red) and this is true for any point in the domain.

An example of a Gaussian Process

If a finite-dimensional random variable follows a Gaussian distribution, this distribution is fully determined by the mean value and the variance. In the case of GPs, the distribution of functions is fully determined by a mean function and a covariance function.

The mean function m(x), is any function defined in the domain (i.e. the inputs) of the GP. The covariance function k(x,y), often called kernel, is a function defined for pairs of inputs and it must have certain properties:

A commonly used kernel function is the squared exponential

k(x,y) = exp [ - ( x - y)² / (2 s²) ]

Let’s see how we can sample functions from the GP defined by this kernel. We first need to generate a set of input values and build the covariance matrix. Then we need to sample a Gaussian distribution with this covariance. To do this last step we can use the function mvnrnd provided by the

x = linspace (0, 1, 100).';
K = exp ( -(x - x.').^2 / 2 / 0.025 );

# Lets take 10 samples from this GP and plot them.
pkg load statistics;
y = mvnrnd (0, K, 10).';
plot (x, y, "linewidth", 2);
axis tight
xlabel ("x");
ylabel ("y");
Plotting an interval function and a classic function

2.1 Calculating Covariance Functions

The covariance function of a distribution of functions is calculated as the expectation value of squared deviations from the mean:

k(x,x’) = cov[ y(x),y(x’) ] = E [ ( y(x) - m(x) ) ⋅ ( y(x’) - m(x’) )ᵀ ]

The expectation is calculated over realization of whatever random parameter defines the distribution. For example, take a distribution of parabolas with random concavities which follow a Gaussian distribution with mean value ã and variance v

y(x) = a ( x - x₀ )² 

a ~ N(ã,v)

The mean function is

m(x) = E[ y(x) ] = ∫ a p(a) ( x - x₀ )² da = ã ( x - x₀ )²

and the covariance function is then

k(x, x’) = ( x - x₀ )²( x’ - x₀ )² v
x = linspace (0, 1, 100).';          ## Inputs

a = -1;                              ## Mean concavity
v = 0.5^2;                           ## Variance of concavity
K =  (x-0.5).^2 .* (x.'-0.5).^2 * v; ## Covariance matrix
m = a * (x.' - 0.5).^2;              ## Mean vector

# Lets take 10 samples from this GP and plot them.
pkg load statistics;
y = mvnrnd (m, K, 10).';
plot (x, y, "linewidth", 2)
axis tight
xlabel ("x");
ylabel ("y");
Parabolas with Gaussian distributed concavities with mean value -1 and variance 0.025

This is how you build GPs from known function distributions. When you do not know the analytic distribution of the functions you can try to estimate it from a large number of samples trajectories. We do this in the example below, where the covariance is obtained from samples of the solutions of a stochastic differential equation. The result is compared with the analytical covariance of this SDE implemented in the function covLTIo1d1. See Covariance Function of a First Order Ordinary Differential Equation with Additive Gaussian Noise

# Comparing covariance calculated from n samples of a
# stochastic differential equation and the theoretical calculation.
nT = 1e3;                         # Time samples for Ito integration
t_ = linspace (0, 1, nT);         # Time vector for Ito integral
dt = t_(2) - t_(1);               # Time step
a = -3;
eA = exp (a * dt); # Evolution matrix

# Ito integral
s       = 1;
n       = 1e3;
dW      = sqrt (s*dt) * randn (nT, n);   # Wiener differential
y_      = zeros (nT,n);
y_(1,:) = 1;
for i = 1:nT-1
 y_(i+1,:) = y_(i,:) * eA + dW(i+1,:);

# Experimental covariance
i0 = nT/2;                              # index at t0 = 0.5
ym =  mean (y_, 2);
dy = y_ - ym;
k  = (dy * dy.'(:,i0)) / n;

# Theoretical covariance
t0 = 0.5;                               # centered in the interval
t  = linspace (t_(1),t_(end), 1e2).';
K  = covLTIo1d1 (t, t0, a, s, 0);

h = plot (t_, k,'-r;Experimental;', t, K,'-k;Exact;');
set (h,'linewidth', 2);
xlabel ('Time');
ylabel ('Covariance')
axis tight

Parabolas with Gaussian distributed concavities with mean value -1 and variance 0.025

2.2 Covariance Function of a First Order Ordinary Differential Equation with Additive Gaussian Noise

A first order linear ordinary differential equation (LODE) with additive Gaussian noise is also known as a linear stochastic differential equation (LSDE)

(t) = A y(t) + u(t) + ξ(t) 

y(0) = y₀ + ξξ(t) ~ N(0,𝛴) , ξ₀ ~ N(0,𝛴₀)

where u(t) is a deterministic function called actuation.

If there were no noise terms and no actuation, this equation has a unique solution (a.k.a. homogeneous solution)

yₕ(t) = exp(At)y

The effect of the actuation can be calculated using the Green’s functions of the LODE. This give raise to the particular solution of the LODE

yₚ(t) = ∫ G(t,r)u(r) dr 

G(t,r) = H(t-r)exp(A(t-r))

Where H(t-r) is the Heaviside function. Note that the particular solution is a linear transformation of the actuation and we can write it as yₚ = Gu. The general solution of the LODE can be written in the form

yₘ(t) = yₕ(t) + yₚ(t)

The noise terms in the LSDE induce a distribution of functions whose mean value coincides with the deterministic solution given above.

To calculate the covariance function of the distribution of solutions we need the concept of the Green’s function of the adjoint differenital operator, which is a big name for something simple.

Gᵀ(t,r) = H(r-t)exp(Aᵀ(r-t))

This will be needed when we calculate the covariance function (See covariance function's formula).

k(t,r) = cov[ y(t),y(r) ] = exp(At) 𝛴₀ exp(Ar) + G 𝛴 G

When the last term is expressed as an integral1 we get

k(t,r) = exp(At) 𝛴₀ exp(Ar) + ∫ G(t,z) 𝛴 Gᵀ(r,z) dz

The function covLTIo1d1 calculates this covariance function for one-dimensional LSDE, i.e. A ∈ ℝ

nT = 100;
t  = linspace (0, 1, nT).';   ## Inputs

a  = -5;                      ## Decay rate
s  = 0.3^2;                   ## Variance of noise
s0 = 0.25^2;                  ## Variance of initial conditions
K  = covLTIo1d1 (t,t,a,s,s0); ## Covariance matrix
m  = [1; -1].* exp (a*t.');   ## Mean functions

subplot (2,2,1)
colormap (cubehelix);
imagesc (t,t,K);

subplot (2,2,2)
plot (t,K(:,50));
axis tight

subplot (2,2,3:4)
# Lets take 10 samples from this GP and plot them.
pkg load statistics;
y          = zeros (nT,10);
y(:,1:5)   = mvnrnd (m(1,:), K, 5).';
y(:,6:end) = mvnrnd (m(2,:), K, 5).';
plot (t, y, "linewidth", 2)
axis tight
xlabel ("t");
ylabel ("y");
covariance functon of a 1D LSDE

The function covLTIo1diag calculates a covariance function for d-dimensional LSDE, i.e. A ∈ ℝᵈˣᵈ, when the system matrix can be diagonalized, i.e. A = VDV⁻¹ where D is a diagonal matrix.

nT = 100;
t = linspace (0,1,nT).';

# Oscillator
w     = 2*pi*5;
M     = [0 1; -(w)^2 -0.2*2*w];
[V d] = eig (M);
d     = diag (d);
Vi    = inv (V);

S  = diag([0 w]);                    # Covariance of noise
S0 = zeros(2);                       # Covariance of inital conditions
K  = covLTIo1diag ({t,{V,d,Vi}}, [], S, S0);
K  = covmat2cell (K, [length(t) length(d)]);

subplot (2,4,[1 2 5 6])
colormap (cubehelix);
imagesc (cell2mat (K));
hold on
line(nT*[0 1; 2 1], nT*[1 0; 1 2]);

ax = [3 7 4 8];
for i=1:4
  subplot (2,4,ax(i));
  plot(t,K{i}(:,nT/2),"linewidth", 2);
  [k,l] = ind2sub([2 2],i);
  title (sprintf("Component %d with %d",k,l));

covariance functon of a 1D LSDE

This function can also be used to calculate the covariance between two d-dimensional LSDE coupled by noise.

nT = 100;
t = linspace (0,1,nT).';

# Oscillator
w     = 2*pi*5;
M     = [0 1; -(w)^2 -0.1*2*w];
[V d] = eig (M);
Vi    = inv (V);
d     = diag (d);

# First order system
a = -3;

# combined eigenspace
V  = blkdiag (1,V);
Vi = blkdiag (1,Vi);
d  = [a;d];

# coupling
S  = ones(3);                        # Covariance of noise
S0 = zeros(3);                       # Covariance of inital conditions
K  = covLTIo1diag ({t,{V,d,Vi}}, [], S, S0);
K  = covmat2cell (K, [length(t) length(d)]);

subplot (3,6,[1:3 7:9 13:15])
colormap (cubehelix);
imagesc (cell2mat (K));
hold on
line(nT*[0 1; 3 1], nT*[1 0; 1 3]);

ax = [4 10 16 5 11 17 6 12 18];
for i=1:9
  subplot (3,6,ax(i));
  [k,l] = ind2sub([3 3],i);
  y = K{i}(:,nT/2);
  plot(t,y,"linewidth", 2);
  title (sprintf("%d-%d",k,l));

covariance functon of a coupled LSDE

2.3 Conditioning on data or What’s the Use of the Covariance Function?

Given a GP and a collection of input-output data (X,Y), i.e. the design data, we can calculate another GP that is consistent with this given data. This is done following Bayes’ theorem for conditional distributions. The conditioned GP’s mean function is used for prediction and the formula reads

y(x) = k(x,X) ( k(X,X) + 𝛌 )⁻¹ ( Y(X) - m(X) ) + m(x)

The covariance function evaluated at the design data inputs k(X,X) is the covariance matrix of the observations. That is, GPs allow the use of linear algebra in Bayesian inference, considerably simplifying computations. This is a important reason for the popularity of GPs.

The predictive mean, the function that we use to predict values, is a linear combination of columns of the covariance matrix k(x,X):

y(x) = k(x,X) 𝛂 + m(x) 

𝛂 = ( k(X,X) + 𝛌 )⁻¹ ( Y(X) - m(X) )

These columns were plotted in all the figures shown before. For example, the columns of the covariance of a first order system (See this figure) is a two sided decaying exponential centered at the data. For a second oder system the columns show oscillatory behavior (See this figure). It s useful when deciding what covariance to use to look at the data and decide what functions would be good for prediction.

The parameter 𝛌 called regularization parameter, it is free and we need to decide its value. It represents the trust we have on the mean function of our GP. If 𝛌 is big, the conditioned GP will pay little attention to the design data and mostly use the mean function for prediction. If 𝛌 is small, the conditioned GP will stick to the design data. The extreme case of 𝛌 = 0 forces the prediction to go exactly through the design data, this is called interpolation. The value of the regularization parameter can be selected based on prediction errors (like cross-validation) or it can be decided on knowledge about the data.

To keep things simple, lets use our stochastic parabolas (See the derivation here) for prediction. First let say we got this data consisting of 10 points (X,Y)

data used for parabolic prediction

Clearly whatever process generated this data, it cannot be described with a parabola. Indeed the data was generated with the function Y = 0.5*(X - 0.5).^2 - 5*(X - 0.5).^4

We build the covariance matrix by evaluating our covariance function at these points.

covariance for parabolic prediction

The figure also shows the columns of the covariance matrix k(x,X), with x ∈ [0,1]. This are the functions that will be used for prediction. As you can see our covariance function has columns that are parabolas of different concavity (which could be guesses by the way we built this covariance matrix). As we know the prediction will be a linear combination of these columns (See this equation). What do we get if we linearly combine parabolas of different concavity?

z(x) = ∑ aₙ (x - x₀)² = A (x - x₀)²

Well just another parabola. This means that our GP can only interpolate parabolas, in other words, the covariance matrix essentially generates one type of function. When this happens, we will always get covariance matrices with rank 1, because the columns are proportional to each other. Covariances with finite ranks are called degenerate covariance functions (or degenerate kernels). This also means that when the rank of the covariance is lower than the number of inputs, we will not be able to set 𝛌 = 0 when conditioning (See this equation).

Whenever you are given a covariance matrix, check its eigenvalue or singular value spectrum to know how powerful it is. In this case

s = svd (KXX).';
round (s * 1e2) * 1e-2
⇒0.05  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00

illustrating the degeneracy of our kernel.

We apply the conditioning formula (See here) using two different mean functions to show the effect they have on the solution. The first mean function we use is

mA(x) = -5 (x - 0.5)⁴

which will give us parabolic residuals and a good solution to the interpolation problem. The second mean functions to be used is

mB(x) = (x - 0.5)²

which is ok, but it does not help our covariance function very much, as the residuals will still be a mix of quartic and parabolic functions. We use this last function to show the effect of two extreme values of the regularization parameter. An solution without regularization 𝛌 = 1e-8 and a strongly regularized solution 𝛌 = 1.

The figure below shows the results.

covariance for parabolic prediction

As you can see the solution using mA (labeled A) is able to interpolate (𝛌 = 1e-8) because mA takes care of the quartic term and leaves parabolic residuals for the covariance, which is able to find the right concavity.

The second mean, mB, is not so useful. The regularized solution (labeled B reg) just sticks to the mean function with little correction. The solution without regularization (labeled B), finds the best parabola it can build. Indeed, the theory shows that this is the same result as doing a least square error regression using a parabola that goes through the origin. This is shown in the plot with the legend "fit".

Next: , Previous: , Up: Top   [Contents]

3 Introduction to Gaussian Processes based emulation

At this point you have a conceptual grasping of what a GP is and how they are built. In particular you know that there is a GP associated to a LSDE (a LODE with Gaussian noise actuation). In this chapter you will learn about GP emulators and how to build one to emulate a deterministic nonlinear system, i.e. the simulator.

The objective of an emulator is to reduce the computational cost needed to obtain a trajectory of a nonlinear complicated simulator, without losing too much accuracy.

To build an emulator you need samples from the simulator and the parameters you used to generate them. You will also need to make some decisions about the structure of the emulator. We discuss this in the subsequent sections.

The emulator will end up being a GP associated to a LSDE, similar to the ones we saw before (See Covariance Function of a First Order Ordinary Differential Equation with Additive Gaussian Noise). The particular way in which we construct this LSDE is what makes the associated GP an emulator of a nonlinear simulator.

The input to the emulator will be the time stamps (as before) and a parameter vector. That is, the emulator will make predictions of the simulator output at a given time and parameter value.

3.1 The Anatomy of an Emulator

3.1.1 Design datasets

Lets say you have a collection of n sample simulations from your simulator, i.e. the nonlinear model to emulate. Each simulation was evaluated at nT time stamps and the signal are d dimensional (that is you observe d signals from the simulator). We call this collection of sample trajectories the design data for the emulator, and we will organize them in a multidimensional array of size [nT n d].

The parameters used in the simulator to generate the n datasets are stored in a matrix of size [n p]. The i-th row of this matrix contains the parameters used to generate the corresponding dataset.

In the following example we generate design data with n = 80 samples from a nonlinear ODE (the Van der Pol oscillator)

 = y 

 = μ ( 1 -  x² ) y - x

That has one parameter (i.e. p = 1) controlling the shape of the oscillations of the system. Note that we remove the transient of the system and take only the data at the limit cycle. Emulating transient behavior is difficult so we skip it for now (See Advanced Examples)

nT  = 300;                          # Number of time samples
n   = 10;                           # Number of design datasets

T   = 30;                           # Time to wait before observation
t   = linspace (0, 18, nT).';       # Time samples
mu  = linspace (0.3, 4, n)';        # The oscillator parameters
d   = 2;                            # Dimension of the simulator signals

vdp = @(y,t, mu) [y(2); ...                 # Simulator is
          mu * (1 - y(1)^2) * y(2) - y(1)];  # Van der Pol oscillator

y0        = [0.5 0];                         # Initialization conditions
simulator = @(t,p) lsode (@(x,t)vdp(x,t,p), ...
        lsode(@(x,t)vdp(x,t,p), y0, [0 T])(end,:), ... # Initial condition

y    = zeros (nT, n, d);      # Matrix to store the design datasets
for i=1:n
  tmp = lsode (@(x,t)vdp(x,t,mu(i)), y0, [0 T]);
  tmp = tmp(end,:);           # Use final position as initial condition
  y(:,i,:) = lsode (@(x,t)vdp(x,t,mu(i)), tmp, t);

p   = 1;          # Dimension of simulator parameter vector
nlp = mu;         # Simulator parameter matrix
Samples from a Van der Pol oscillator

The figure shows the phase space with some of the trajectories on the limit cycle of the Van der Pol oscillator. To produce the figure we used a lot more time stamps than in the actual design data (nT = 40) and it only shows some of the n = 80 datasets.

For different parameter values, the shape of the cycle changes, the ellipse observed at low parameter values (system behaves more linearly) gets deformed as the parameter increases (nonlinear regime). The lateral plots show the two components of the data.

3.1.2 Prior Model

The way we calculate the covariance matrix is valid only for linear systems. Therefore, to build the emulator we propose a linear model. This linear model is built as an aggregation of independent linear submodels. Each of these submodels is called replica (when all submodels have the same structure) or proxy. We build a proxy for each dataset and as before we will organize the parameters of the proxy in a matrix of size [n q] (Note that number of parameters can change, i.e. q can be different from p).

Since we know the parameters of the dataset (this is part of the dataset as explained before) we need to define a way to obtain the parameters of the proxy. There is no fixed rule to do this, we can define and ad-hoc formula or we can fit our linear model to the i-th dataset. It all depends on how much information about the simulator you can use.

When we use the emulator to generate predictions at a previously unseen parameter vector, we will need the parameter vector for the corresponding (unseen) proxy. This means that we will need a mapping from the simulator parameters to the proxy parameters. If you used an ad-hoc formula to do this mapping, then you will use this formula for this as well. If you have fitted the datasets you will need to build an interpolant, which can be as simple as a linear interpolation or as complex as a GP.

In the following example we fit the design data See previous example and create a simulator-emulator parameter mapping with a linear interpolation. The proxies are all of the form

y(t) = A cos ( w t ) + (B/w) sin ( w t )

Which is the solution of the harmonic oscillator equation (second order LODE) with frequency w, initial deformation A and initial speed B. It will be important later to know that these functions are the first component of the solutions of the dynamical system

 = z 

 = -w² y

with y(0) = A and (0) = z(0) = B.

Since we use the same proxy structure for all datasets we could call them replicas. Once we have estimated the values of A,B,w for each data set generated with paramter μ, we generate a linear interpolation to map from any μ to the corresponding A,B,w.

##  Proxy calibration  ##
proxy = @(t,p) p(2) * cos (p(1) * t) + p(3)/p(1) * sin (p(1) * t);
q     = 3;                                  # Parameter vector size

z     = zeros (nT, n, 1);                   # Matrix to store the proxies
Q     = zeros (n, q);                       # Matrix of proxy parameters

tmp2 = [1 0 0.5];
for i=1:n
  tmp2   = sqp (tmp2, @(x)sumsq(y(:,i,1)-proxy(t,x)));
  Q(i,:) = tmp2;
  z(:,i) = proxy (t, tmp2);

##  Parameter mapping  ##
pp     = interp1 (mu, Q, "pp");
parmap = @(x) ppval(pp, x);
Proxies for samples from a Van der Pol oscillator

The top plot shows some design data and the fitted proxies. The bottom plot shows several trajectories that are not in the design data and the proxies we obtain with the parameter mapping.

The aggregated linear model built from the proxies can be used for the mean function of the emulator.

Note that so far the proxies are independent, while the datasets might not be. Therefore, we need to couple these models and we are going to use the covariance function to do this (See the equation).

3.1.3 Parameter Coupling

In an emulator the independent proxies, are coupled through a noisy actuation term. If two proxies are highly coupled, their noisy actuation will be very correlated.

Assuming that our simulator is deterministic implies that for each parameter vector there is a unique dataset. It is natural then to propose that if two datasets were generated with similar parameter vectors, then their proxies must be actuated with correlated noise. This all means that the noise correlation matrix 𝚺 (See the equation) depends on the parameter vectors. The closer the parameter vectors the higher the entry in the correlation matrix.

3.2 The gpemulator Object

At this point all the ingredients are in place. Our emulator is a GP associated with a system of linear models (independent proxies) which are actuated with correlated noise (the coupling).

The package includes a simple object called gpemulator, that simplifies the use of an emulator. Currently it does not reduce the work required for each of the items described in the previous section, but it helps putting them all together.

The interface of the constructor @gpemulator/gpemulator tries to follow the style of the GPML Toolbox. So instantiation of the object looks like

# Covariance function definition
covfunc  = {@covFunc_emu};

# Mean function definition
meanfunc = {@meanFun_emu};

emu = gpemulator ('covFunc', covfunc,'meanFunc', meanfunc);

The last line is the actual construction of the emulator and we need to pass a covariance function and a mean function. The previous lines define the mean and the covariance functions, which are just containers of function handles, that we will discuss in detail in the examples below.

To train (or condition) the emulator on some given datasets Y we call the method @gpemulator/designdata

emu = designdata (emu, Y);

After this the emulator is able to make predictions at inputs X using the method @gpemulator/emulate

y = emulate (emu, X);

3.3 Emulating the Van der Pol Oscillator Limit Cycle

The full code for this example can be found in the script file emulator_VdP.m inside the examples folder of the package.

We continue with the example from the previous sections (See The Anatomy of an Emulator). We assume that we have the design data and the parameter mapping. The script to emulate will then have three parts. In the first part we load the design data and define some helper functions. In the second part we use the design data to condition our emulator. In the third part we emulate at unseen parameter values and time stamps.

#### Design data ####
fname = "vanderpol.dat";
ddata = load (fname);

vdp = @(y,t, mu) [y(2); ...                   # Simulator is
           mu * (1 - y(1)^2) * y(2) - y(1)];   # Van der Pol oscillator

y0        = [0.5 0];
simulator = @(t,p) cell2mat(
                    arrayfun (
  @(r)lsode (@(x,t)vdp(x,t,r), ...
  lsode(@(x,t)vdp(x,t,r), y0, [0 ddata.T])(end,:), ... # Initial condition
                    t)(:,1), ...
                    p(:).', 'unif', 0));

proxy = @(x) x.ic(:,1).' .* cos (x.lp.' .* x.time) + ...
              (x.ic(:,2)./x.lp).' .* sin (x.lp.' .* x.time); # Proxy

#### Emulator conditioning ####
pkg load gpml
pkg load gpemulation

covfunc  = {@covFunc_emu, {@covMaterniso,1},{@covSEiso}};
meanfunc = {proxy};

mdp     = log (min (tril (abs (ddata.nlp - ddata.nlp.'),-1,'pack')));
hyp.S   = [mdp+log(2) 0];
hyp.S0  = [mdp+log(1) 0];

emu     = gpemulator ('covFunc', {covfunc{:},hyp}, ...
                      'meanFunc', meanfunc, ...
                      'regParam', 0);

#### Condition the emulator on design data
par = ddata.parmap (ddata.nlp);
x   = struct ('time', ddata.t, ...
              'nlp', ddata.nlp, ...
              'lp', par(:,1), ...
              'ic', par(:,2:3), ...
              'idx', 1);

emu = designdata (emu, x, ddata.y(:,:,1));

#### Emulation ####
nT  = 2 * length (ddata.t);
t   = linspace (ddata.t(1), ddata.t(end)+2, nT).';
N   = 5;
nls = [min(ddata.nlp) max(ddata.nlp)];
nlp = linspace (nls(1)-exp(mdp), nls(2)+exp(mdp), N).';
Y   = simulator (t, nlp);

par = ddata.parmap (nlp);
x   = struct ('time', t, ...
              'nlp', nlp, ...
              'lp', par(:,1), ...
              'ic', par(:,2:3), ...
              'idx', 1);
yp  = emulate (emu, x);

In the first part, we define the function simulator that takes an array of parameters for the Van der Pol oscillator and returns trajectories on the limit cycle for each parameter value. Also we use a function proxy that takes an array of parameters and evaluates the proxy. The proxy function expects and input structure x with the fields time, ic and lp. The field time contains the time stamps, ic the initial conditions and lp the three parameters that define the proxy.

In the second part we define the emulator using proxy as a mean function and the function covFunc_emu for the covariance. This function is described in detail below. At this point it is relevant to know that among its inputs there are two other covariance functions that it will use to couple the proxies. In this case we are using the Matern kernel for the parameter coupling and the squared exponential for the initial conditions, both functions are from the the GPML Toolbox. The parameters for these functions are given in the hyp structure, fields hyp.S and hyp.S0, respectively.

At the end of this second part we condition the emulator on the design data by calling the method

x   = struct ('time', t, ...
              'nlp', nlp, ...
              'lp', par(:,1), ...
              'ic', par(:,2:3), ...
              'idx', 1);
emu = designdata (emu, x, ddata.y(:,:,1));

You might have noticed that the input structure x has also the field idx. This field contains the component of our proxy that we will use to emulate the design data. In this case we use only the first component of our emulator for the first component of the design data.

Finally, in the third part of the script we emulate at unseen parameter values and time stamps. The field x.time covers the observed time span plus an extra time window so we can see how good our emulator extrapolates outside of the data region. It also has the double number of time stamps, so we are intrapolating. The field x.nlp contains the unseen parameters at which we are emulating. This field contains values between the design data values, so again we are intrapolating. It also contains values outside the seen region of parameter values, i.e. extrapolation.

The extrapolation at low parameter values presents no difficulty for the emulator because the oscillator behaves more linearly there (see the top panel of the figure below). However, when we extrapolate at higher parameter values, the results get worse and you can see this in the bottom panel of the figure.

Emulation results

These are results using n = 80 trajectories of the oscillator and nT = 40 time samples.

3.3.1 Construction of the Covariance function of the emulator

In the example shown above the code for the actual covariance of our emulator was not shown. Here we show and describe the details of the implementation.

We want to use the function covLTIo1diag which we described before (See Covariance Function of a LSDE). The inputs argument to this function are the timestamps and the diagonalized representation of the emulator system matrix A = VDV⁻¹, that is we need to pass V, V⁻¹ and the diagonal of D. Therefore we first write a function that builds these matrices from a given vector of parameters. As we said before, our mean function

y(t) = A cos ( w t ) + (B/w) sin ( w t )

gives us the paramters for the system

 = z 

 = -w² y

with y(0) = A and (0) = z(0) = B.

The columns of the matrix V are

Vₙᵀ = [ 1 Dₙₙ]

where the diagonal of D is

Dₙₙ = (-1)ⁿiw

with i the imaginary unit.

The construction of these matrices is implemented in the function param2matrix that you see in the following code.

function A = param2matrix (lp)
  # Egiven-values and -vector functions
  Df  =@(l) l.*[-J; J];                 # Eigenvalues
  Vf  =@(l) [ones(1,2); Df(l).'];       # Eigenvectors
  Vif =@(l) 0.5 * [ones(2,1) 1./Df(l)]; # Inverse of Eigenvectors

  # Eigenspace
  V    = arrayfun (Vf, lp, "unif", 0);
  V(1) = sparse (V{1});
  V    = blkdiag (V{:});
  # Inverse eigenspace
  Vi    = arrayfun (Vif, lp, "unif", 0);
  Vi(1) = sparse (Vi{1});
  Vi    = blkdiag (Vi{:});
  # Eigenvalues
  D = cell2mat (arrayfun (Df, lp, "unif", 0));

  A = {V,D,Vi};

function k = covFunc_emu (covS, covS0, h, x, y=[])

  n = length(x.nlp);
  X = {x.time, param2matrix(x.lp)};
  idxx = [];
  if ~isempty (x.idx)
    idxx = cell2mat (arrayfun ( ...
           @(i)(i-1)*2 + x.idx, 1:n, ...
           "unif", 0));

  if isempty (y) || strcmp (y, 'diag')

    S  = feval (covS{:}, h.S, x.nlp);
    S0 = feval (covS0{:}, h.S0, x.ic);
    S  = kron (S, ones(2));
    S0 = kron (S0, ones(2));
    k  = covLTIo1diag (X, y, S ,S0, idxx);

    ny = length(y.nlp);
    Y  = {y.time, param2matrix(y.lp)};
    idxy = [];
    if ~isempty (y.idx)
      idxy = cell2mat (arrayfun ( ...
            @(i)(i-1)*2 + y.idx, 1:ny, ...
            "unif", 0));

    S  = feval (covS{:}, h.S, x.nlp, y.nlp);
    S0 = feval (covS0{:}, h.S0, x.ic, y.ic);
    S  = kron (S, ones(2));
    S0 = kron (S0, ones(2));
    k  = covLTIo1diag (X, Y, S ,S0, idxx, idxy);


The second function implements the covariance function using covLTIo1diag.

Next: , Previous: , Up: Top   [Contents]

4 Advanced Emulation

Next: , Previous: , Up: Top   [Contents]

Appendix A Attribution-ShareAlike 4.0 International license

Attribution-ShareAlike 4.0 International


Creative Commons Corporation ("Creative Commons") is not a law firm and
does not provide legal services or legal advice. Distribution of
Creative Commons public licenses does not create a lawyer-client or
other relationship. Creative Commons makes its licenses and related
information available on an "as-is" basis. Creative Commons gives no
warranties regarding its licenses, any material licensed under their
terms and conditions, or any related information. Creative Commons
disclaims all liability for damages resulting from their use to the
fullest extent possible.

Using Creative Commons Public Licenses

Creative Commons public licenses provide a standard set of terms and
conditions that creators and other rights holders may use to share
original works of authorship and other material subject to copyright
and certain other rights specified in the public license below. The
following considerations are for informational purposes only, are not
exhaustive, and do not form part of our licenses.

     Considerations for licensors: Our public licenses are
     intended for use by those authorized to give the public
     permission to use material in ways otherwise restricted by
     copyright and certain other rights. Our licenses are
     irrevocable. Licensors should read and understand the terms
     and conditions of the license they choose before applying it.
     Licensors should also secure all rights necessary before
     applying our licenses so that the public can reuse the
     material as expected. Licensors should clearly mark any
     material not subject to the license. This includes other CC-
     licensed material, or material used under an exception or
     limitation to copyright. More considerations for licensors:

     Considerations for the public: By using one of our public
     licenses, a licensor grants the public permission to use the
     licensed material under specified terms and conditions. If
     the licensor's permission is not necessary for any reason--for
     example, because of any applicable exception or limitation to
     copyright--then that use is not regulated by the license. Our
     licenses grant only permissions under copyright and certain
     other rights that a licensor has authority to grant. Use of
     the licensed material may still be restricted for other
     reasons, including because others have copyright or other
     rights in the material. A licensor may make special requests,
     such as asking that all changes be marked or described.
     Although not required by our licenses, you are encouraged to
     respect those requests where reasonable. More_considerations
     for the public:


Creative Commons Attribution-ShareAlike 4.0 International Public

By exercising the Licensed Rights (defined below), You accept and agree
to be bound by the terms and conditions of this Creative Commons
Attribution-ShareAlike 4.0 International Public License ("Public
License"). To the extent this Public License may be interpreted as a
contract, You are granted the Licensed Rights in consideration of Your
acceptance of these terms and conditions, and the Licensor grants You
such rights in consideration of benefits the Licensor receives from
making the Licensed Material available under these terms and

Section 1 -- Definitions.

  a. Adapted Material means material subject to Copyright and Similar
     Rights that is derived from or based upon the Licensed Material
     and in which the Licensed Material is translated, altered,
     arranged, transformed, or otherwise modified in a manner requiring
     permission under the Copyright and Similar Rights held by the
     Licensor. For purposes of this Public License, where the Licensed
     Material is a musical work, performance, or sound recording,
     Adapted Material is always produced where the Licensed Material is
     synched in timed relation with a moving image.

  b. Adapter's License means the license You apply to Your Copyright
     and Similar Rights in Your contributions to Adapted Material in
     accordance with the terms and conditions of this Public License.

  c. BY-SA Compatible License means a license listed at, approved by Creative
     Commons as essentially the equivalent of this Public License.

  d. Copyright and Similar Rights means copyright and/or similar rights
     closely related to copyright including, without limitation,
     performance, broadcast, sound recording, and Sui Generis Database
     Rights, without regard to how the rights are labeled or
     categorized. For purposes of this Public License, the rights
     specified in Section 2(b)(1)-(2) are not Copyright and Similar

  e. Effective Technological Measures means those measures that, in the
     absence of proper authority, may not be circumvented under laws
     fulfilling obligations under Article 11 of the WIPO Copyright
     Treaty adopted on December 20, 1996, and/or similar international

  f. Exceptions and Limitations means fair use, fair dealing, and/or
     any other exception or limitation to Copyright and Similar Rights
     that applies to Your use of the Licensed Material.

  g. License Elements means the license attributes listed in the name
     of a Creative Commons Public License. The License Elements of this
     Public License are Attribution and ShareAlike.

  h. Licensed Material means the artistic or literary work, database,
     or other material to which the Licensor applied this Public

  i. Licensed Rights means the rights granted to You subject to the
     terms and conditions of this Public License, which are limited to
     all Copyright and Similar Rights that apply to Your use of the
     Licensed Material and that the Licensor has authority to license.

  j. Licensor means the individual(s) or entity(ies) granting rights
     under this Public License.

  k. Share means to provide material to the public by any means or
     process that requires permission under the Licensed Rights, such
     as reproduction, public display, public performance, distribution,
     dissemination, communication, or importation, and to make material
     available to the public including in ways that members of the
     public may access the material from a place and at a time
     individually chosen by them.

  l. Sui Generis Database Rights means rights other than copyright
     resulting from Directive 96/9/EC of the European Parliament and of
     the Council of 11 March 1996 on the legal protection of databases,
     as amended and/or succeeded, as well as other essentially
     equivalent rights anywhere in the world.

  m. You means the individual or entity exercising the Licensed Rights
     under this Public License. Your has a corresponding meaning.

Section 2 -- Scope.

  a. License grant.

       1. Subject to the terms and conditions of this Public License,
          the Licensor hereby grants You a worldwide, royalty-free,
          non-sublicensable, non-exclusive, irrevocable license to
          exercise the Licensed Rights in the Licensed Material to:

            a. reproduce and Share the Licensed Material, in whole or
               in part; and

            b. produce, reproduce, and Share Adapted Material.

       2. Exceptions and Limitations. For the avoidance of doubt, where
          Exceptions and Limitations apply to Your use, this Public
          License does not apply, and You do not need to comply with
          its terms and conditions.

       3. Term. The term of this Public License is specified in Section

       4. Media and formats; technical modifications allowed. The
          Licensor authorizes You to exercise the Licensed Rights in
          all media and formats whether now known or hereafter created,
          and to make technical modifications necessary to do so. The
          Licensor waives and/or agrees not to assert any right or
          authority to forbid You from making technical modifications
          necessary to exercise the Licensed Rights, including
          technical modifications necessary to circumvent Effective
          Technological Measures. For purposes of this Public License,
          simply making modifications authorized by this Section 2(a)
          (4) never produces Adapted Material.

       5. Downstream recipients.

            a. Offer from the Licensor -- Licensed Material. Every
               recipient of the Licensed Material automatically
               receives an offer from the Licensor to exercise the
               Licensed Rights under the terms and conditions of this
               Public License.

            b. Additional offer from the Licensor -- Adapted Material.
               Every recipient of Adapted Material from You
               automatically receives an offer from the Licensor to
               exercise the Licensed Rights in the Adapted Material
               under the conditions of the Adapter's License You apply.

            c. No downstream restrictions. You may not offer or impose
               any additional or different terms or conditions on, or
               apply any Effective Technological Measures to, the
               Licensed Material if doing so restricts exercise of the
               Licensed Rights by any recipient of the Licensed

       6. No endorsement. Nothing in this Public License constitutes or
          may be construed as permission to assert or imply that You
          are, or that Your use of the Licensed Material is, connected
          with, or sponsored, endorsed, or granted official status by,
          the Licensor or others designated to receive attribution as
          provided in Section 3(a)(1)(A)(i).

  b. Other rights.

       1. Moral rights, such as the right of integrity, are not
          licensed under this Public License, nor are publicity,
          privacy, and/or other similar personality rights; however, to
          the extent possible, the Licensor waives and/or agrees not to
          assert any such rights held by the Licensor to the limited
          extent necessary to allow You to exercise the Licensed
          Rights, but not otherwise.

       2. Patent and trademark rights are not licensed under this
          Public License.

       3. To the extent possible, the Licensor waives any right to
          collect royalties from You for the exercise of the Licensed
          Rights, whether directly or through a collecting society
          under any voluntary or waivable statutory or compulsory
          licensing scheme. In all other cases the Licensor expressly
          reserves any right to collect such royalties.

Section 3 -- License Conditions.

Your exercise of the Licensed Rights is expressly made subject to the
following conditions.

  a. Attribution.

       1. If You Share the Licensed Material (including in modified
          form), You must:

            a. retain the following if it is supplied by the Licensor
               with the Licensed Material:

                 i. identification of the creator(s) of the Licensed
                    Material and any others designated to receive
                    attribution, in any reasonable manner requested by
                    the Licensor (including by pseudonym if

                ii. a copyright notice;

               iii. a notice that refers to this Public License;

                iv. a notice that refers to the disclaimer of

                 v. a URI or hyperlink to the Licensed Material to the
                    extent reasonably practicable;

            b. indicate if You modified the Licensed Material and
               retain an indication of any previous modifications; and

            c. indicate the Licensed Material is licensed under this
               Public License, and include the text of, or the URI or
               hyperlink to, this Public License.

       2. You may satisfy the conditions in Section 3(a)(1) in any
          reasonable manner based on the medium, means, and context in
          which You Share the Licensed Material. For example, it may be
          reasonable to satisfy the conditions by providing a URI or
          hyperlink to a resource that includes the required

       3. If requested by the Licensor, You must remove any of the
          information required by Section 3(a)(1)(A) to the extent
          reasonably practicable.

  b. ShareAlike.

     In addition to the conditions in Section 3(a), if You Share
     Adapted Material You produce, the following conditions also apply.

       1. The Adapter's License You apply must be a Creative Commons
          license with the same License Elements, this version or
          later, or a BY-SA Compatible License.

       2. You must include the text of, or the URI or hyperlink to, the
          Adapter's License You apply. You may satisfy this condition
          in any reasonable manner based on the medium, means, and
          context in which You Share Adapted Material.

       3. You may not offer or impose any additional or different terms
          or conditions on, or apply any Effective Technological
          Measures to, Adapted Material that restrict exercise of the
          rights granted under the Adapter's License You apply.

Section 4 -- Sui Generis Database Rights.

Where the Licensed Rights include Sui Generis Database Rights that
apply to Your use of the Licensed Material:

  a. for the avoidance of doubt, Section 2(a)(1) grants You the right
     to extract, reuse, reproduce, and Share all or a substantial
     portion of the contents of the database;

  b. if You include all or a substantial portion of the database
     contents in a database in which You have Sui Generis Database
     Rights, then the database in which You have Sui Generis Database
     Rights (but not its individual contents) is Adapted Material,

     including for purposes of Section 3(b); and
  c. You must comply with the conditions in Section 3(a) if You Share
     all or a substantial portion of the contents of the database.

For the avoidance of doubt, this Section 4 supplements and does not
replace Your obligations under this Public License where the Licensed
Rights include other Copyright and Similar Rights.

Section 5 -- Disclaimer of Warranties and Limitation of Liability.



  c. The disclaimer of warranties and limitation of liability provided
     above shall be interpreted in a manner that, to the extent
     possible, most closely approximates an absolute disclaimer and
     waiver of all liability.

Section 6 -- Term and Termination.

  a. This Public License applies for the term of the Copyright and
     Similar Rights licensed here. However, if You fail to comply with
     this Public License, then Your rights under this Public License
     terminate automatically.

  b. Where Your right to use the Licensed Material has terminated under
     Section 6(a), it reinstates:

       1. automatically as of the date the violation is cured, provided
          it is cured within 30 days of Your discovery of the
          violation; or

       2. upon express reinstatement by the Licensor.

     For the avoidance of doubt, this Section 6(b) does not affect any
     right the Licensor may have to seek remedies for Your violations
     of this Public License.

  c. For the avoidance of doubt, the Licensor may also offer the
     Licensed Material under separate terms or conditions or stop
     distributing the Licensed Material at any time; however, doing so
     will not terminate this Public License.

  d. Sections 1, 5, 6, 7, and 8 survive termination of this Public

Section 7 -- Other Terms and Conditions.

  a. The Licensor shall not be bound by any additional or different
     terms or conditions communicated by You unless expressly agreed.

  b. Any arrangements, understandings, or agreements regarding the
     Licensed Material not stated herein are separate from and
     independent of the terms and conditions of this Public License.

Section 8 -- Interpretation.

  a. For the avoidance of doubt, this Public License does not, and
     shall not be interpreted to, reduce, limit, restrict, or impose
     conditions on any use of the Licensed Material that could lawfully
     be made without permission under this Public License.

  b. To the extent possible, if any provision of this Public License is
     deemed unenforceable, it shall be automatically reformed to the
     minimum extent necessary to make it enforceable. If the provision
     cannot be reformed, it shall be severed from this Public License
     without affecting the enforceability of the remaining terms and

  c. No term or condition of this Public License will be waived and no
     failure to comply consented to unless expressly agreed to by the

  d. Nothing in this Public License constitutes or may be interpreted
     as a limitation upon, or waiver of, any privileges and immunities
     that apply to the Licensor or You, including from the legal
     processes of any jurisdiction or authority.


Creative Commons is not a party to its public
licenses. Notwithstanding, Creative Commons may elect to apply one of
its public licenses to material it publishes and in those instances
will be considered the “Licensor.” The text of the Creative Commons
public licenses is dedicated to the public domain under the CC0 Public
Domain Dedication. Except for the limited purpose of indicating that
material is shared under a Creative Commons public license or as
otherwise permitted by the Creative Commons policies published at, Creative Commons does not authorize the
use of the trademark "Creative Commons" or any other trademark or logo
of Creative Commons without its prior written consent including,
without limitation, in connection with any unauthorized modifications
to any of its public licenses or any other arrangements,
understandings, or agreements concerning use of licensed material. For
the avoidance of doubt, this paragraph does not form part of the
public licenses.

Creative Commons may be contacted at

Previous: , Up: Top   [Contents]

Appendix B GNU General Public License

Version 3, 29 June 2007
Copyright © 2007 Free Software Foundation, Inc.

Everyone is permitted to copy and distribute verbatim copies of this
license document, but changing it is not allowed.


The GNU General Public License is a free, copyleft license for software and other kinds of works.

The licenses for most software and other practical works are designed to take away your freedom to share and change the works. By contrast, the GNU General Public License is intended to guarantee your freedom to share and change all versions of a program—to make sure it remains free software for all its users. We, the Free Software Foundation, use the GNU General Public License for most of our software; it applies also to any other work released this way by its authors. You can apply it to your programs, too.

When we speak of free software, we are referring to freedom, not price. Our General Public Licenses are designed to make sure that you have the freedom to distribute copies of free software (and charge for them if you wish), that you receive source code or can get it if you want it, that you can change the software or use pieces of it in new free programs, and that you know you can do these things.

To protect your rights, we need to prevent others from denying you these rights or asking you to surrender the rights. Therefore, you have certain responsibilities if you distribute copies of the software, or if you modify it: responsibilities to respect the freedom of others.

For example, if you distribute copies of such a program, whether gratis or for a fee, you must pass on to the recipients the same freedoms that you received. You must make sure that they, too, receive or can get the source code. And you must show them these terms so they know their rights.

Developers that use the GNU GPL protect your rights with two steps: (1) assert copyright on the software, and (2) offer you this License giving you legal permission to copy, distribute and/or modify it.

For the developers’ and authors’ protection, the GPL clearly explains that there is no warranty for this free software. For both users’ and authors’ sake, the GPL requires that modified versions be marked as changed, so that their problems will not be attributed erroneously to authors of previous versions.

Some devices are designed to deny users access to install or run modified versions of the software inside them, although the manufacturer can do so. This is fundamentally incompatible with the aim of protecting users’ freedom to change the software. The systematic pattern of such abuse occurs in the area of products for individuals to use, which is precisely where it is most unacceptable. Therefore, we have designed this version of the GPL to prohibit the practice for those products. If such problems arise substantially in other domains, we stand ready to extend this provision to those domains in future versions of the GPL, as needed to protect the freedom of users.

Finally, every program is threatened constantly by software patents. States should not allow patents to restrict development and use of software on general-purpose computers, but in those that do, we wish to avoid the special danger that patents applied to a free program could make it effectively proprietary. To prevent this, the GPL assures that patents cannot be used to render the program non-free.

The precise terms and conditions for copying, distribution and modification follow.


  1. Definitions.

    “This License” refers to version 3 of the GNU General Public License.

    “Copyright” also means copyright-like laws that apply to other kinds of works, such as semiconductor masks.

    “The Program” refers to any copyrightable work licensed under this License. Each licensee is addressed as “you”. “Licensees” and “recipients” may be individuals or organizations.

    To “modify” a work means to copy from or adapt all or part of the work in a fashion requiring copyright permission, other than the making of an exact copy. The resulting work is called a “modified version” of the earlier work or a work “based on” the earlier work.

    A “covered work” means either the unmodified Program or a work based on the Program.

    To “propagate” a work means to do anything with it that, without permission, would make you directly or secondarily liable for infringement under applicable copyright law, except executing it on a computer or modifying a private copy. Propagation includes copying, distribution (with or without modification), making available to the public, and in some countries other activities as well.

    To “convey” a work means any kind of propagation that enables other parties to make or receive copies. Mere interaction with a user through a computer network, with no transfer of a copy, is not conveying.

    An interactive user interface displays “Appropriate Legal Notices” to the extent that it includes a convenient and prominently visible feature that (1) displays an appropriate copyright notice, and (2) tells the user that there is no warranty for the work (except to the extent that warranties are provided), that licensees may convey the work under this License, and how to view a copy of this License. If the interface presents a list of user commands or options, such as a menu, a prominent item in the list meets this criterion.

  2. Source Code.

    The “source code” for a work means the preferred form of the work for making modifications to it. “Object code” means any non-source form of a work.

    A “Standard Interface” means an interface that either is an official standard defined by a recognized standards body, or, in the case of interfaces specified for a particular programming language, one that is widely used among developers working in that language.

    The “System Libraries” of an executable work include anything, other than the work as a whole, that (a) is included in the normal form of packaging a Major Component, but which is not part of that Major Component, and (b) serves only to enable use of the work with that Major Component, or to implement a Standard Interface for which an implementation is available to the public in source code form. A “Major Component”, in this context, means a major essential component (kernel, window system, and so on) of the specific operating system (if any) on which the executable work runs, or a compiler used to produce the work, or an object code interpreter used to run it.

    The “Corresponding Source” for a work in object code form means all the source code needed to generate, install, and (for an executable work) run the object code and to modify the work, including scripts to control those activities. However, it does not include the work’s System Libraries, or general-purpose tools or generally available free programs which are used unmodified in performing those activities but which are not part of the work. For example, Corresponding Source includes interface definition files associated with source files for the work, and the source code for shared libraries and dynamically linked subprograms that the work is specifically designed to require, such as by intimate data communication or control flow between those subprograms and other parts of the work.

    The Corresponding Source need not include anything that users can regenerate automatically from other parts of the Corresponding Source.

    The Corresponding Source for a work in source code form is that same work.

  3. Basic Permissions.

    All rights granted under this License are granted for the term of copyright on the Program, and are irrevocable provided the stated conditions are met. This License explicitly affirms your unlimited permission to run the unmodified Program. The output from running a covered work is covered by this License only if the output, given its content, constitutes a covered work. This License acknowledges your rights of fair use or other equivalent, as provided by copyright law.

    You may make, run and propagate covered works that you do not convey, without conditions so long as your license otherwise remains in force. You may convey covered works to others for the sole purpose of having them make modifications exclusively for you, or provide you with facilities for running those works, provided that you comply with the terms of this License in conveying all material for which you do not control copyright. Those thus making or running the covered works for you must do so exclusively on your behalf, under your direction and control, on terms that prohibit them from making any copies of your copyrighted material outside their relationship with you.

    Conveying under any other circumstances is permitted solely under the conditions stated below. Sublicensing is not allowed; section 10 makes it unnecessary.

  4. Protecting Users’ Legal Rights From Anti-Circumvention Law.

    No covered work shall be deemed part of an effective technological measure under any applicable law fulfilling obligations under article 11 of the WIPO copyright treaty adopted on 20 December 1996, or similar laws prohibiting or restricting circumvention of such measures.

    When you convey a covered work, you waive any legal power to forbid circumvention of technological measures to the extent such circumvention is effected by exercising rights under this License with respect to the covered work, and you disclaim any intention to limit operation or modification of the work as a means of enforcing, against the work’s users, your or third parties’ legal rights to forbid circumvention of technological measures.

  5. Conveying Verbatim Copies.

    You may convey verbatim copies of the Program’s source code as you receive it, in any medium, provided that you conspicuously and appropriately publish on each copy an appropriate copyright notice; keep intact all notices stating that this License and any non-permissive terms added in accord with section 7 apply to the code; keep intact all notices of the absence of any warranty; and give all recipients a copy of this License along with the Program.

    You may charge any price or no price for each copy that you convey, and you may offer support or warranty protection for a fee.

  6. Conveying Modified Source Versions.

    You may convey a work based on the Program, or the modifications to produce it from the Program, in the form of source code under the terms of section 4, provided that you also meet all of these conditions:

    1. The work must carry prominent notices stating that you modified it, and giving a relevant date.
    2. The work must carry prominent notices stating that it is released under this License and any conditions added under section 7. This requirement modifies the requirement in section 4 to “keep intact all notices”.
    3. You must license the entire work, as a whole, under this License to anyone who comes into possession of a copy. This License will therefore apply, along with any applicable section 7 additional terms, to the whole of the work, and all its parts, regardless of how they are packaged. This License gives no permission to license the work in any other way, but it does not invalidate such permission if you have separately received it.
    4. If the work has interactive user interfaces, each must display Appropriate Legal Notices; however, if the Program has interactive interfaces that do not display Appropriate Legal Notices, your work need not make them do so.

    A compilation of a covered work with other separate and independent works, which are not by their nature extensions of the covered work, and which are not combined with it such as to form a larger program, in or on a volume of a storage or distribution medium, is called an “aggregate” if the compilation and its resulting copyright are not used to limit the access or legal rights of the compilation’s users beyond what the individual works permit. Inclusion of a covered work in an aggregate does not cause this License to apply to the other parts of the aggregate.

  7. Conveying Non-Source Forms.

    You may convey a covered work in object code form under the terms of sections 4 and 5, provided that you also convey the machine-readable Corresponding Source under the terms of this License, in one of these ways:

    1. Convey the object code in, or embodied in, a physical product (including a physical distribution medium), accompanied by the Corresponding Source fixed on a durable physical medium customarily used for software interchange.
    2. Convey the object code in, or embodied in, a physical product (including a physical distribution medium), accompanied by a written offer, valid for at least three years and valid for as long as you offer spare parts or customer support for that product model, to give anyone who possesses the object code either (1) a copy of the Corresponding Source for all the software in the product that is covered by this License, on a durable physical medium customarily used for software interchange, for a price no more than your reasonable cost of physically performing this conveying of source, or (2) access to copy the Corresponding Source from a network server at no charge.
    3. Convey individual copies of the object code with a copy of the written offer to provide the Corresponding Source. This alternative is allowed only occasionally and noncommercially, and only if you received the object code with such an offer, in accord with subsection 6b.
    4. Convey the object code by offering access from a designated place (gratis or for a charge), and offer equivalent access to the Corresponding Source in the same way through the same place at no further charge. You need not require recipients to copy the Corresponding Source along with the object code. If the place to copy the object code is a network server, the Corresponding Source may be on a different server (operated by you or a third party) that supports equivalent copying facilities, provided you maintain clear directions next to the object code saying where to find the Corresponding Source. Regardless of what server hosts the Corresponding Source, you remain obligated to ensure that it is available for as long as needed to satisfy these requirements.
    5. Convey the object code using peer-to-peer transmission, provided you inform other peers where the object code and Corresponding Source of the work are being offered to the general public at no charge under subsection 6d.

    A separable portion of the object code, whose source code is excluded from the Corresponding Source as a System Library, need not be included in conveying the object code work.

    A “User Product” is either (1) a “consumer product”, which means any tangible personal property which is normally used for personal, family, or household purposes, or (2) anything designed or sold for incorporation into a dwelling. In determining whether a product is a consumer product, doubtful cases shall be resolved in favor of coverage. For a particular product received by a particular user, “normally used” refers to a typical or common use of that class of product, regardless of the status of the particular user or of the way in which the particular user actually uses, or expects or is expected to use, the product. A product is a consumer product regardless of whether the product has substantial commercial, industrial or non-consumer uses, unless such uses represent the only significant mode of use of the product.

    “Installation Information” for a User Product means any methods, procedures, authorization keys, or other information required to install and execute modified versions of a covered work in that User Product from a modified version of its Corresponding Source. The information must suffice to ensure that the continued functioning of the modified object code is in no case prevented or interfered with solely because modification has been made.

    If you convey an object code work under this section in, or with, or specifically for use in, a User Product, and the conveying occurs as part of a transaction in which the right of possession and use of the User Product is transferred to the recipient in perpetuity or for a fixed term (regardless of how the transaction is characterized), the Corresponding Source conveyed under this section must be accompanied by the Installation Information. But this requirement does not apply if neither you nor any third party retains the ability to install modified object code on the User Product (for example, the work has been installed in ROM).

    The requirement to provide Installation Information does not include a requirement to continue to provide support service, warranty, or updates for a work that has been modified or installed by the recipient, or for the User Product in which it has been modified or installed. Access to a network may be denied when the modification itself materially and adversely affects the operation of the network or violates the rules and protocols for communication across the network.

    Corresponding Source conveyed, and Installation Information provided, in accord with this section must be in a format that is publicly documented (and with an implementation available to the public in source code form), and must require no special password or key for unpacking, reading or copying.

  8. Additional Terms.

    “Additional permissions” are terms that supplement the terms of this License by making exceptions from one or more of its conditions. Additional permissions that are applicable to the entire Program shall be treated as though they were included in this License, to the extent that they are valid under applicable law. If additional permissions apply only to part of the Program, that part may be used separately under those permissions, but the entire Program remains governed by this License without regard to the additional permissions.

    When you convey a copy of a covered work, you may at your option remove any additional permissions from that copy, or from any part of it. (Additional permissions may be written to require their own removal in certain cases when you modify the work.) You may place additional permissions on material, added by you to a covered work, for which you have or can give appropriate copyright permission.

    Notwithstanding any other provision of this License, for material you add to a covered work, you may (if authorized by the copyright holders of that material) supplement the terms of this License with terms:

    1. Disclaiming warranty or limiting liability differently from the terms of sections 15 and 16 of this License; or
    2. Requiring preservation of specified reasonable legal notices or author attributions in that material or in the Appropriate Legal Notices displayed by works containing it; or
    3. Prohibiting misrepresentation of the origin of that material, or requiring that modified versions of such material be marked in reasonable ways as different from the original version; or
    4. Limiting the use for publicity purposes of names of licensors or authors of the material; or
    5. Declining to grant rights under trademark law for use of some trade names, trademarks, or service marks; or
    6. Requiring indemnification of licensors and authors of that material by anyone who conveys the material (or modified versions of it) with contractual assumptions of liability to the recipient, for any liability that these contractual assumptions directly impose on those licensors and authors.

    All other non-permissive additional terms are considered “further restrictions” within the meaning of section 10. If the Program as you received it, or any part of it, contains a notice stating that it is governed by this License along with a term that is a further restriction, you may remove that term. If a license document contains a further restriction but permits relicensing or conveying under this License, you may add to a covered work material governed by the terms of that license document, provided that the further restriction does not survive such relicensing or conveying.

    If you add terms to a covered work in accord with this section, you must place, in the relevant source files, a statement of the additional terms that apply to those files, or a notice indicating where to find the applicable terms.

    Additional terms, permissive or non-permissive, may be stated in the form of a separately written license, or stated as exceptions; the above requirements apply either way.

  9. Termination.

    You may not propagate or modify a covered work except as expressly provided under this License. Any attempt otherwise to propagate or modify it is void, and will automatically terminate your rights under this License (including any patent licenses granted under the third paragraph of section 11).

    However, if you cease all violation of this License, then your license from a particular copyright holder is reinstated (a) provisionally, unless and until the copyright holder explicitly and finally terminates your license, and (b) permanently, if the copyright holder fails to notify you of the violation by some reasonable means prior to 60 days after the cessation.

    Moreover, your license from a particular copyright holder is reinstated permanently if the copyright holder notifies you of the violation by some reasonable means, this is the first time you have received notice of violation of this License (for any work) from that copyright holder, and you cure the violation prior to 30 days after your receipt of the notice.

    Termination of your rights under this section does not terminate the licenses of parties who have received copies or rights from you under this License. If your rights have been terminated and not permanently reinstated, you do not qualify to receive new licenses for the same material under section 10.

  10. Acceptance Not Required for Having Copies.

    You are not required to accept this License in order to receive or run a copy of the Program. Ancillary propagation of a covered work occurring solely as a consequence of using peer-to-peer transmission to receive a copy likewise does not require acceptance. However, nothing other than this License grants you permission to propagate or modify any covered work. These actions infringe copyright if you do not accept this License. Therefore, by modifying or propagating a covered work, you indicate your acceptance of this License to do so.

  11. Automatic Licensing of Downstream Recipients.

    Each time you convey a covered work, the recipient automatically receives a license from the original licensors, to run, modify and propagate that work, subject to this License. You are not responsible for enforcing compliance by third parties with this License.

    An “entity transaction” is a transaction transferring control of an organization, or substantially all assets of one, or subdividing an organization, or merging organizations. If propagation of a covered work results from an entity transaction, each party to that transaction who receives a copy of the work also receives whatever licenses to the work the party’s predecessor in interest had or could give under the previous paragraph, plus a right to possession of the Corresponding Source of the work from the predecessor in interest, if the predecessor has it or can get it with reasonable efforts.

    You may not impose any further restrictions on the exercise of the rights granted or affirmed under this License. For example, you may not impose a license fee, royalty, or other charge for exercise of rights granted under this License, and you may not initiate litigation (including a cross-claim or counterclaim in a lawsuit) alleging that any patent claim is infringed by making, using, selling, offering for sale, or importing the Program or any portion of it.

  12. Patents.

    A “contributor” is a copyright holder who authorizes use under this License of the Program or a work on which the Program is based. The work thus licensed is called the contributor’s “contributor version”.

    A contributor’s “essential patent claims” are all patent claims owned or controlled by the contributor, whether already acquired or hereafter acquired, that would be infringed by some manner, permitted by this License, of making, using, or selling its contributor version, but do not include claims that would be infringed only as a consequence of further modification of the contributor version. For purposes of this definition, “control” includes the right to grant patent sublicenses in a manner consistent with the requirements of this License.

    Each contributor grants you a non-exclusive, worldwide, royalty-free patent license under the contributor’s essential patent claims, to make, use, sell, offer for sale, import and otherwise run, modify and propagate the contents of its contributor version.

    In the following three paragraphs, a “patent license” is any express agreement or commitment, however denominated, not to enforce a patent (such as an express permission to practice a patent or covenant not to sue for patent infringement). To “grant” such a patent license to a party means to make such an agreement or commitment not to enforce a patent against the party.

    If you convey a covered work, knowingly relying on a patent license, and the Corresponding Source of the work is not available for anyone to copy, free of charge and under the terms of this License, through a publicly available network server or other readily accessible means, then you must either (1) cause the Corresponding Source to be so available, or (2) arrange to deprive yourself of the benefit of the patent license for this particular work, or (3) arrange, in a manner consistent with the requirements of this License, to extend the patent license to downstream recipients. “Knowingly relying” means you have actual knowledge that, but for the patent license, your conveying the covered work in a country, or your recipient’s use of the covered work in a country, would infringe one or more identifiable patents in that country that you have reason to believe are valid.

    If, pursuant to or in connection with a single transaction or arrangement, you convey, or propagate by procuring conveyance of, a covered work, and grant a patent license to some of the parties receiving the covered work authorizing them to use, propagate, modify or convey a specific copy of the covered work, then the patent license you grant is automatically extended to all recipients of the covered work and works based on it.

    A patent license is “discriminatory” if it does not include within the scope of its coverage, prohibits the exercise of, or is conditioned on the non-exercise of one or more of the rights that are specifically granted under this License. You may not convey a covered work if you are a party to an arrangement with a third party that is in the business of distributing software, under which you make payment to the third party based on the extent of your activity of conveying the work, and under which the third party grants, to any of the parties who would receive the covered work from you, a discriminatory patent license (a) in connection with copies of the covered work conveyed by you (or copies made from those copies), or (b) primarily for and in connection with specific products or compilations that contain the covered work, unless you entered into that arrangement, or that patent license was granted, prior to 28 March 2007.

    Nothing in this License shall be construed as excluding or limiting any implied license or other defenses to infringement that may otherwise be available to you under applicable patent law.

  13. No Surrender of Others’ Freedom.

    If conditions are imposed on you (whether by court order, agreement or otherwise) that contradict the conditions of this License, they do not excuse you from the conditions of this License. If you cannot convey a covered work so as to satisfy simultaneously your obligations under this License and any other pertinent obligations, then as a consequence you may not convey it at all. For example, if you agree to terms that obligate you to collect a royalty for further conveying from those to whom you convey the Program, the only way you could satisfy both those terms and this License would be to refrain entirely from conveying the Program.

  14. Use with the GNU Affero General Public License.

    Notwithstanding any other provision of this License, you have permission to link or combine any covered work with a work licensed under version 3 of the GNU Affero General Public License into a single combined work, and to convey the resulting work. The terms of this License will continue to apply to the part which is the covered work, but the special requirements of the GNU Affero General Public License, section 13, concerning interaction through a network will apply to the combination as such.

  15. Revised Versions of this License.

    The Free Software Foundation may publish revised and/or new versions of the GNU General Public License from time to time. Such new versions will be similar in spirit to the present version, but may differ in detail to address new problems or concerns.

    Each version is given a distinguishing version number. If the Program specifies that a certain numbered version of the GNU General Public License “or any later version” applies to it, you have the option of following the terms and conditions either of that numbered version or of any later version published by the Free Software Foundation. If the Program does not specify a version number of the GNU General Public License, you may choose any version ever published by the Free Software Foundation.

    If the Program specifies that a proxy can decide which future versions of the GNU General Public License can be used, that proxy’s public statement of acceptance of a version permanently authorizes you to choose that version for the Program.

    Later license versions may give you additional or different permissions. However, no additional obligations are imposed on any author or copyright holder as a result of your choosing to follow a later version.

  16. Disclaimer of Warranty.


  17. Limitation of Liability.


  18. Interpretation of Sections 15 and 16.

    If the disclaimer of warranty and limitation of liability provided above cannot be given local legal effect according to their terms, reviewing courts shall apply local law that most closely approximates an absolute waiver of all civil liability in connection with the Program, unless a warranty or assumption of liability accompanies a copy of the Program in return for a fee.


How to Apply These Terms to Your New Programs

If you develop a new program, and you want it to be of the greatest possible use to the public, the best way to achieve this is to make it free software which everyone can redistribute and change under these terms.

To do so, attach the following notices to the program. It is safest to attach them to the start of each source file to most effectively state the exclusion of warranty; and each file should have at least the “copyright” line and a pointer to where the full notice is found.

one line to give the program's name and a brief idea of what it does.
Copyright (C) year name of author

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see

Also add information on how to contact you by electronic and paper mail.

If the program does terminal interaction, make it output a short notice like this when it starts in an interactive mode:

program Copyright (C) year name of author
This program comes with ABSOLUTELY NO WARRANTY; for details type ‘show w’.
This is free software, and you are welcome to redistribute it
under certain conditions; type ‘show c’ for details.

The hypothetical commands ‘show w’ and ‘show c’ should show the appropriate parts of the General Public License. Of course, your program’s commands might be different; for a GUI interface, you would use an “about box”.

You should also get your employer (if you work as a programmer) or school, if any, to sign a “copyright disclaimer” for the program, if necessary. For more information on this, and how to apply and follow the GNU GPL, see

The GNU General Public License does not permit incorporating your program into proprietary programs. If your program is a subroutine library, you may consider it more useful to permit linking proprietary applications with the library. If this is what you want to do, use the GNU Lesser General Public License instead of this License. But first, please read



Note that the noise terms are uncorrelated in time, i.e. cov[ξ(t), ξ(r)] = 𝛴𝛿(t-r).