%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Kevin Fries (2016)
%This is an example script for learning the hyperparameters of a Gaussian
%Process Regression model and then making estimates at unobserved locations
%by using the GPML toolbox created by Rasmussen
%and Williams. To use this script, you MUST have the GPML toolbox in your
%MATLAB path or added to your path via your script.

%This script uses Lake Erie Summer Air Temperatures for 2013 as an example.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all, close all

year = 2013;
file = sprintf('./collocated ship and model data/erie%g.mat',year);
load(file);

errorTemp = shipTemp - modelTemp; %create zero mean process to be modeled

a = find(abs(errorTemp)<20); %eliminate extreme outliers
errorTemp = errorTemp(a);
LatTemp = LatTemp(a);
LonTemp = LonTemp(a);
modelTemp = modelTemp(a);
shipTemp = shipTemp(a);
TimeTemp = TimeTemp(a);

mar31 = datenum(year,3,31,23,59,59);
jun30 = datenum(year,6,30,23,59,59);
sep30 = datenum(year,9,30,23,59,59);
dec31 = datenum(year,12,31,23,59,59);

winter = find(TimeTemp < mar31);
if isempty(winter)
    winter = 1;
end
spring = find(TimeTemp(winter(end):end)<jun30);
spring = spring + winter(end)-1;
summer = find(TimeTemp(spring(end):end)<sep30);
summer = summer + spring(end)-1;
fall = find(TimeTemp(summer(end):end)<dec31);
fall = fall + summer(end)-1;

% normalize data to lie between 0 and 1 such that length scales can be
% compared
minLat = min(LatTemp);
maxLat = max(LatTemp);
minLon = min(LonTemp);
maxLon = max(LonTemp);
minTemp = min(modelTemp);
maxTemp = max(modelTemp);

LatTemp = (LatTemp - minLat)./(maxLat - minLat);
LonTemp = (LonTemp - minLon)./(maxLon - minLon);
Temps = (modelTemp - minTemp)./(maxTemp - minTemp);


%% summer

%define training and testing sets
train_summer = randsample(summer,fix(length(summer)/2));
test_summer = setdiff(summer,train_summer);

x = [LatTemp(train_summer),LonTemp(train_summer),Temps(train_summer)]; %input to GP
y = errorTemp(train_summer); %output for GP
x_star = [LatTemp(test_summer),LonTemp(test_summer),Temps(test_summer)]; %testing set

covfunc = @covSEard_nosf;
% we use a modified Squared Exponential kernel that eliminates the
% scale factor. This makes comparison of length scales between
% different random restarts possible. This covariance function is
% provided in this repository

likfunc = @likGauss; % use a Gaussian likelihood function (i.e. assume Gaussian noise on measurements)
    
%initialize hyperparameters
hyps_summer.cov = [rand() rand() rand()];
hyps_summer.lik = rand()+1;

%learn hyperparameters. We only perform 1 random start here, in practice
%it is good to perform multiple restarts and choose the best
%hyperparameters from the set of random initializations
hyps_summer = minimize(hyps_summer, @gp, -100, @infExact, [], covfunc, likfunc, x, y);

%use learned hyperparameters to predict values at "unobserved" locations.
%In reality we have obs here, but we did not use them to learn the model
%and therefore it is equivalent to being unobserved. The outputs here are
%'m' which is the error (i.e. the value to be added to the GLCFS initial
%guess) and 's2' which is the variance in that estimate (i.e. confidence
%level)
[m,s2] = gp(hyps_summer, @infExact, [], covfunc, likfunc, x, y, x_star);


%% generate figure
figure('units','normalized','outerposition',[0 0 1 1])
subplot(2,2,1)
scatter(shipTemp(test_summer),modelTemp(test_summer),'r')
xlabel('Ship Measurement (degrees C)')
ylabel('GLCFS Estimate (degrees C)')
refline(1,0)


subplot(2,2,2)
scatter(shipTemp(test_summer),modelTemp(test_summer) + m,'r')
refline(1,0)
xlabel('Ship Measurement (degrees C)')
ylabel('GP Estimate (degrees C)')


subplot(2,2,3)
hist(shipTemp(test_summer) - modelTemp(test_summer))
mean1 = mean(shipTemp(test_summer) - modelTemp(test_summer));
std1 = std(shipTemp(test_summer) - modelTemp(test_summer));
xaxis = sprintf('Mean = %g, Standard Deviation = %g',mean1,std1);
xlabel(xaxis)

subplot(2,2,4)
hist(shipTemp(test_summer) - modelTemp(test_summer) + m)
mean2 = mean(shipTemp(test_summer) - (modelTemp(test_summer) + m));
std2 = std(shipTemp(test_summer) - (modelTemp(test_summer) + m));
xaxis = sprintf('Mean = %g, Standard Deviation = %g',mean2,std2);
xlabel(xaxis)

titles = sprintf('Erie Summer Temp %g \r\n Hyperparameters: Lat=%g,Lon=%g,Model=%g,Noise=%g'...
    ,year,exp(hyps_summer.cov(1)),exp(hyps_summer.cov(2)),exp(hyps_summer.cov(3)),exp(hyps_summer.lik));
suptitle(titles)