function MLmethods_multimodal_exp2
% Test Gaussian process regression and ridge regression on combined sEMG 
% andTMG data of experiment 2. Six movements are performed for 3 
% repetetions by 9 non-amputated participants. For GPR, RBF kernels with 
% combinations of KL divergence, logSPD and Euclidean distance (TMG) and
% SPD, logSPD and Euclidean distance (sEMG) are tested.
% Training data are taken as zero and full activations of all movements,
% while testing data are taken as transition states corresponding to the
% whole range of intermediate activation.
% Results are registered in results/MULTIMODALresults%d.txt (RMSE values) 
% and predictions are registered in MULTIMODALpredictions%d.mat wit d=1:9.
%
%
%
% Writing code takes time. Polishing it and making it available to others takes longer! 
% If some parts of the code were useful for your research of for a better understanding 
% of the algorithms, please cite the related publications.
%
% @article{Jaquier17MDPI,
%   author="Jaquier, N. and Connan, M. and Castellini, C. and Calinon, S.",
%   title="Combining electro- and tactile myography to improve hand and 
%		wrist activity detection in prostheses",
%   year="2017",
%	journal = "Technologies",
%	volume = "",
%	number = "",
%	article number = "",
% }
% 
% Copyright (c) 2017 Idiap Research Institute, http://idiap.ch/
% Written by Noémie Jaquier
% 
% This code is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 3 as
% published by the Free Software Foundation.
% 
% This code is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this code. If not, see <http://www.gnu.org/licenses/>.

for subjectId = 1:9
	%% Parameters
	
	filestr = sprintf('results/MULTIMODALresults%d.txt',subjectId);
	predstr = sprintf('results/MULTIMODALpredictions%d.mat',subjectId);
	
	nbIterCV = 10;
	testProp = 0.6;

	%% Load and preprocess training data
	fprintf('Loading and preprocessing data %d...',subjectId);
	load(sprintf('data/exp2/data_S%d.mat',subjectId));
	
	% Covariance features for EMG data
	% Sliding window for spatial covariance computation
	timeWindow = 0.4; %s
	freq = Fs; %Hz
	numberOfSamples = round(timeWindow*freq);
	
	N = size(emg,1);
	D = size(emg,2);
	emgCov = zeros(N,D,D);
	for i = numberOfSamples:size(emg,1) % the first samples are zero, but there are not used
	emgCov(i,:,:) = reshape( cov(emg(i-numberOfSamples+1:i,:))+1E-12*eye(D), 1,D,D);
	end
	
	% Extract data corresponding to activation index 1 (measurement time) 
	% for rest and the six movements -> training data
	activeIdx = find(stim(:,12)==1);

	tmgTrainVec = brac(activeIdx,:);
	emgTrain = emg(activeIdx,:);
	emgCovTrain = emgCov(activeIdx,:,:);
	yTrain = stim(activeIdx,1:9);
	
	multTrainVec = [tmgTrainVec,emgTrain];
	
	% Training input TMG data in matrix form (tactile bracelet of 8 x 40 sensors)
	N = size(tmgTrainVec,1);
	tmgTrain = zeros(N,8,40);
	for i=1:size(tmgTrainVec,1)
		tmgTrain(i,:,:) = reshape(tmgTrainVec(i,:),8,40);
	end
	
	% Extract data corresponding to intermediate results with activation
	% index 0 -> testing data
	activeIdx = find(~mod(stim(:,10),1) & ~mod(stim(:,11),1) & stim(:,12)==0);

	tmgTestVec = brac(activeIdx,:);
	emgTest = emg(activeIdx,:);
	emgCovTest = emgCov(activeIdx,:,:);
	yTest = stim(activeIdx,:);
	
	multTestVec = [tmgTestVec,emgTest];
	
	% Testing input TMG data in matrix form
	N = size(tmgTestVec,1);
	tmgTest = zeros(N,8,40);
	for i=1:size(tmgTestVec,1)
		tmgTest(i,:,:) = reshape(tmgTestVec(i,:),8,40);
	end
	
	disp(' -> Done!');
	
	%% Compute distances for EMG
	fprintf('Compute training kernels for EMG... ');
	kerdistTrTr_emg_SPD = covmatCovDist(emgCovTrain,emgCovTrain,0);
	kerdistTrTr_emg_logSPD = covmatCovDist(emgCovTrain,emgCovTrain,1);
	kerdistTrTr_emg_Eucl = compute_Eucldist(emgTrain,emgTrain);
	fprintf('Done.\n');
	
	fprintf('Compute testing kernels for EMG... ');
	kerdistTeTr_emg_SPD = covmatCovDist(emgCovTest,emgCovTrain,0);
	kerdistTeTr_emg_logSPD = covmatCovDist(emgCovTest,emgCovTrain,1);
	kerdistTeTr_emg_Eucl = compute_Eucldist(emgTest,emgTrain);
	fprintf('Done.\n');
	
	%% Compute distances for TMG
	fprintf('Compute training kernel for TMG (KL)... ');
	tic
	kerdistTrTr_tmg_KL = tensorKLdivGauss(tmgTrain,tmgTrain);
	fprintf('Done in %.2fs.\n',toc);
	
	fprintf('Compute testing kernel for TMG (KL)... ');
	tic
	kerdistTeTr_tmg_KL = tensorKLdivGauss(tmgTest,tmgTrain);
	fprintf('Done in %.2fs.\n',toc);
	
	fprintf('Compute training kernel for TMG (SPD)... ');
	tic
	kerdistTrTr_tmg_SPD = tensorCovDist(tmgTrain,tmgTrain,0);
	fprintf('Done in %.2fs.\n',toc);
	
	fprintf('Compute testing kernel for TMG (SPD)... ');
	tic
	kerdistTeTr_tmg_SPD = tensorCovDist(tmgTest,tmgTrain,0);
	fprintf('Done in %.2fs.\n',toc);
	
	fprintf('Compute training kernel for TMG (KL)... ');
	tic
	kerdistTrTr_tmg_Eucl = compute_Eucldist(tmgTrain,tmgTrain);
	fprintf('Done in %.2fs.\n',toc);
	
	fprintf('Compute testing kernel for TMG (KL)... ');
	tic
	kerdistTeTr_tmg_Eucl = compute_Eucldist(tmgTest,tmgTrain);
	fprintf('Done in %.2fs.\n',toc);
	

	%% Tensor Gaussian process TMG/EMG: KL divergence/SPD
	disp('Tensor Gaussian Process TMG (KL) - EMG (SPD)');
	
	kerdistTrTr = cat(3, kerdistTrTr_tmg_KL, kerdistTrTr_emg_SPD); 
	
	fprintf('Cross validation to find best parameters...');
	betaParam_KL_SPD{1} = logspace(-1,3,10);
	betaParam_KL_SPD{2} = logspace(-1,3,10);
	betaParam_KL_SPD{3} = logspace(-2,2,10);
	predictor_KL_SPD = @(XTr, YTr, XTe, idxTr, idxTe, param) tensorGPR_cv(XTr, YTr, XTe, idxTr, idxTe, param, kerdistTrTr);
	[ ~, teRMSE ] = crossvalidation( tmgTrain, yTrain, nbIterCV, testProp, predictor_KL_SPD, betaParam_KL_SPD);
	
	meanRMSE_KL_SPD = mean(teRMSE,ndims(teRMSE));
	[minTeRMSE_KL_SPD, minId] = min(meanRMSE_KL_SPD(:));
	[idTeRMSE_KL_SPD(1), idTeRMSE_KL_SPD(2), idTeRMSE_KL_SPD(3)] = ind2sub(size(meanRMSE_KL_SPD),minId);
	
	beta_KL_SPD(1,1,1) = betaParam_KL_SPD{1}(idTeRMSE_KL_SPD(1));
	beta_KL_SPD(1,1,2) = betaParam_KL_SPD{2}(idTeRMSE_KL_SPD(2));
	beta_KL_SPD(1,1,3) = betaParam_KL_SPD{3}(idTeRMSE_KL_SPD(3));
	
	fprintf('Min. training RMSE GP KL-SPD: %f for beta1 %.2f, beta2 %.2f and beta3 %.2f\n', minTeRMSE_KL_SPD, beta_KL_SPD(1,1,1), beta_KL_SPD(1,1,2), beta_KL_SPD(1,1,3));
	
	fprintf('Testing on intermediate data... ')
	tic
	kerdistTeTr = cat(3, kerdistTeTr_tmg_KL, kerdistTeTr_emg_SPD);
	yhat_tmgemg_KL_SPD = tensorGPR(tmgTrain, yTrain, tmgTest, beta_KL_SPD, kerdistTrTr, kerdistTeTr);
	fprintf('Done.\n');
	
	%% Tensor Gaussian process TMG/EMG: KL divergence/logSPD
	disp('Gaussian process TMG (KL) - EMG (logSPD)');
	
	kerdistTrTr = cat(3, kerdistTrTr_tmg_KL, kerdistTrTr_emg_logSPD); 
	
	fprintf('Cross validation to find best parameters...');
	betaParam_KL_logSPD{1} = logspace(-1,3,10);
	betaParam_KL_logSPD{2} = logspace(-1,3,10);
	betaParam_KL_logSPD{3} = logspace(-2,2,10);
	predictor_KL_logSPD = @(XTr, YTr, XTe, idxTr, idxTe, param) tensorGPR_cv(XTr, YTr, XTe, idxTr, idxTe, param, kerdistTrTr);
	[ ~, teRMSE ] = crossvalidation( tmgTrain, yTrain, nbIterCV, testProp, predictor_KL_logSPD, betaParam_KL_logSPD);
	
	meanRMSE_KL_logSPD = mean(teRMSE,ndims(teRMSE));
	[minTeRMSE_KL_logSPD, minId] = min(meanRMSE_KL_logSPD(:));
	[idTeRMSE_KL_logSPD(1), idTeRMSE_KL_logSPD(2), idTeRMSE_KL_logSPD(3)] = ind2sub(size(meanRMSE_KL_logSPD),minId);
	
	beta_KL_logSPD(1,1,1) = betaParam_KL_logSPD{1}(idTeRMSE_KL_logSPD(1));
	beta_KL_logSPD(1,1,2) = betaParam_KL_logSPD{2}(idTeRMSE_KL_logSPD(2));
	beta_KL_logSPD(1,1,3) = betaParam_KL_logSPD{3}(idTeRMSE_KL_logSPD(3));
	
	fprintf('Min. training RMSE GP KL-logSPD: %f for beta1 %.2f, beta2 %.2f and beta3 %.2f\n', minTeRMSE_KL_logSPD, beta_KL_logSPD(1,1,1), beta_KL_logSPD(1,1,2), beta_KL_logSPD(1,1,3));
	
	fprintf('Testing on intermediate data... ')
	tic
	kerdistTeTr = cat(3, kerdistTeTr_tmg_KL, kerdistTeTr_emg_logSPD);
	yhat_tmgemg_KL_logSPD = tensorGPR(tmgTrain, yTrain, tmgTest, beta_KL_logSPD, kerdistTrTr, kerdistTeTr);
	fprintf('Done.\n');

	%% Tensor Gaussian process TMG/EMG: KL divergence/Eucl
	disp('Gaussian process TMG (KL) - EMG (Eucl)');
	
	kerdistTrTr = cat(3, kerdistTrTr_tmg_KL, kerdistTrTr_emg_Eucl); 
	
	fprintf('Cross validation to find best parameters...');
	betaParam_KL_Eucl{1} = logspace(-1,3,10);
	betaParam_KL_Eucl{2} = logspace(-1,3,10);
	betaParam_KL_Eucl{3} = logspace(-1,1,10);
	predictor_KL_Eucl = @(XTr, YTr, XTe, idxTr, idxTe, param) tensorGPR_cv(XTr, YTr, XTe, idxTr, idxTe, param, kerdistTrTr);
	[ ~, teRMSE ] = crossvalidation( tmgTrain, yTrain, nbIterCV, testProp, predictor_KL_Eucl, betaParam_KL_Eucl);
	
	meanRMSE_KL_Eucl = mean(teRMSE,ndims(teRMSE));
	[minTeRMSE_KL_Eucl, minId] = min(meanRMSE_KL_Eucl(:));
	[idTeRMSE_KL_Eucl(1), idTeRMSE_KL_Eucl(2), idTeRMSE_KL_Eucl(3)] = ind2sub(size(meanRMSE_KL_Eucl),minId);
	
	beta_KL_Eucl(1,1,1) = betaParam_KL_Eucl{1}(idTeRMSE_KL_Eucl(1));
	beta_KL_Eucl(1,1,2) = betaParam_KL_Eucl{2}(idTeRMSE_KL_Eucl(2));
	beta_KL_Eucl(1,1,3) = betaParam_KL_Eucl{3}(idTeRMSE_KL_Eucl(3));
	
	fprintf('Min. training RMSE GP KL-Eucl: %f for beta1 %.2f, beta2 %.2f and beta3 %.2f\n', minTeRMSE_KL_Eucl, beta_KL_Eucl(1,1,1), beta_KL_Eucl(1,1,2), beta_KL_Eucl(1,1,3));
	
	fprintf('Testing on intermediate data... ')
	tic
	kerdistTeTr = cat(3, kerdistTeTr_tmg_KL, kerdistTeTr_emg_Eucl);
	yhat_tmgemg_KL_Eucl = tensorGPR(tmgTrain, yTrain, tmgTest, beta_KL_Eucl, kerdistTrTr, kerdistTeTr);
	fprintf('Done.\n');

	%% Tensor Gaussian process TMG/EMG: SPD/SPD
	disp('Gaussian process TMG (SPD) - EMG (SPD)');
	
	kerdistTrTr = cat(3, kerdistTrTr_tmg_SPD, kerdistTrTr_emg_SPD); 
	
	fprintf('Cross validation to find best parameters...');
	betaParam_SPD_SPD{1} = logspace(-1,1.5,10);
	betaParam_SPD_SPD{2} = logspace(-1,1.5,10);
	betaParam_SPD_SPD{3} = logspace(-2,2,10);
	predictor_SPD_SPD = @(XTr, YTr, XTe, idxTr, idxTe, param) tensorGPR_cv(XTr, YTr, XTe, idxTr, idxTe, param, kerdistTrTr);
	[ ~, teRMSE ] = crossvalidation( tmgTrain, yTrain, nbIterCV, testProp, predictor_SPD_SPD, betaParam_SPD_SPD);
	
	meanRMSE_SPD_SPD = mean(teRMSE,ndims(teRMSE));
	[minTeRMSE_SPD_SPD, minId] = min(meanRMSE_SPD_SPD(:));
	[idTeRMSE_SPD_SPD(1), idTeRMSE_SPD_SPD(2), idTeRMSE_SPD_SPD(3)] = ind2sub(size(meanRMSE_SPD_SPD),minId);
	
	beta_SPD_SPD(1,1,1) = betaParam_SPD_SPD{1}(idTeRMSE_SPD_SPD(1));
	beta_SPD_SPD(1,1,2) = betaParam_SPD_SPD{2}(idTeRMSE_SPD_SPD(2));
	beta_SPD_SPD(1,1,3) = betaParam_SPD_SPD{3}(idTeRMSE_SPD_SPD(3));
	
	fprintf('Min. training RMSE GP SPD-SPD: %f for beta1 %.2f, beta2 %.2f and beta3 %.2f\n', minTeRMSE_SPD_SPD, beta_SPD_SPD(1,1,1), beta_SPD_SPD(1,1,2), beta_SPD_SPD(1,1,3));
	
	fprintf('Testing on intermediate data... ')
	tic
	kerdistTeTr = cat(3, kerdistTeTr_tmg_SPD, kerdistTeTr_emg_SPD);
	yhat_tmgemg_SPD_SPD = tensorGPR(tmgTrain, yTrain, tmgTest, beta_SPD_SPD, kerdistTrTr, kerdistTeTr);
	fprintf('Done.\n');

	%% Tensor Gaussian process TMG/EMG: SPD/logSPD
	disp('Gaussian process TMG (SPD) - EMG (logSPD)');
	
	kerdistTrTr = cat(3, kerdistTrTr_tmg_SPD, kerdistTrTr_emg_logSPD); 
	
	fprintf('Cross validation to find best parameters...');
	betaParam_SPD_logSPD{1} = logspace(-1,1.5,10);
	betaParam_SPD_logSPD{2} = logspace(-1,1.5,10);
	betaParam_SPD_logSPD{3} = logspace(-2,2,10);
	predictor_SPD_logSPD = @(XTr, YTr, XTe, idxTr, idxTe, param) tensorGPR_cv(XTr, YTr, XTe, idxTr, idxTe, param, kerdistTrTr);
	[ ~, teRMSE ] = crossvalidation( tmgTrain, yTrain, nbIterCV, testProp, predictor_SPD_logSPD, betaParam_SPD_logSPD);
	
	meanRMSE_SPD_logSPD = mean(teRMSE,ndims(teRMSE));
	[minTeRMSE_SPD_logSPD, minId] = min(meanRMSE_SPD_logSPD(:));
	[idTeRMSE_SPD_logSPD(1), idTeRMSE_SPD_logSPD(2), idTeRMSE_SPD_logSPD(3)] = ind2sub(size(meanRMSE_SPD_logSPD),minId);
	
	beta_SPD_logSPD(1,1,1) = betaParam_SPD_logSPD{1}(idTeRMSE_SPD_logSPD(1));
	beta_SPD_logSPD(1,1,2) = betaParam_SPD_logSPD{2}(idTeRMSE_SPD_logSPD(2));
	beta_SPD_logSPD(1,1,3) = betaParam_SPD_logSPD{2}(idTeRMSE_SPD_logSPD(3));
	
	fprintf('Min. training RMSE GP SPD-logSPD: %f for beta1 %.2f, beta2 %.2f and beta3 %.2f\n', minTeRMSE_SPD_logSPD, beta_SPD_logSPD(1,1,1), beta_SPD_logSPD(1,1,2), beta_SPD_logSPD(1,1,3));
	
	fprintf('Testing on intermediate data... ')
	tic
	kerdistTeTr = cat(3, kerdistTeTr_tmg_SPD, kerdistTeTr_emg_logSPD);
	yhat_tmgemg_SPD_logSPD = tensorGPR(tmgTrain, yTrain, tmgTest, beta_SPD_logSPD, kerdistTrTr, kerdistTeTr);
	fprintf('Done.\n');
	
	%% Tensor Gaussian process TMG/EMG: SPD/Eucl
	disp('Gaussian process TMG (SPD) - EMG (Eucl)');
	
	kerdistTrTr = cat(3, kerdistTrTr_tmg_SPD, kerdistTrTr_emg_Eucl); 
	
	fprintf('Cross validation to find best parameters...');
	betaParam_SPD_Eucl{1} = logspace(-1,1.5,10);
	betaParam_SPD_Eucl{2} = logspace(-1,1.5,10);
	betaParam_SPD_Eucl{3} = logspace(-1,1,10);
	predictor_Eucl_Eucl = @(XTr, YTr, XTe, idxTr, idxTe, param) tensorGPR_cv(XTr, YTr, XTe, idxTr, idxTe, param, kerdistTrTr);
	[ ~, teRMSE ] = crossvalidation( tmgTrain, yTrain, nbIterCV, testProp, predictor_Eucl_Eucl, betaParam_SPD_Eucl);
	
	meanRMSE_SPD_Eucl = mean(teRMSE,ndims(teRMSE));
	[minTeRMSE_SPD_Eucl, minId] = min(meanRMSE_SPD_Eucl(:));
	[idTeRMSE_SPD_Eucl(1), idTeRMSE_SPD_Eucl(2), idTeRMSE_SPD_Eucl(3)] = ind2sub(size(meanRMSE_SPD_Eucl),minId);
	
	beta_SPD_Eucl(1,1,1) = betaParam_SPD_Eucl{1}(idTeRMSE_SPD_Eucl(1));
	beta_SPD_Eucl(1,1,2) = betaParam_SPD_Eucl{2}(idTeRMSE_SPD_Eucl(2));
	beta_SPD_Eucl(1,1,3) = betaParam_SPD_Eucl{3}(idTeRMSE_SPD_Eucl(3));
	
	fprintf('Min. training RMSE GP SPD-Eucl: %f for beta1 %.2f, beta2 %.2f and beta3 %.2f\n', minTeRMSE_SPD_Eucl, beta_SPD_Eucl(1,1,1), beta_SPD_Eucl(1,1,2), beta_SPD_Eucl(1,1,3));
	
	fprintf('Testing on intermediate data... ')
	tic
	kerdistTeTr = cat(3, kerdistTeTr_tmg_SPD, kerdistTeTr_emg_Eucl);
	yhat_tmgemg_SPD_Eucl = tensorGPR(tmgTrain, yTrain, tmgTest, beta_SPD_Eucl, kerdistTrTr, kerdistTeTr);
	fprintf('Done.\n');
	
	%% Tensor Gaussian process TMG/EMG: Eucl/SPD
	disp('Gaussian process TMG (Eucl) - EMG (SPD)');
	
	kerdistTrTr = cat(3, kerdistTrTr_tmg_Eucl, kerdistTrTr_emg_SPD); 
	
	fprintf('Cross validation to find best parameters...');
	betaParam_Eucl_SPD{1} = logspace(-1,1,10);
	betaParam_Eucl_SPD{2} = logspace(-2,2,10);
	predictor_Eucl_SPD = @(XTr, YTr, XTe, idxTr, idxTe, param) tensorGPR_cv(XTr, YTr, XTe, idxTr, idxTe, param, kerdistTrTr);
	[ ~, teRMSE ] = crossvalidation( tmgTrain, yTrain, nbIterCV, testProp, predictor_Eucl_SPD, betaParam_Eucl_SPD);
	
	meanRMSE_Eucl_SPD = mean(teRMSE,ndims(teRMSE));
	[minTeRMSE_Eucl_SPD, minId] = min(meanRMSE_Eucl_SPD(:));
	[idTeRMSE_Eucl_SPD(1), idTeRMSE_Eucl_SPD(2)] = ind2sub(size(meanRMSE_Eucl_SPD),minId);
	
	beta_Eucl_SPD(1,1,1) = betaParam_Eucl_SPD{1}(idTeRMSE_Eucl_SPD(1));
	beta_Eucl_SPD(1,1,2) = betaParam_Eucl_SPD{2}(idTeRMSE_Eucl_SPD(2));
	
	fprintf('Min. training RMSE GP Eucl-SPD: %f for beta1 %.2f and beta2 %.2f\n', minTeRMSE_Eucl_SPD, beta_Eucl_SPD(1,1,1), beta_Eucl_SPD(1,1,2));
	
	fprintf('Testing on intermediate data... ')
	tic
	kerdistTeTr = cat(3, kerdistTeTr_tmg_Eucl, kerdistTeTr_emg_SPD);
	yhat_tmgemg_Eucl_SPD = tensorGPR(tmgTrain, yTrain, tmgTest, beta_Eucl_SPD, kerdistTrTr, kerdistTeTr);
	fprintf('Done.\n');
	
	%% Tensor Gaussian process TMG/EMG: Eucl/logSPD
	disp('Gaussian process TMG (Eucl) - EMG (logSPD)');
	
	kerdistTrTr = cat(3, kerdistTrTr_tmg_Eucl, kerdistTrTr_emg_logSPD); 
	
	fprintf('Cross validation to find best parameters...');
	betaParam_Eucl_logSPD{1} = logspace(-1,1,10);
	betaParam_Eucl_logSPD{2} = logspace(-2,2,10);
	predictor_Eucl_logSPD = @(XTr, YTr, XTe, idxTr, idxTe, param) tensorGPR_cv(XTr, YTr, XTe, idxTr, idxTe, param, kerdistTrTr);
	[ ~, teRMSE ] = crossvalidation( tmgTrain, yTrain, nbIterCV, testProp, predictor_Eucl_logSPD, betaParam_Eucl_logSPD);
	
	meanRMSE_Eucl_logSPD = mean(teRMSE,ndims(teRMSE));
	[minTeRMSE_Eucl_logSPD, minId] = min(meanRMSE_Eucl_logSPD(:));
	[idTeRMSE_Eucl_logSPD(1), idTeRMSE_Eucl_logSPD(2)] = ind2sub(size(meanRMSE_Eucl_logSPD),minId);
	
	beta_Eucl_logSPD(1,1,1) = betaParam_Eucl_logSPD{1}(idTeRMSE_Eucl_logSPD(1));
	beta_Eucl_logSPD(1,1,2) = betaParam_Eucl_logSPD{2}(idTeRMSE_Eucl_logSPD(2));
	
	fprintf('Min. training RMSE GP Eucl-logSPD: %f for beta1 %.2f and beta2 %.2f\n', minTeRMSE_Eucl_logSPD, beta_Eucl_logSPD(1,1,1), beta_Eucl_logSPD(1,1,2));
	
	fprintf('Testing on intermediate data... ')
	tic
	kerdistTeTr = cat(3, kerdistTeTr_tmg_Eucl, kerdistTeTr_emg_logSPD);
	yhat_tmgemg_Eucl_logSPD = tensorGPR(tmgTrain, yTrain, tmgTest, beta_Eucl_logSPD, kerdistTrTr, kerdistTeTr);
	fprintf('Done.\n');
	
	%% Gaussian process TMG/EMG: Eucl/Eucl
	disp('Tensor Gaussian Process TMG (Eucl) - EMG (Eucl)');
	
	kerdistTrTr = cat(3, kerdistTrTr_tmg_Eucl, kerdistTrTr_emg_Eucl); 
	
	fprintf('Cross validation to find best parameters...');
	betaParam_Eucl_Eucl{1} = logspace(-1,1,10);
	betaParam_Eucl_Eucl{2} = logspace(-1,1,10);
	predictor_Eucl_Eucl = @(XTr, YTr, XTe, idxTr, idxTe, param) tensorGPR_cv(XTr, YTr, XTe, idxTr, idxTe, param, kerdistTrTr);
	[ ~, teRMSE ] = crossvalidation( tmgTrain, yTrain, nbIterCV, testProp, predictor_Eucl_Eucl, betaParam_Eucl_Eucl);
	
	meanRMSE_Eucl_Eucl = mean(teRMSE,ndims(teRMSE));
	[minTeRMSE_Eucl_Eucl, minId] = min(meanRMSE_Eucl_Eucl(:));
	[idTeRMSE_Eucl_Eucl(1), idTeRMSE_Eucl_Eucl(2)] = ind2sub(size(meanRMSE_Eucl_Eucl),minId);
	
	beta_Eucl_Eucl(1,1,1) = betaParam_Eucl_Eucl{1}(idTeRMSE_Eucl_Eucl(1));
	beta_Eucl_Eucl(1,1,2) = betaParam_Eucl_Eucl{2}(idTeRMSE_Eucl_Eucl(2));
	
	fprintf('Min. training RMSE GP Eucl-Eucl: %f for beta1 %.2f and beta2 %.2f\n', minTeRMSE_Eucl_Eucl, beta_Eucl_Eucl(1,1,1), beta_Eucl_Eucl(1,1,2));
	
	fprintf('Testing on intermediate data... ')
	tic
	kerdistTeTr = cat(3, kerdistTeTr_tmg_Eucl, kerdistTeTr_emg_Eucl);
	yhat_tmgemg_Eucl_Eucl = tensorGPR(tmgTrain, yTrain, tmgTest, beta_Eucl_Eucl, kerdistTrTr, kerdistTeTr);
	fprintf('Done.\n');
	
	%% Ridge Regression
	disp('Ridge Regression');
	
	lambdaParam = logspace(-4,-1,20);
	predictorRR = @(XTr, YTr, XTe, idxTr, idxTe, param) RR(XTr, YTr, XTe, idxTr, idxTe, param);
	[ ~, teRMSE ] = crossvalidation( multTrainVec, yTrain, nbIterCV, testProp, predictorRR, lambdaParam);
	
	meanRMSE_RR = mean(teRMSE,ndims(teRMSE));
	
	[minTeRMSE_RR, idTeRMSE_RR] = min(meanRMSE_RR(:));
	
	lambda = lambdaParam(idTeRMSE_RR);
	
	fprintf('Min. training RMSE RR %f for lambda %f\n', minTeRMSE_RR, lambda);
	
	fprintf('Testing on intermediate data... ')
	tic
	yhat_tmgemg_rr = RR(multTrainVec, yTrain, multTestVec, 1, 1, lambda);
	fprintf('Done in %.2fs.\n',toc);
	
	%% Save predictions
	save(predstr,'tmgTrain','tmgTest','emgCovTrain','emgCovTest','emgTrain','emgTest','yTrain','yTest',...
		'yhat_tmgemg_KL_SPD','yhat_tmgemg_KL_logSPD', 'yhat_tmgemg_KL_Eucl','yhat_tmgemg_SPD_SPD', ...
		'yhat_tmgemg_SPD_logSPD','yhat_tmgemg_SPD_Eucl', 'yhat_tmgemg_Eucl_SPD','yhat_tmgemg_Eucl_logSPD',...
		'yhat_tmgemg_Eucl_Eucl','yhat_tmgemg_rr');

	
	%% Compute errors
	err_KL_SPD = computeRMSE(yTest(:,1:9),yhat_tmgemg_KL_SPD);
	err_KL_logSPD = computeRMSE(yTest(:,1:9),yhat_tmgemg_KL_logSPD);
	err_KL_Eucl = computeRMSE(yTest(:,1:9),yhat_tmgemg_KL_Eucl);

	err_SPD_SPD = computeRMSE(yTest(:,1:9),yhat_tmgemg_SPD_SPD);
	err_SPD_logSPD = computeRMSE(yTest(:,1:9),yhat_tmgemg_SPD_logSPD);
	err_SPD_Eucl = computeRMSE(yTest(:,1:9),yhat_tmgemg_SPD_Eucl);

	err_Eucl_SPD = computeRMSE(yTest(:,1:9),yhat_tmgemg_Eucl_SPD);
	err_Eucl_logSPD = computeRMSE(yTest(:,1:9),yhat_tmgemg_Eucl_logSPD);
	err_Eucl_Eucl = computeRMSE(yTest(:,1:9),yhat_tmgemg_Eucl_Eucl);

	err_RR = computeRMSE(yTest(:,1:9),yhat_tmgemg_rr);
	%% Write in txt file
	fileID = fopen(filestr,'w');
		
	fprintf(fileID,'KL-SPD      %2.2f for beta1 %2f, beta2 %2f and beta3 %2f\n',err_KL_SPD, beta_KL_SPD(1,1,1), beta_KL_SPD(1,1,2), beta_KL_SPD(1,1,3));
	fprintf(fileID,'KL-logSPD   %2.2f for beta1 %2f, beta2 %2f and beta3 %2f\n',err_KL_logSPD, beta_KL_logSPD(1,1,1), beta_KL_logSPD(1,1,2), beta_KL_logSPD(1,1,3));
	fprintf(fileID,'KL-Eucl     %2.2f for beta1 %2f, beta2 %2f and beta3 %2f\n',err_KL_Eucl, beta_KL_Eucl(1,1,1), beta_KL_Eucl(1,1,2), beta_KL_Eucl(1,1,3));
	
	fprintf(fileID,'SPD-SPD     %2.2f for beta1 %2f, beta2 %2f and beta3 %2f\n',err_SPD_SPD, beta_SPD_SPD(1,1,1), beta_SPD_SPD(1,1,2), beta_SPD_SPD(1,1,3));
	fprintf(fileID,'SPD-logSPD  %2.2f for beta1 %2f, beta2 %2f and beta3 %2f\n',err_SPD_logSPD, beta_SPD_logSPD(1,1,1), beta_SPD_logSPD(1,1,2), beta_SPD_logSPD(1,1,3));
	fprintf(fileID,'SPD-Eucl    %2.2f for beta1 %2f, beta2 %2f and beta3 %2f\n',err_SPD_Eucl, beta_SPD_Eucl(1,1,1), beta_SPD_Eucl(1,1,2), beta_SPD_Eucl(1,1,3));
	
	fprintf(fileID,'Eucl-SPD    %2.2f for beta1 %2f and beta2 %2f\n', err_Eucl_SPD, beta_Eucl_SPD(1,1,1), beta_Eucl_SPD(1,1,2));
	fprintf(fileID,'Eucl-logSPD %2.2f for beta1 %2f and beta2 %2f\n', err_Eucl_logSPD, beta_Eucl_logSPD(1,1,1), beta_Eucl_logSPD(1,1,2));
	fprintf(fileID,'Eucl-Eucl   %2.2f for beta1 %2f and beta2 %2f\n', err_Eucl_Eucl, beta_Eucl_Eucl(1,1,1), beta_Eucl_Eucl(1,1,2));
	
	fprintf(fileID,'RR         %2.2f for lambda %f\n',err_RR,lambda);
	
	fclose(fileID);
end

end

function [ trRMSE, teRMSE ] = crossvalidation( X, Y, iterationsNb, testRatio, predict, param)
%crossvalidation Perform crossvalidation to return the train and test RMSE
%for the method predict and all values of parameters param (up to 4 dim.)
%
% X             Input N x D1 x .. x D_I 
% y             Output N x d
% iterationsNb  Number of random shuffles
% testRatio		Proportion of test data
% predict		Machine learning method used to find Y from X
% param			Values of the parameter(s) of the method predict
% trRMSE		Training RMSE for each parameter and iteration
% teRMSE		Testing RMSE for each parameter and iteration

randomStream = RandStream.getGlobalStream; % using globalStream

trRMSE = [];
teRMSE = [];

getDivision = divideDataSet( X, Y, iterationsNb, testRatio, randomStream );

for iteration = 1:iterationsNb
	if ~mod(iteration,10)
		fprintf('%d out of %d...\n',iteration, iterationsNb);
	end
	% Separe X and y in train and test parts
	[ XTr, XTe, YTr, YTe, idxTr, idxTe ] = getDivision(iteration);
	
	if ~iscell(param)
		for p = 1:length(param)
			% Make prediction for test data
			yTrPred = predict( XTr, YTr, XTr, idxTr, idxTr, param(p) );
			yTePred = predict( XTr, YTr, XTe, idxTr, idxTe, param(p) );
			
			% Compute RMSE
			trRMSE(p,iteration) = computeRMSE(YTr, yTrPred);
			teRMSE(p,iteration) = computeRMSE(YTe, yTePred);
		end
	else
		if numel(param) == 2
			for p1 = 1:length(param{1})
				for p2 = 1:length(param{2})
					p(1,1,1) = param{1}(p1);
					p(1,1,2) = param{2}(p2);
					% Make prediction for test data
					yTrPred = predict( XTr, YTr, XTr, idxTr, idxTr, p );
					yTePred = predict( XTr, YTr, XTe, idxTr, idxTe, p );
					
					% Compute RMSE
					trRMSE(p1,p2,iteration) = computeRMSE(YTr, yTrPred);
					teRMSE(p1,p2,iteration) = computeRMSE(YTe, yTePred);
				end
			end
		else
			if numel(param) == 3
				for p1 = 1:length(param{1})
					for p2 = 1:length(param{2})
						for p3 = 1:length(param{3})
						p(1,1,1) = param{1}(p1);
						p(1,1,2) = param{2}(p2);
						p(1,1,3) = param{2}(p3);
						% Make prediction for test data
						yTrPred = predict( XTr, YTr, XTr, idxTr, idxTr, p );
						yTePred = predict( XTr, YTr, XTe, idxTr, idxTe, p );
						
						% Compute RMSE
						trRMSE(p1,p2,p3,iteration) = computeRMSE(YTr, yTrPred);
						teRMSE(p1,p2,p3,iteration) = computeRMSE(YTe, yTePred);
					end
					end
				end
			else
				if numel(param) == 4
					for p1 = 1:length(param{1})
						for p2 = 1:length(param{2})
							for p3 = 1:length(param{3})
								for p4 = 1:length(param{4})
									p(1,1,1) = param{1}(p1);
									p(1,1,2) = param{2}(p2);
									p(1,1,3) = param{2}(p3);
									p(1,1,4) = param{2}(p4);
									% Make prediction for test data
									yTrPred = predict( XTr, YTr, XTr, idxTr, idxTr, p );
									yTePred = predict( XTr, YTr, XTe, idxTr, idxTe, p );

									% Compute RMSE
									trRMSE(p1,p2,p3,p4,iteration) = computeRMSE(YTr, yTrPred);
									teRMSE(p1,p2,p3,p4,iteration) = computeRMSE(YTe, yTePred);
								end
							end
						end
					end
					
				end
			end
		end
	end
end
end

function [ getDivision ] = divideDataSet( X, Y, iterationsNb, testRatio, randomStream )
%divideDataSet Divide the dataset in train and test parts
%
% X             Input N x D1 x .. x D_I 
% y             Output N x d
% iterationsNb  Number of random shuffles
% randomStream  Random stream used to divide the dataset
% getDivision   Function from iteration number to [XTr, XTe, yTr, yTe]

dataSize = size(Y, 1);

trainSize = floor(dataSize*(1-testRatio));
idxPerIt = zeros(iterationsNb, dataSize);

for iteration = 1:iterationsNb
    idxPerIt(iteration, :) = randomStream.randperm(dataSize);
end


    function [XTr, XTe, YTr, YTe, idxTr, idxTe] = i_getDivision(iterationNumber)
        % Get first part of shuffled samples for train, the rest for test
        idxTe = idxPerIt(iterationNumber, trainSize+1:end);
        idxTr = idxPerIt(iterationNumber, 1:trainSize);
        
        dim = ndims(X);
        indsTe = repmat({':'},1,dim);
        indsTe{1} = idxTe;
        YTe = Y(idxTe, :);
        XTe = X(indsTe{:});
        
        indsTr = repmat({':'},1,dim);
        indsTr{1} = idxTr;
        YTr = Y(idxTr, :);
        XTr = X(indsTr{:});
        
    end

getDivision = @i_getDivision;

end

function [ error ] = computeRMSE( Y, Yhat )
%computeRMSE Compute the root mean square error
%
% Y       Expected output N x d
% YHat    Predicted output N x d

N = length(Y);

diff_square = sum((Y-Yhat).^2,2);

error = sum(diff_square);

error = sqrt(error./N);
end

function Yhat = tensorGPR_cv(XTr, YTr, XTe, idxTr, idxTe, beta, kerdist)
%tensorGPR_cv Predict a new output with tensor Gaussian process regression 
%(used for cross-validation only)
%
% XTr		Input N x D1 x .. x D_I
% YTr		Output N x d
% XTe		New input Nte x D1 x .. x D_I
% idxTr		Training indexes
% idxTe		Testing indexes
% beta		TGP parameter
% kerdist	Kernel distance computed for all inputs (N+Nte) x (N+Nte) x I
% Yhat		Predicted output Nte x d

diagRegFact = 1E-4;

M = ndims(kerdist);

kerdistTrTr = kerdist(idxTr,idxTr,:);
kerdistTeTr = kerdist(idxTe,idxTr,:);

kerdistTrTr = kerdistTrTr./repmat((2.*beta.^2).*ones(1,1,M+1-ndims(beta)),size(XTr,1),size(XTr,1),1);
KTrTr = prod(exp(-kerdistTrTr),3);

kerdistTeTr = kerdistTeTr./repmat((2.*beta.^2).*ones(1,1,M+1-ndims(beta)),size(XTe,1),size(XTr,1),1);
KTeTr = prod(exp(-kerdistTeTr),3);

Yhat = KTeTr*((KTrTr+diagRegFact.*eye(size(XTr,1)))\YTr);
end

function Yhat = tensorGPR(XTr, YTr, XTe, beta, kerdistTrTr, kerdistTeTr)
%tensorGPR Predict a new output with tensor Gaussian process regression 
%
% XTr			Input
% YTr			Output
% XTe			New input
% beta			TGP parameter
% kerdistTrTr	Kernel distance between input data
% kerdistTeTr	Kernel distance between new input and input data
% Yhat			Predicted output

diagRegFact = 1E-4;

M = ndims(kerdistTrTr);

kerdistTrTr = kerdistTrTr./repmat((2.*beta.^2).*ones(1,1,M+1-ndims(beta)),size(XTr,1),size(XTr,1),1);
KTrTr = prod(exp(-kerdistTrTr),3);

kerdistTeTr = kerdistTeTr./repmat((2.*beta.^2).*ones(1,1,M+1-ndims(beta)),size(XTe,1),size(XTr,1),1);
KTeTr = prod(exp(-kerdistTeTr),3);

Yhat = KTeTr*((KTrTr+diagRegFact.*eye(size(XTr,1)))\YTr);
end

function Yhat = RR(XTr, YTr, XTe, idxTr, idxTe, lambda)
%RR Predict a new output with Ridge regression
%
% XTr		Input N x D
% YTr		Output N x d
% XTe		New input Nte x D
% idxTr		Training indexes (for cross-validation)
% idxTe		Testing indexes (for cross-validation)
% lambda	RR parameter
% Yhat		Predicted output Nte x d

% Define lambda matrix
lambdaMatrix = lambda*eye(size(XTr,2));

% Compute beta
beta = (XTr'*XTr + lambdaMatrix)\XTr'*YTr;

% Prediction
Yhat = XTe*beta;
end

function KLdivTens = tensorKLdivGauss(X1,X2)
%tensorKLdivGauss Computes the KL divergence KLdivTens between X1 and X2 by
%treating each matrix data as a generative model

N1 = size(X1,1);
N2 = size(X2,1);

M1 = ndims(X1);
M2 = ndims(X2);
if M1~=M2
	error('Inputs must have the same number of dimensions');
end
M = M1;

% Compute mean and covariance for each data
diagRegFact = 1E-4;
m1 = cell(N1,M-1);
S1 = cell(N1,M-1);
m2 = cell(N2,M-1);
S2 = cell(N2,M-1);

for n = 1:N1
	for m = 2:M
		mat = tens2mat(X1(n,:,:),m,[1:m-1,m+1:M]);
		m1{n,m-1} = mean(mat,2);
		S1{n,m-1} = cov(mat')+diagRegFact.*eye(size(m1{n,m-1},1));
	end
end

for n = 1:N2
	for m = 2:M
		mat = tens2mat(X2(n,:,:),m,[1:m-1,m+1:M]);
		m2{n,m-1} = mean(mat,2);
		S2{n,m-1} = cov(mat')+diagRegFact.*eye(size(m2{n,m-1},1));
	end
end

KLdivTens = zeros(N1,N2,M-1);

symKernel = false;
n2init = 1;
% Check if the kernel is symmetric (if yes, we compute the similar elements
% only once.
if N1 == N2
	if X1 == X2
		symKernel = true;
	end
end

% Compute the kernel
for n1 = 1:N1
	if ~mod(n1,100)
		fprintf('.');
	end
	if symKernel
		n2init = n1;
	end
	for n2 = n2init:N2
		for m = 2:M
			KLdivTens(n1,n2,m-1) = computeKLdivGauss(m1{n1,m-1},S1{n1,m-1},m2{n2,m-1},S2{n2,m-1});
			if symKernel
				KLdivTens(n2,n1,m-1) = KLdivTens(n1,n2,m-1);
			end
		end
	end
end
end

function KLdiv = computeKLdivGauss(m1,S1,m2,S2)
%KLdiv Compute the Kullback-Leibler divergence KLdiv between distributions
%N(m1,S1) and N(m2,S2)

k = size(m1,1);
KLdiv = 0.5*(trace(S1\S2)+((m1-m2)'/S1)*(m1-m2)-k+log(det(S1)/det(S2)));
end

function kerdist = compute_Eucldist(X1,X2)
%kerdist Compute the Euclidean distance kerdist between X1 and X2

N1 = size(X1,1);
N2 = size(X2,1);

symKernel = false;
n2init = 1;
% Check if the kernel is symmetric (if yes, we compute the similar elements
% only once.
if N1 == N2
	if X1 == X2
		symKernel = true;
	end
end

% Compute the kernel
kerdist = zeros(N1,N2);
for n1 = 1:N1
	if ~mod(n1,100)
		fprintf('.');
	end
	if symKernel
		n2init = n1;
	end
	for n2 = n2init:N2
		kerdist(n1,n2) = norm(X1(n1,:)-X2(n2,:),'fro');
		if symKernel
			kerdist(n2,n1) = kerdist(n1,n2);
		end
	end
end

end

function covDistTens = tensorCovDist(X1,X2,logDist)
%tensorCovDist Compute the covariance or log covariance (if logDist = 1)
%distance covDistTens between X1 and X2 (data are not SPD matrices)

N1 = size(X1,1);
N2 = size(X2,1);

M1 = ndims(X1);
M2 = ndims(X2);
if M1~=M2
	error('Inputs must have the same number of dimensions');
end
M = M1;

diagRegFact = 1E-4;

S1 = cell(N1,M-1);
S2 = cell(N2,M-1);

for n = 1:N1
	for m = 2:M
		mat = tens2mat(X1(n,:,:),m,[1:m-1,m+1:M]);
		S1{n,m-1} = cov(mat')+diagRegFact.*eye(size(mat,1));
	end
end

for n = 1:N2
	for m = 2:M
		mat = tens2mat(X2(n,:,:),m,[1:m-1,m+1:M]);
		S2{n,m-1} = cov(mat')+diagRegFact.*eye(size(mat,1));
	end
end

covDistTens = zeros(N1,N2,M-1);

symKernel = false;
n2init = 1;
% Check if the kernel is symmetric (if yes, we compute the similar elements
% only once.
if N1 == N2
	if X1 == X2
		symKernel = true;
	end
end

% Compute the kernel
for n1 = 1:N1
	if ~mod(n1,100)
		fprintf('.');
	end
	if symKernel
		n2init = n1;
	end
	for n2 = n2init:N2
		for m = 2:M
			if(logDist)
				covDistTens(n1,n2,m-1) = computeLogCovDist(S1{n1,m-1},S2{n2,m-1});
			else
				covDistTens(n1,n2,m-1) = computeCovDist(S1{n1,m-1},S2{n2,m-1});
			end
			if symKernel
				covDistTens(n2,n1,m-1) = covDistTens(n1,n2,m-1);
			end
		end
	end
end
end

function covDistTens = covmatCovDist(X1,X2,logDist)
%tensorCovDist Compute the covariance or log covariance (if logDist = 1)
%distance covDistTens between X1 and X2 (data are SPD matrices)

N1 = size(X1,1);
N2 = size(X2,1);

M1 = size(X1,2);
M2 = size(X1,3);

M3 = size(X2,2);
M4 = size(X2,3);

if M1~=M2 || M3~=M4 || M1~=M3
	error('Inputs must have the same size and be covariance matrices.');
end
M = M1;

diagRegFact = 1E-6;

covDistTens = zeros(N1,N2);

symKernel = false;
n2init = 1;
% Check if the kernel is symmetric (if yes, we compute the similar elements
% only once.
if N1 == N2
	if X1 == X2
		symKernel = true;
	end
end

% Compute the kernel
for n1 = 1:N1
	if ~mod(n1,100)
		fprintf('.');
	end
	if symKernel
		n2init = n1;
	end
	for n2 = n2init:N2
		S1 = permute(X1(n1,:,:),[2,3,1]) + diagRegFact.*eye(M);
		S2 = permute(X2(n2,:,:),[2,3,1]) + diagRegFact.*eye(M);
			if(logDist)
				covDistTens(n1,n2) = computeLogCovDist(S1,S2);
			else
				covDistTens(n1,n2) = computeCovDist(S1,S2);
			end
			if symKernel
				covDistTens(n2,n1) = covDistTens(n1,n2);
			end
		end
	end
end

function covDist = computeCovDist(S1,S2)
%computeCovDist Compute the non-geodesic SPD distance between S1 and S2

covDist = norm(S1-S2,'fro');
end

function covDist = computeLogCovDist(S1,S2)
%computeCovDist Compute the non-geodesic SPD distance between S1 and S2

covDist = norm(logm(S1)-logm(S2),'fro');
end

function M = tens2mat(T, rows, cols)
%tens2mat Matricisation of a tensor T such that dimensions contained in
%rows and cols become the rows, respectively columns of the new matrix M

if nargin <=2
	cols = [];
end

sizeT = size(T);
N = ndims(T);

if isempty(rows)
	rows = 1:N;
	rows(cols) = [];
end
if isempty(cols)
	cols = 1:N;
	cols(rows) = [];
end

if any(rows(:)' ~= 1:length(rows)) || any(cols(:)' ~= length(rows)+(1:length(cols)))
	T = permute(T,[rows(:)' cols(:)']);
end

M = reshape(T,prod(sizeT(rows)), prod(sizeT(cols)));

end