%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Kevin Fries (2016)
%This is an example script for learning the hyperparameters of a Gaussian
%Process Regression model 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 Air Temperatures as an example of how to
%generate the FINAL hyperparameters (i.e. after verifying the model works).
%To see how to then use the hyperparameters to generate updated estimates
%at unobserved locations, see the second part of generate_estimates_example.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all, close all
%initialize variables
bestparams_winter=[];
bestliks_winter = [];
params_winter = [];
liks_winter = [];
minlike_winter = [];

bestparams_spring=[];
bestliks_spring = [];
params_spring = [];
liks_spring = [];
minlike_spring = [];

bestparams_summer=[];
bestliks_summer = [];
params_summer = [];
liks_summer = [];
minlike_summer = [];

bestparams_fall=[];
bestliks_fall = [];
params_fall = [];
liks_fall = [];
minlike_fall = [];

for year = 2006:2014
    % load in file
    file = sprintf('./collocated ship and model data/erie%g.mat',year);
    load(file);

    errorTemp = shipTemp - modelTemp; %create zero-mean process you seek to estimate
    a = find(abs(errorTemp)<20); %eliminate extreme outliers
    LatTemp = LatTemp(a);
    LonTemp = LonTemp(a);
    modelTemp = modelTemp(a);
    shipTemp = shipTemp(a);
    TimeTemp = TimeTemp(a);

    %break data into seasons
    jan1 = datenum(year,1,1);
    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);

    a = find(TimeTemp>=jan1);
    b = find(TimeTemp<=mar31);
    winter = a(find(ismember(a,b)));

    a = find(TimeTemp>mar31);
    b = find(TimeTemp<=jun30);
    spring = a(find(ismember(a,b)));

    a = find(TimeTemp>jun30);
    b = find(TimeTemp<=sep30);
    summer = a(find(ismember(a,b)));

    a = find(TimeTemp>sep30);
    b = find(TimeTemp<=dec31);
    fall = a(find(ismember(a,b)));
    
    errorTemp = shipTemp - modelTemp;

    minLat = min(LatTemp);
    maxLat = max(LatTemp);
    minLon = min(LonTemp);
    maxLon = max(LonTemp);
    minTemp = min(modelTemp);
    maxTemp = max(modelTemp);

    %normalize data between 0 and 1 to enable comparison of length scales
    LatTemp = (LatTemp - minLat)./(maxLat - minLat);
    LonTemp = (LonTemp - minLon)./(maxLon - minLon);
    Temps = (modelTemp - minTemp)./(maxTemp - minTemp);
    
    %initialize hyperparameters
    hyp_winter.lik = -10000;
    hyp_winter.cov = [0 0 0];
    hyp_spring.lik = -10000;
    hyp_spring.cov = [0 0 0];
    hyp_summer.lik = -10000;
    hyp_summer.cov = [0 0 0];
    hyp_fall.lik = -10000;
    hyp_fall.cov = [0 0 0];
    
    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)
    
    %perform 5 random restarts and choose hyperparameters with largest
    %measurement noise (hyp.lik) to yield the least overfitted model. This is
    %the method chosen for this study, other options include choosing the
    %hyperparameters with the largest likelihood (smallest neg log
    %likelihood) by executing the command:
    % nlml = gp(hyp, @infExact, [], covfunc, likfunc, x, y)
    for i = 1:5
        %define inputs
        x_winter = [LatTemp(winter),LonTemp(winter),Temps(winter)];
        %define outputs
        y_winter = errorTemp(winter);

        %randomly initialize hyperparameters
        fprintf('Year=%g \r\n Iteration=%g \r\n winter \r\n',year,i)
        hyps_winter.cov = [rand() rand() rand()];
        hyps_winter.lik = rand()+1;
        hyps_winter = minimize(hyps_winter, @gp, -100, @infExact, [], covfunc, likfunc, x_winter, y_winter);
        
        %check if noise is larger than the current largest noise value
        if hyps_winter.lik>hyp_winter.lik
           hyp_winter = hyps_winter;
        end
        
        %store hyperparameters
        params_winter = [params_winter;year,hyps_winter.cov];
        liks_winter = [liks_winter;year,hyps_winter.lik];
        
        x_spring = [LatTemp(spring),LonTemp(spring),Temps(spring)];%,TimeTemp(train_range)];
        y_spring = errorTemp(spring);

        fprintf('Year=%g \r\n Iteration=%g \r\n spring\r\n ',year,i)
        hyps_spring.cov = [rand() rand() rand()];% rand()];
        hyps_spring.lik = rand()+1;
        hyps_spring = minimize(hyps_spring, @gp, -100, @infExact, [], covfunc, likfunc, x_spring, y_spring);
        
        if hyps_spring.lik>hyp_spring.lik
           hyp_spring = hyps_spring;
        end
        params_spring = [params_spring;year,hyps_spring.cov];
        liks_spring = [liks_spring;year,hyps_spring.lik];
        
        x_summer = [LatTemp(summer),LonTemp(summer),Temps(summer)];%,TimeTemp(train_range)];
        y_summer = errorTemp(summer);

        fprintf('Year=%g \r\n Iteration=%g \r\n  summer\r\n',year,i)
        hyps_summer.cov = [rand() rand() rand()];% rand()];
        hyps_summer.lik = rand()+1;
        hyps_summer = minimize(hyps_summer, @gp, -100, @infExact, [], covfunc, likfunc, x_summer, y_summer);
        
        if hyps_summer.lik>hyp_summer.lik
           hyp_summer = hyps_summer;
        end
        params_summer = [params_summer;year,hyps_summer.cov];
        liks_summer = [liks_summer;year,hyps_summer.lik];
        
        x_fall = [LatTemp(fall),LonTemp(fall),Temps(fall)];%,TimeTemp(train_range)];
        y_fall = errorTemp(fall);

        fprintf('Year=%g \r\n Iteration=%g \r\n  fall\r\n',year,i)
        hyps_fall.cov = [rand() rand() rand()];% rand()];
        hyps_fall.lik = rand()+1;
        hyps_fall = minimize(hyps_fall, @gp, -100, @infExact, [], covfunc, likfunc, x_fall, y_fall);
        
        if hyps_fall.lik>hyp_fall.lik
           hyp_fall = hyps_fall;
        end
        params_fall = [params_fall;year,hyps_fall.cov];
        liks_fall = [liks_fall;year,hyps_fall.lik];
    end


    bestparams_winter = [bestparams_winter;year,hyp_winter.cov];
    bestliks_winter = [bestliks_winter;year,hyp_winter.lik];
    bestparams_spring = [bestparams_spring;year,hyp_spring.cov];
    bestliks_spring = [bestliks_spring;year,hyp_spring.lik];
    bestparams_summer = [bestparams_summer;year,hyp_summer.cov];
    bestliks_summer = [bestliks_summer;year,hyp_summer.lik];
    bestparams_fall = [bestparams_fall;year,hyp_fall.cov];
    bestliks_fall = [bestliks_fall;year,hyp_fall.lik];
end