%% Multi-Output Gaussian Stochastic Processes, part 1: Representations, realizations, principal components (Karhunen-Loeve eigenfunctions)
% © 2024, Antonio Sala, Universitat Politècnica de València, Spain. All rights reserved.
% Objectives: illustrate the concept and representation of a 2-output GP, via an example.
% Presentations in video:
% http://personales.upv.es/asala/YT/V/mimoGP1EN.html
% http://personales.upv.es/asala/YT/V/mimogp2EN.html
% http://personales.upv.es/asala/YT/V/mimogp3pcaEN.html
%
%
% Table of Contents
% Basic idea, motivation
% Representations of multi-output GP
%     As an (n+1)-dimensional input , single output GP
%     As an "independent component/factor" representation
% Samples, realizations (example)
%     Representation 1: full covariance 2x2 matrix kernel
%     Representation 2: from component realizations
% Karhunen-Loeve eigenfunctions (principal component analysis if finite number of samples)
%
%% Basic idea, motivation
% In some applications we may have two "functions" of a given argument real x, or even multidimensional input, and functions are f_1(x) and f_2(x) which somehow are "correlated".
%
% Examples:
% \bullet1\bullet Weather sensor network (input would be latitude-longitude)...
% Temperature and Humidity at "same" points. Isotopic.
% Heterotopic: Temperature at some points, humidity at other different points.
%
% \bullet2\bullet Other interesting cases in engineering: multioutput " f_1(t)=position, f_2(t)=velocity", in a 2nd order mechanical system subject to white-noise acceleration, say.
%
%
%
%
% If we consider a statistical GP model of that setup, we may wish to generate realizations, or maybe we have samples f_1(x_k), f_2(x_k) and we wish to "interpolate" them; it will be better exploiting the "correlation" between f_1 and f_2 to improve our estimations. For example, we will generalize the idea of estimating speed from position measurements, or something like that.
%
%
%% Representations of multi-output GP
% Multi-output (say, 2 outputs) may be consider either a function f:\mathbb R\mapsto \mathbb R^2 or a function f:\mathbb R\times \{1,2\}\mapsto \mathbb R.
% We may have "independent" GPs k_1(x,x'), k_2(x,x'), or "correlated" ones k(x,x') being a 2x2 matrix... independent case: diagonal matrix, we would just think in f_1:\mathbb R\mapsto \mathbb R, f_2:\mathbb R\mapsto \mathbb R, separately.
% *Note: The ODE or PDE/based representation (stochastic differential equations SDE, signal processing, filtering) will NOT be addressed here: this PDF is just introductory/motivational, and the topic of SDE  has room for a whole course in itself.
%
%
% Let us detail, assuming zero mean everywhere to avoid notational clutter and to keep things simple.
%% As an (n+1)-dimensional input , single output GP
% If we consider f:\mathbb R\times \{1,2\}\mapsto \mathbb R, then we need a covariance kernel k:\{R\times \{1,2\}\}^2\mapsto \mathbb R, expressed as, say, k((x,i),(x',i')) or, well, just k(x,i,x',i').
% Once we think in the "(n+1)-dimensional input, single-output" GP, there is NOTHING NEW UNDER THE HOOD: everything from single-ouput GPs applies to multi-ouput GPs (prediction, principal components, etc.). So theory ends here... The rest is "plotting commands for better intuition", but nothing "fundamental" from a conceptual point of view.
% *Note: when one output is related to the other one via an ODE (position--speed, filtering, etc.) there IS more theory under the hood. These issues are NOT in the scope of the present material.
%
%
%
% *Two ways of "stacking" when building covariance matrices of a finite set of data: first "i", then "x", or first "x" then "i"... f:\mathbb R\times \{1,2\}\mapsto \mathbb R or f: \{1,2\}\times\mathbb R\mapsto \mathbb R. It's the same, anyway, but I just mention it to avoid confusion when interpreting the Matlab code below.
% How to obtain these covariance matrices? We cannot "invent" them, as they must be positive-definite for any finite set of data... well, we may consider (x,i) as an "extended" 1+dimensional vector, and use any usual kernel, say
% k(x,i,x',j)=Me^{-\alpha\cdot (x-x')^2-\beta\cdot (i-j)^2. It would be OK, but there is another option, see below...
%
%% As an "independent component/factor" representation
% We consider a "component" being sampled from "independent factor GPs u_i(x)", each one with a different kernel, if so wished. Later on, obtain f_1(x) and f_2(x) via a linear transformation  (a matrix C) of those independent factors u_i(x).
% Thus, we will assume, for instance: f_{2\times 1}(x)=C_{2\times 3}\cdot u_{3\times 1}(x), being u_i(x) realizations of three independent GPs.

U=cell(3,1);
sgu1=3;sgu2=0.75; sgu3=.3; %stationary variances
dd1=2;  dd2=0.65; dd3=0.2; %distance parameters
U{1}=@(x1,x2) sgu1^2*exp(-norm(x2-x1)^2/dd1^2);%kernel for component 1
U{2}=@(x1,x2) sgu2^2*exp(-norm(x2-x1)^2/dd2^2); %kernel 2
U{3}=@(x1,x2) sgu3^2*exp(-norm(x2-x1)^2/dd3^2);%kernel 3
C=[1 1 0;1 0 1]

% The "factors" may be considered "latent" variables, which may be worth estimating if so wished (see below), if there is any "physical" interpretation of interest for them.
% Note that we are NOT discussing "how" or "why" we arrived to this statistical GP model... it comes from whoever studied the data-based or theoretical modelling phase.
% Covariance matrix computation: So, f_{2\times 1}(x)=C_{2\times 3}\cdot u_{3\times 1}(x) would entail that cov(f(x),f(x'))=E[f(x)\cdot f(x')^T]=C\cdot E[u(x)\cdot u(x')^T]\cdot C^T= C\cdot K(x,x')\cdot C^T
% being K(x,x') a diagonal matrix with the covariance kernels of each of the "independent" components".

ComponentCovariance=@(x1,x2) diag([U{1}(x1,x2) U{2}(x1,x2) U{3}(x1,x2)]);
FkernelMatrix=@(x1,x2) C*ComponentCovariance(x1,x2)*C'; %MATRIX output [UNUSED]

% It can be easily seen, developing the matrix product expressions, that cov(f_i(x),f_j(x'))=E[f(x)\cdot f(x')^T]=C_{i,:}\cdot E[u(x)\cdot u(x')^T]\cdot C^T= C\cdot K(x,x')\cdot (C_{j,:})^T
% Let's program that as a function for later use:

Fkernel=@(x1,i1,x2,i2) C(i1,:)*ComponentCovariance(x1,x2)*C(i2,:)'; %Scalar output of Fkernel

%
%% Samples, realizations (example)

%Xtest=(-5:0.05:5)';
Xtest=(-4:0.05:4)';
Ntest=length(Xtest)

%% Representation 1: full covariance 2x2 matrix kernel
% We will stack input samples as "x index" in rows [Xtest \ 1; Xtest\ 2; \cdots ]:

XtestTwice=[Xtest; Xtest];
XtestLabels=kron([1;2],ones(size(Xtest)));

% We compute the full covariance matrix:

Ktest=K(XtestTwice,XtestLabels,XtestTwice,XtestLabels, Fkernel);
size(Ktest)
Samples=mvnrnd(zeros(length(XtestTwice),1),Ktest);
plot(Xtest,[Samples(1:Ntest); Samples((Ntest+1):end)],LineWidth=2.5), grid on
legend("f1","f2",Location="best"), title("Sample realizations")

%% Representation 2: from component realizations
% Of course, we might have preferred obtaining three separate SISO realizations for each of the u1(x) ,u2(x), u3(x) components and later on multiply by C. Results would be equivalent, clearly. Let's do it:

PriorVarU1=K_siso(Xtest,Xtest,U{1});
PriorVarU2=K_siso(Xtest,Xtest,U{2});
PriorVarU3=K_siso(Xtest,Xtest,U{3});
U1sample=mvnrnd(zeros(length(Xtest),1),PriorVarU1); size(U1sample)
U2sample=mvnrnd(zeros(length(Xtest),1),PriorVarU2); size(U2sample)
U3sample=mvnrnd(zeros(length(Xtest),1),PriorVarU3); size(U3sample)
plot(Xtest,[U1sample;U2sample;U3sample],LineWidth=2), grid on
legend("u1(x)","u2(x)","u3(x)",Location="best"), title("Latent Components U")
Fsample=C*[U1sample;U2sample;U3sample];
size(Fsample)
plot(Xtest,Fsample,LineWidth=2.5), grid on, title("Resulting outputs F=CU")
hold on
plot(Xtest,U1sample,"-.",LineWidth=1.5)
hold off, legend("f1","f2","component u1, common",Location="best")

%
%% Karhunen-Loeve eigenfunctions (principal component analysis if finite number of samples)
% Let's analyze eigenvectors and eigenvalues of the covariance matrix, taking into account the 2-output structure:

[V,D]=eig(Ktest);
PrincStd=sqrt(diag(D)+3e-14);
[PrincStdSorted,idx]=sort(PrincStd,1,"descend");
bar(PrincStdSorted), grid on , xline(9)
PrincCom=V(:,idx)*diag(PrincStdSorted); %eigenvectors, reordered
for CompNumber=[1 58]
    plot(Xtest,PrincCom(1:Ntest,CompNumber),'b' )
    hold on
    plot(Xtest,PrincCom((Ntest+1):end,CompNumber),'r')
end
hold off, grid on, 
ylim([-2.3 2.6]), yline(0), legend("f1","f2")
title("KL Principal Components")

%
%

function M=K(x1,i1,x2,i2,KernelFunction)
%Kernelfunction returns a SCALAR
%data samples in x1 and x2 are ROWS
N1=size(x1,1);
N2=size(x2,1);
M=zeros(N1,N2);
for i=1:N1
    for j=1:N2
        M(i,j)=KernelFunction(x1(i,:), i1(i), x2(j,:), i2(j));
    end
end
end

function M=K_siso(x1,x2,KernelFunction)
%Kernelfunction returns a SCALAR
%data samples in x1 and x2 are ROWS
N1=size(x1,1);
N2=size(x2,1);
M=zeros(N1,N2);
for i=1:N1
    for j=1:N2
        M(i,j)=KernelFunction(x1(i,:), x2(j,:));
    end
end
end
