Commit 21c1e90d authored by Jordi's avatar Jordi
Browse files

Scripts for the 1st revision experiments, and some fixes.

parent ab7db222

.gitignore

+1 −0
Original line number Diff line number Diff line 
*.pkl

*.mat

*.json

.vscode

.ipynb_*

__pycache__

matlab/demSPARC_LFM_rev1.m

0 → 100644
+193 −0
Original line number Diff line number Diff line 
% MULTIGP SPARC, two LFs, individual GPs

clc, clear

rng('default'), rng(0)



% DONE: 2 LFs, 2 GPs, one for LAI, other for FAPAR

% TODO:

% - LMC with 1 LF, all variables, for comparing vs convolution GG kernel

% - If possible GG kernel + periodic kernel



% modelType = 'sim';  % First order differential equation kernel (see simMeanCreate.m)

modelType = 'gg';   % Gaussian kernel (for convolution) and Gaussian covariance

% modelType = 'lmc';   % Linear model of corregionalization



nlf = 2;                  % number of latent forces

missing_rate = 0;         % set to 0.9 for the 90% missing rate experiment

gf_experiment = true;     % set to true for the gap-filling experiment

gf_range = [2009,2014];   % setup range (from 2019 to 2010-2014)

var = 'both';             % LAI, FAPAR or both



%% Set needed paths

% addpath(genpath('../../../../Software/SheffieldML/GPmat'))

% addpath('../../../../Software/SheffieldML/netlab')

% addpath('../../../../Software/SheffieldML/multigp/matlab')

addpath(genpath('GPmat'))

addpath('netlab')

addpath('multigp/matlab')



%% Load data

load ../SPARC/timeseries_1pixel_3countries



% Named wavelength, is time actually

% _MOD comes from MODIS, I guess

tt = wavelength_MOD;

% tt = wavelength_VGT;

% Center inputs

% mt = mean(tt); % raises problems with cov. matrices

% mt = 2003; % ad-hoc normalization

mt = 0; % set mean to zero to avoid numerical problems

tt = tt - mt;



% X / Y training sets

idx{1} = IT_LAI_MOD > 0;

idx{2} = ES_LAI_MOD > 0;

idx{3} = IT_FAPAR_MOD > 0;

idx{4} = ES_FAPAR_MOD > 0;

% idx{1} = IT_LAI_VGT > 0;

% idx{2} = ES_LAI_VGT > 0;

% idx{3} = IT_FAPAR_VGT > 0;

% idx{4} = ES_FAPAR_VGT > 0;



is = 0;

if strcmp(var, 'LAI') || strcmp(var, 'both')

    Xtrain{is+1} = tt(idx{1});

    Xtrain{is+2} = tt(idx{2});

    % LAI

    Ytrain{is+1} = IT_LAI_MOD(idx{1});

    % Ytrain{is+1} = IT_LAI_VGT(idx{1});

    my(is+1) = 0;  % my(is+1) = mean(Ytrain{is+1}); Ytrain{is+1} = Ytrain{is+1} - my(is+1);

    Ytrain{is+2} = ES_LAI_MOD(idx{2});

    % Ytrain{is+2} = ES_LAI_VGT(idx{2});

    my(2) = 0;  % my(2) = mean(Ytrain{2}); Ytrain{2} = Ytrain{2} - my(2);

    is = is + 2;

end

if strcmp(var, 'FAPAR') || strcmp(var, 'both')

    Xtrain{is+1} = tt(idx{3});

    Xtrain{is+2} = tt(idx{4});

    % FAPAR

    Ytrain{is+1} = IT_FAPAR_MOD(idx{3});

    % Ytrain{is+1} = IT_FAPAR_VGT(idx{3});

    my(is+1) = 0;  % my(is+1) = mean(Ytrain{1}); Ytrain{is+1} = Ytrain{1} - my(is+1);

    Ytrain{is+2} = ES_FAPAR_MOD(idx{4});

    % Ytrain{is+2} = ES_FAPAR_VGT(idx{4});

    my(is+2) = 0;  % my(is+2) = mean(Ytrain{2}); Ytrain{is+2} = Ytrain{2} - my(is+2);

end



if strcmp(var, 'both'), var = 'LAI_FAPAR'; end



clear *MOD *VGT



%% Apply missing rate

if missing_rate > 0 && missing_rate < 1

    fprintf('Missing rate %0.2f\n', missing_rate)

    Xtest = cell(size(Xtrain));

    Ytest = cell(size(Ytrain));

    for i = 1:length(Xtrain)

        Xtest{i} = Xtrain{i};

        Ytest{i} = Ytrain{i};

        n = fix(length(Xtrain{i}) * (1-missing_rate));

        ridx = randperm(length(Xtrain{i}));

        Xtrain{i} = Xtrain{i}(ridx(1:n));

        Ytrain{i} = Ytrain{i}(ridx(1:n));

    end

end



%% Gap-filling

if gf_experiment

    fprintf('Gap-filling, removed range: %d - %d\n', gf_range(1), gf_range(2))

    Xtest = cell(size(Xtrain));

    Ytest = cell(size(Ytrain));

    for i = 1:length(Xtrain)

        Xtest{i} = Xtrain{i};

        Ytest{i} = Ytrain{i};

        % Left ES_LAI (or ES_FAPAR) untouched

        if i == 2, continue, end

        % Keep only range data

        ridx = ((Xtrain{i}+mt) < gf_range(1)) | ((Xtrain{i}+mt) > gf_range(2));

        Xtrain{i} = Xtrain{i}(ridx);

        Ytrain{i} = Ytrain{i}(ridx);

    end

end



% Show training set sizes

disp(Ytrain)



% Filename for model

mdate = datestr(now, 'yyyymmdd_hhMM');

switch lower(modelType)

    case {'sim', 'gg'}

        fname = sprintf('SPARC_model_rev1_%s_%s_nlf_%d_mr_%0.2f_%s.mat', ...

            var, upper(modelType), nlf, missing_rate, mdate);

    case 'lmc'

        fname = sprintf('SPARC_model_rev1_%s_%s_nlf_%d_rank_%d_mr_%0.2f_%s.mat', ...

            var, upper(modelType), nlf, size(Ytrain,2), missing_rate, mdate);

end

if gf_experiment

    % fname = strrep(fname, '_SPARC_', '_SPARC_gapfilling_');

    srange = sprintf('_model_range_%d_%d_', gf_range(1), gf_range(2));

    fname = strrep(fname, '_model_', srange);

end

% Filename for results

rname = strrep(fname, '_model_', '_results_');



% Show output filenames

disp(fname)

disp(rname)



%% Options for multiGP model

% Common options

options = multigpOptions('ftc');

options.optimiser = 'scg';



% Specific options

options.kernType = modelType;

switch lower(modelType)

    case 'sim'

        options.nlf = nlf;

    case 'gg'

        % Options for GG Full Cov. Latent Force model

        options.nlf = nlf; % size(Ytrain, 2); % <= 3 latent forces in this case

    case 'lmc'

        % Options for linear model of coregionalization (LMC).

        % Here the number of effective latent forces are given by

        % nlf * rankCorregMatrix. As far as I understand, we should set

        % nlf = 1 always and play only with the rank of the coreg. matrix,

        % which defines the number of underlying GP models (i.e., the number of

        % underlying latent forces.

        % Despite of that, a quite interesting result is obtained setting

        % nlf = 2 and rankCorregMatrix = 1. The 1st LF seems related to the low

        % frequency cycle, whereas the 2nd LF is the high freq. terms between

        % pics.

        options.nlf = nlf;

        options.rankCorregMatrix = size(Ytrain, 2);

        options.kern.nout = size(Ytrain, 2);

        options.kern.rankCorregMatrix = options.rankCorregMatrix;

    otherwise

        error(['Unkwnon model type ' modelType])

end



% Build and train model

model = trainLFMGP(Xtrain, Ytrain, options);



% Save results

save(['models_rev1/' fname], 'model')



% Model info

modelDisplay(model)



% Make predictions and save results

[mu, varsigma] = multigpPosteriorMeanVar(model, tt);



% Add means to all variables

for i = 1:length(Xtrain), Xtrain{i} = Xtrain{i} + mt; end

for i = 1:length(Ytrain), Ytrain{i} = Ytrain{i} + my(i); end

for i = nlf+1:length(mu), mu{i} = mu{i} + my(i-nlf); end



save(['results_rev1/' rname], 'Xtrain', 'Ytrain', 'idx', 'mt', 'my', 'tt', 'mu', 'varsigma', 'var')



if exist('Xtest', 'var')

    save(['results_rev1/' rname], 'Xtest', 'Ytest', '-append')

end



%% Show results

% ggResults(model, Xtrain, Ytrain, Xtest, Ytest);

matlab/figures_SPARC_rev1.m

0 → 100644
+471 −0
Original line number Diff line number Diff line 
clear, clc



nicefig

% figures = 'lai-fapar-2'; % 'full', 'gf', 'gf_gg', 'mr', 'mr_gg', 'lai-fapar', 'lai-fapar-2', 'numerical';

% figures = 'smos';



% Rev. 1, experiment 1

% figures = 'gf_gg';

% var = 'LAI';  % LAI or FAPAR



% Rev. 2, experiment 2

figures = 'lmc';

var = 'LAI_FAPAR';



% For both experiments

range = 2014;



switch figures

    case {'full', 'lai-fapar', 'lai-fapar-2'}

        % Using _MOD data

        load results/SPARC_results_GG_nlf_1_mr_0.00_20200622_2006.mat

        % Using _VGT data

        % load results/SPARC_results_GG_nlf_1_mr_0.00_20200625_1224.mat

    case 'gf'

        % SIM results

        % load ../results/predictions_gap_filling.mat

    case 'gf_gg'

        % New GG results

        % load results/SPARC_results_range_2009_2014_GG_nlf_1_mr_0.00_20200626_0920.mat

        fname = dir(sprintf('results_rev1/SPARC_results_range_2009_%d_rev1_%s_GG_*.mat', range, var));

        fname = ['results_rev1/' fname(1).name];

        load(fname)

    case 'lmc'

        fname = dir(sprintf('results_rev1/SPARC_results_range_2009_%d_rev1_%s_LMC_*.mat', range, var));

        fname = ['results_rev1/' fname(1).name];

        load(fname)        

    case 'mr'

        % load ../results/predictions_missing_rate_0.90.mat

    case 'mr_gg'

        load results/SPARC_results_GG_nlf_1_mr_0.90_20200626_0928.mat

    case 'numerical'

        % Look at the end of this script

        load results/SPARC_results_range_2009_2010_GG_nlf_1_mr_0.00_20200626_1544.mat

    case 'smos'

        % load results/SMOS_results_SIM_site_DAHRA_nlf_1_sensors_3_20200720_1703.mat

        % load results/SMOS_results_SIM_site_DAHRA_nlf_3_sensors_3_20200720_1729.mat

        load '../../../2020 SMOS/AutocorKerData/gapFilling_REMEDHUS.mat' ...

            utc_ref SM_both SM_ascat SM_amsr2

        % input/output

        X = utc_ref;

        Y = [SM_both SM_ascat SM_amsr2];

        my = nanmean(Y);

        Y = Y - repmat(my, size(Y,1), 1);

        load results/SMOS_results_SIM_site_REMEDHUS_nlf_3_sensors_3_20200720_1942.mat

    otherwise

        error(['Unknown figures mode ' figures])

end



nlf = length(mu) - length(Ytrain);

xl = [2003, 2014];



countries = ['Italy', 'Spain'];

sets = cell(1,2);

for c = 1:length(countries)

    sets{c} = sprintf('%s (%s)', var, countries(c));

end

pv = @(t,y,v) fill([t ; flipud(t)], [y+2*sqrt(v) ; flipud(y-2*sqrt(v))], ...

        0.4*ones(1,3), 'FaceAlpha', 0.2, 'EdgeColor', 'none');



%% Make figures

switch figures

    case 'full'

        % LF

        id = 1;

        figure(1), clf, hold on

        pv(tt, mu{id}, varsigma{id})

        plot(tt, mu{id}, 'k')

        xlabel('Year'), ylabel('Latent Force')

        xlim(xl), ylim([-1,3])

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigRes1LF1_gg.eps



        % LAI Spain

        id = 2;

        figure(2), clf, hold on

        pv(tt, mu{id+nlf}, varsigma{id+nlf})

        plot(Xtrain{id}, Ytrain{id}, '.', 'MarkerSize', 16)

        plot(tt, mu{id+nlf}, 'b')

        xlabel('Year'), ylabel('LAI [m^2/M^2] (Spain)')

        xlim(xl)

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigRes1SpLAI_gg.eps



        % fAPAR Spain

        id = 4;

        figure(3), clf, hold on

        pv(tt, mu{id+nlf}, varsigma{id+nlf})

        plot(Xtrain{id}, Ytrain{id}, '.', 'MarkerSize', 16)

        plot(tt, mu{id+nlf}, 'b')

        xlabel('Year'), ylabel('fAPAR (Spain)')

        xlim(xl)

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigRes1SpfAPAR_gg.eps



        % Numerical results

        for i = 1:length(Ytrain)

            mse = mean((Ytrain{i} - mu{i+nlf}(idx{i})).^2);

            nmse = mse / mean(Ytrain{i}.^2) * 100;

            fprintf('%-14s MSE: %0.4f, NMSE: %0.2f\n', sets{i}, mse, nmse)

        end



    case 'gf'

        % LAI Italy

        id = 1;

        figure(1), clf, hold on

        pv(tt, Yp(id,:)', Vp(id,:)')

        plot(T{id}, Y{id}, '.', 'MarkerSize', 16)

        plot(Ttest{id}, Ytest{id}, 'k--')

        plot(tt, Yp(id,:), 'b')

        xlabel('Year'), ylabel('LAI [m^2/M^2] (Italy)')

        xlim(xl)

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigRes3ItLAI.eps



        % fAPAR Italy

        id = 3;

        figure(2), clf, hold on

        pv(tt, Yp(id,:)', Vp(id,:)')

        plot(T{id}, Y{id}, '.', 'MarkerSize', 16)

        plot(Ttest{id}, Ytest{id}, 'k--')

        plot(tt, Yp(id,:), 'b')

        xlabel('Year'), ylabel('fAPAR (Italy)')

        xlim(xl), ylim([-0.2, 1.2])

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigRes3ItfAPAR.eps



        % fAPAR Spain

        id = 4;

        figure(3), clf, hold on

        pv(tt, Yp(id,:)', Vp(id,:)')

        plot(T{id}, Y{id}, '.', 'MarkerSize', 16)

        plot(Ttest{id}, Ytest{id}, 'k--')

        plot(tt, Yp(id,:), 'b')

        xlabel('Year'), ylabel('fAPAR (Spain)')

        xlim(xl), ylim([0, 1.2])

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigRes3SpfAPAR.eps



    case 'gf_gg'

        switch var

            case 'LAI'

                % LAI Italy

                id = 1;

                figure(1), clf, hold on

                pv(tt, mu{id+nlf}, varsigma{id+nlf})

                plot(Xtrain{id}, Ytrain{id}, '.', 'MarkerSize', 16)

                plot(Xtest{id}, Ytest{id}, 'k--')

                plot(tt, mu{id+nlf}, 'b')

                xlabel('Year'), ylabel('LAI [m^2/M^2] (Italy)')

                xlim(xl)

                grid on

                hold off

                print -depsc2 ./figures_rev1/FigRes3ItLAI_gg.eps



                % LAI Spain

                id = 2;

                figure(2), clf, hold on

                pv(tt, mu{id+nlf}, varsigma{id+nlf})

                plot(Xtrain{id}, Ytrain{id}, '.', 'MarkerSize', 16)

                plot(Xtest{id}, Ytest{id}, 'k--')

                plot(tt, mu{id+nlf}, 'b')

                xlabel('Year'), ylabel('LAI [m^2/M^2] (Spain)')

                xlim(xl)

                grid on

                hold off

                print -depsc2 ./figures_rev1/FigRes3SpLAI_gg.eps



            case 'FAPAR'

                % fAPAR Italy

                id = 1;

                figure(3), clf, hold on

                pv(tt, mu{id+nlf}, varsigma{id+nlf})

                plot(Xtrain{id}, Ytrain{id}, '.', 'MarkerSize', 16)

                plot(Xtest{id}, Ytest{id}, 'k--')

                plot(tt, mu{id+nlf}, 'b')

                xlabel('Year'), ylabel('fAPAR (Italy)')

                xlim(xl), ylim([-0.2, 1.2])

                grid on

                hold off

                print -depsc2 ./figures_rev1/FigRes3ItfAPAR_gg.eps



                % fAPAR Spain

                id = 2;

                figure(4), clf, hold on

                pv(tt, mu{id+nlf}, varsigma{id+nlf})

                plot(Xtrain{id}, Ytrain{id}, '.', 'MarkerSize', 16)

                plot(Xtest{id}, Ytest{id}, 'k--')

                plot(tt, mu{id+nlf}, 'b')

                xlabel('Year'), ylabel('fAPAR (Spain)')

                xlim(xl), ylim([-0.2, 1.2])

                grid on

                hold off

                print -depsc2 ./figures_rev1/FigRes3SpfAPAR_gg.eps



            otherwise

                error(['Unkwnown var ' var])

        end



    case 'lmc'

        % LAI Italy

        id = 1;

        figure(1), clf, hold on

        pv(tt, mu{id+nlf}, varsigma{id+nlf})

        plot(Xtrain{id}, Ytrain{id}, '.', 'MarkerSize', 16)

        plot(Xtest{id}, Ytest{id}, 'k--')

        plot(tt, mu{id+nlf}, 'b')

        xlabel('Year'), ylabel('LAI [m^2/M^2] (Italy)')

        xlim(xl)

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigResRev1ItLAI_lmc.eps



        % LAI Spain

        id = 2;

        figure(2), clf, hold on

        pv(tt, mu{id+nlf}, varsigma{id+nlf})

        plot(Xtrain{id}, Ytrain{id}, '.', 'MarkerSize', 16)

        plot(Xtest{id}, Ytest{id}, 'k--')

        plot(tt, mu{id+nlf}, 'b')

        xlabel('Year'), ylabel('LAI [m^2/M^2] (Spain)')

        xlim(xl)

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigResRev1SpLAI_lmc.eps



        % fAPAR Italy

        id = 3;

        figure(3), clf, hold on

        pv(tt, mu{id+nlf}, varsigma{id+nlf})

        plot(Xtrain{id}, Ytrain{id}, '.', 'MarkerSize', 16)

        plot(Xtest{id}, Ytest{id}, 'k--')

        plot(tt, mu{id+nlf}, 'b')

        xlabel('Year'), ylabel('fAPAR (Italy)')

        xlim(xl), ylim([-0.2, 1.2])

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigResRev1ItfAPAR_lmc.eps



        % fAPAR Spain

        id = 4;

        figure(4), clf, hold on

        pv(tt, mu{id+nlf}, varsigma{id+nlf})

        plot(Xtrain{id}, Ytrain{id}, '.', 'MarkerSize', 16)

        plot(Xtest{id}, Ytest{id}, 'k--')

        plot(tt, mu{id+nlf}, 'b')

        xlabel('Year'), ylabel('fAPAR (Spain)')

        xlim(xl), ylim([-0.2, 1.2])

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigResRev1SpfAPAR_lmc.eps



    case 'mr'

        % LF

        figure(1), clf, hold on

        pv(tt, lf, vf)

        plot(tt, lf, 'k')

        xlabel('Year'), ylabel('Latent Force')

        xlim(xl), ylim([-1,3])

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigRes4LF1.eps



        % fAPAR Italy

        id = 3;

        figure(2), clf, hold on

        pv(tt, Yp(id,:)', Vp(id,:)')

        plot(T{id}, Y{id}, '.', 'MarkerSize', 16)

        plot(Ttest{id}, Ytest{id}, 'k--')

        plot(tt, Yp(id,:), 'b')

        xlabel('Year'), ylabel('fAPAR (Italy)')

        xlim(xl), ylim([-0.4, 1.4])

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigRes4ItfAPAR.eps



        % fAPAR Spain

        id = 4;

        figure(3), clf, hold on

        pv(tt, Yp(id,:)', Vp(id,:)')

        plot(T{id}, Y{id}, '.', 'MarkerSize', 16)

        plot(Ttest{id}, Ytest{id}, 'k--')

        plot(tt, Yp(id,:), 'b')

        xlabel('Year'), ylabel('fAPAR (Spain)')

        xlim(xl), ylim([-0.4, 1.4])

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigRes4SpfAPAR.eps



    case 'mr_gg'

        % LF

        figure(1), clf, hold on

        pv(tt, mu{nlf}, varsigma{nlf})

        plot(tt, mu{nlf}, 'k')

        xlabel('Year'), ylabel('Latent Force')

        xlim(xl), ylim([-1.5,3.5])

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigRes4LF1_gg.eps



        % fAPAR Italy

        id = 3;

        figure(2), clf, hold on

        pv(tt, mu{id+nlf}, varsigma{id+nlf})

        plot(Xtrain{id}, Ytrain{id}, '.', 'MarkerSize', 16)

        plot(Xtest{id}, Ytest{id}, 'k--')

        plot(tt, mu{id+nlf}, 'b')

        xlabel('Year'), ylabel('fAPAR (Italy)')

        xlim(xl), ylim([-0.4, 1.4])

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigRes4ItfAPAR_gg.eps



        % fAPAR Spain

        id = 4;

        figure(3), clf, hold on

        pv(tt, mu{id+nlf}, varsigma{id+nlf})

        plot(Xtrain{id}, Ytrain{id}, '.', 'MarkerSize', 16)

        plot(Xtest{id}, Ytest{id}, 'k--')

        plot(tt, mu{id+nlf}, 'b')

        xlabel('Year'), ylabel('fAPAR (Spain)')

        xlim(xl), ylim([-0.4, 1.4])

        grid on

        hold off

        print -depsc2 ./figures_rev1/FigRes4SpfAPAR_gg.eps



    case 'lai-fapar'

        id = [1,3];  % id = [2,4] for Spain

        id = id + nlf;

        % data_lai = Ytrain{id(1)-nlf}; data_fapar = Ytrain{id(2)-nlf};

        lai = mu{id(1)}; fapar = mu{id(2)};

        vp_lai = varsigma{id(1)}; vp_fapar = varsigma{id(2)};

        % Sort variables

        [lai, ids] = sort(lai);

        fapar = fapar(ids);

        % Linear model

        alpha = (log(1-fapar)'*lai) / (lai'*lai);

        model_lai = linspace(lai(1), lai(end));

        model_fapar = 1 - exp(alpha * model_lai);

        figure(4), clf, hold on

        % Draw 2 stds as ellipses

        t = 0:pi/12:2*pi;

        sint = sin(t); cost = cos(t);

        for i = 1:length(lai)

            xc = lai(i) + 2*sqrt(vp_lai(i)) * sint;

            yc = fapar(i) + 2*sqrt(vp_fapar(i)) * cost;

            fill(xc, yc, [0.4, 0.4, 0.4], 'LineStyle', 'none', 'FaceAlpha', 0.01);

        end

        % Draw data, prediction and linear model

        % plot(data_lai, data_fapar, 'x', 'MarkerSize', 16)

        plot(lai, fapar, 'k.', 'MarkerSize', 16)

        plot(model_lai, model_fapar, 'r', 'LineWidth', 2)

        xlim([-2,6.5]), ylim([-0.2,1.2])

        xlabel('LAI [m^2/m^2] (Spain)'), ylabel('fAPAR (Spain)')

        grid on



        print -depsc2 ./figures_rev1/FigRes2_gg.eps

        print -dpng ./figures_rev1/FigRes2_gg.png



    case 'lai-fapar-2'

        id = [1,3];  % id = [2,4] for Spain

        id = id + nlf;

        % data_lai = Ytrain{id(1)-nlf}; data_fapar = Ytrain{id(2)-nlf};

        lai = mu{id(1)}; fapar = mu{id(2)};

        vp_lai = varsigma{id(1)}; vp_fapar = varsigma{id(2)};

        % Sort variables

        [lai, ids] = sort(lai);

        fapar = fapar(ids);



        % Linear model 1

        range = lai >=0 & lai <= 2;

        model1 = regress(fapar(range), [lai(range), ones(size(lai(range)))]);

        model1_x = linspace(-0.2, 2.0);

        model1_y = model1_x * model1(1) + model1(2);



        % Linear model 2

        range = lai >=2 & lai <= 5;

        model2 = regress(fapar(range), [lai(range), ones(size(lai(range)))]);

        model2_x = linspace(1.2, 5.5);

        model2_y = model2_x * model2(1) + model2(2);



        % Linear model 3

        %range = lai >=4 & lai <= 6;

        %model3 = regress(fapar(range), [lai(range), ones(size(lai(range)))]);

        %model3_x = linspace(3.5, 5.5);

        %model3_y = model3_x * model3(1) + model3(2);



        figure(4), clf, hold on

        % Draw 2 stds as ellipses

        t = 0:pi/12:2*pi;

        sint = sin(t); cost = cos(t);

        for i = 1:length(lai)

            xc = lai(i) + 2*sqrt(vp_lai(i)) * sint;

            yc = fapar(i) + 2*sqrt(vp_fapar(i)) * cost;

            %fill(xc, yc, [0.6, 0.6, 0.6], 'LineStyle', 'none', 'FaceAlpha', 0.01);

            fill(xc, yc, [1, 1, 1], 'LineStyle', '-', 'FaceAlpha', 0.01, ...

                'EdgeColor', [0, 0, 0.6], 'EdgeAlpha', 0.05);

        end

        % Draw data, prediction and linear model

        % plot(data_lai, data_fapar, 'x', 'MarkerSize', 16)

        plot(lai, fapar, '.', 'Color', [0, 0, 0.6], 'MarkerSize', 16)

        plot(model1_x, model1_y, 'r--', 'LineWidth', 2)

        plot(model2_x, model2_y, 'r--', 'LineWidth', 2)

        %plot(model3_x, model3_y, 'r-', 'LineWidth', 2)

        xlim([-1,6]), ylim([-0.1,1.2])

        xlabel('LAI [m^2/m^2] (Spain)'), ylabel('fAPAR (Spain)')

        grid on



        fprintf('alpha1 %f\n', model1(1))

        fprintf('alpha2 %f\n', model2(1))



        print -depsc2 ./figures_rev1/FigRes2_gg2.eps

        print -dpng ./figures_rev1/FigRes2_gg2.png



    case 'numerical'

        files = dir('../matlab/SPARC_results_range_2009_*_GG_*.mat');

        for f = files'

            % disp(f.name)

            load(sprintf('../matlab/%s', f.name))

            ii = strfind(f.name, '_2009');

            year = str2double(f.name(ii+6:ii+9));

            fprintf('2009--%d', year-1)

            for i = 1:length(Ytest)

                if i == 2, continue, end  % Skip Spain

                mse = mean((Ytest{i} - mu{i+nlf}(idx{i})).^2);

                nmse = mse / mean(Ytest{i}.^2) * 100;

                fprintf(' & %0.4f (%0.2f \\%%)', mse, nmse)

            end

            fprintf(' \\\\\n');

        end



    case 'smos'

        figure(1), clf

        cm = get(gca, 'colororder');

        for i = 1:3

            subplot(3,1,i), hold on

            plot(Xtrain{i}, Ytrain{i}+my(i), '.', 'MarkerSize', 12, 'Color', cm(2,:))

            plot(Xtest{i}, mu{nlf+i}+my(i), 'k-')

            pv(Xtest{i}, mu{i+nlf}+my(i), varsigma{i+nlf})

            axis tight, grid on

            datetick keeplimits, hold off

        end

        figure(2), clf

        for i = 1:3

            subplot(3,1,i), hold on

            plot(Xtrain{i}, Ytrain{i}+my(i), '.', 'MarkerSize', 12, 'Color', cm(2,:))

            plot(Xtest{i}, mu{nlf+i}+my(i), 'k-')

            pv(Xtest{i}, mu{i+nlf}+my(i), varsigma{i+nlf})

            axis tight

            xlim([datenum('2012-05-01'), datenum('2012-10-31')]), grid on

            datetick keeplimits, hold off

        end

        legend('Data', 'Prediction')

end



shg

matlab/results_SPARC.m

+38 −28
Original line number Diff line number Diff line 
clear, clc



nlf = 1;

vars = {'LAI (IT)', 'LAI (ES)', 'fAPAR (IT)', 'fAPAR (ES)'};



% load model_SPARC_results_range_2009_2010_SIM_nlf_1_mr_0.00_20200622_1433.mat
@@ -12,16 +11,27 @@ vars = {'LAI (IT)', 'LAI (ES)', 'fAPAR (IT)', 'fAPAR (ES)'}; 
% end



disp('LaTex table')

files = dir('model_SPARC_results_range_2009_*_SIM*');

% files = dir('model_SPARC_results_range_2009_*_GG*');

% files = dir('model_SPARC_results_range_2009_*_SIM*.mat');

% files = dir('model_SPARC_results_range_2009_*_GG*.mat');

% files = dir('results/SPARC_results_range_2009_*_GG*.mat');



% LAI IT

% files = dir('results_rev1/SPARC_results_range_2009_*_rev1_LAI_GG*.mat');

% FAPAR IT

% files = dir('results_rev1/SPARC_results_range_2009_*_rev1_FAPAR_GG*.mat');

% LAI & FAPAR using LMC and 2 LFs

% files = dir('results_rev1/SPARC_results_range_2009_*_rev1_LAI_FAPAR_LMC*.mat');

% LAI & FAPAR using GG and 2 LFs

files = dir('results_rev1/SPARC_results_range_2009_*_rev1_LAI_FAPAR_GG*.mat');

for f = files'

    y = str2double(f.name(32:35));

    y = str2double(f.name(26:29));

    fprintf('2009-%d', y-1)

    load(f.name)

    load(['results_rev1/' f.name])

    nlf = length(mu) - length(Xtrain);

    for i = 1:length(Ytest)

        mse = mean((Ytest{i} - mu{i+nlf}(idx{i})).^2);

        nmse = mse / mean(Ytest{i}.^2) * 100;

        if i == 2, continue, end  % Skip LAI (ES)

        if i == 2, continue, end  % Skip LAI (ES) or FAPAR (ES)

        fprintf(' & %0.4f (%0.2f \\%%)', mse, nmse)

    end

    fprintf(' \\\\\n')

matlab/trainLFMGP.m

+4 −2

File changed.

Preview size limit exceeded, changes collapsed.