clear all
close all
clc

%load data from .csv
dataset=readtable('dataset.csv');

%%choose one of the following options
%compare 2 sites
% site1='GN';
% site2='GL';
% rowsok1=find(strcmp(site1,dataset.site));
% rowsok2=find(strcmp(site2,dataset.site));

% compare 1 site against the others
site1='AL';
rowsok1=find(strcmp(site1,dataset.site));
rowsok2=find(~strcmp(site1,dataset.site));

%compare 1 site against the entire dataset (including tested site)
% site1='AL';
% rowsok1=find(strcmp(site1,dataset.site));
% rowsok2=1:1:size(dataset,1);

%compare 2 years
% yr1=2019;
% yr2=2020;
% rowsok1=table2array(dataset(:,4))==yr1;
% rowsok2=table2array(dataset(:,4))==yr2;

%compare 1 year against the others
% yr1=2021;
% rowsok1=table2array(dataset(:,4))==yr1;
% rowsok2=table2array(dataset(:,4))~=yr1;

%compare 1 year against the entire dataset (including tested year)
% yr1=2021;
% rowsok1=table2array(dataset(:,4))==yr1;
% rowsok2=1:1:size(dataset,1);

data1=table2array(dataset(rowsok1,2:end));
data2=table2array(dataset(rowsok2,2:end));

er1=data1(:,13);
gpp1=data1(:,14);
ta1=data1(:,7);
rs1=data1(:,6);
pr1=data1(:,9);
vwc1=data1(:,11);
doy1=data1(:,4);

er2=data2(:,13);
gpp2=data2(:,14);
ta2=data2(:,7);
rs2=data2(:,6);
pr2=data2(:,9);
vwc2=data2(:,11);
doy2=data2(:,4);

%output file
fileID = fopen('comparisons.txt','w+');

%% ER

fprintf(fileID,'ER model \n');

modeler='y~exp(b1*x1)*(b2+b3*x2+b4*x3+b5*x4)';
good1=(~isnan(er1) & ~isnan(ta1) & ~isnan(vwc1) & ~isnan(pr1) & ~isnan(doy1));
good2=(~isnan(er2) & ~isnan(ta2) & ~isnan(vwc2) & ~isnan(pr2) & ~isnan(doy2));

%starting values for nonlinear fit, to be estimated with 'BuildERModel.m'
a01=[0.03,8.6,0.06,-0.002,-0.02];
a02=[0.06,25.0,0.04,-0.03,-0.01];

z1=[ta1(good1),vwc1(good1),pr1(good1),doy1(good1)];
z2=[ta2(good2),vwc2(good2),pr2(good2),doy2(good2)];

% f1=fitnlm(z1,er1(good1),modeler,a01);
% f2=fitnlm(z2,er2(good2),modeler,a02);

%compute P-values using the shuffling method: function in 'shuffER.m'
[db0,da0,da1,da2,da3,pb0,pa0,pa1,pa2,pa3]=shuffER(z1,z2,er1(good1),er2(good2),a01,a02,100);

%print results on output file
fprintf(fileID,'probabilities: b_0 a_0 a_1 a_2 a_3 \n');
fprintf(fileID,[repmat('%g ',1,5),'\n'],pb0,pa0,pa1,pa2,pa3);


%% GPP

clear good1 z1 good2 z2 a01 a02

fprintf(fileID,' \n');
fprintf(fileID,'GPP model \n');

modelgpp='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2+b5*x3)';
good1=(~isnan(gpp1) & ~isnan(rs1) & ~isnan(vwc1) & ~isnan(doy1));
good2=(~isnan(gpp2) & ~isnan(rs2) & ~isnan(vwc2) & ~isnan(doy2));

%starting values for nonlinear fit, to be estimated with 'BuildGPPModel.m'
a01=[-4.8,-0.05,5.8,0.04,-0.02]; 
a02=[-10.4,-0.02,3.1,0.02,-0.01];
% a01=a02;

z1=[rs1(good1),vwc1(good1),doy1(good1)];
z2=[rs2(good2),vwc2(good2),doy2(good2)];

% f1=fitnlm(z1,gpp1(good1),modelgpp,a01);
% f2=fitnlm(z2,gpp2(good2),modelgpp,a02);

%compute P-values using the shuffling method: function in 'shuffGPP.m'
[dF0,dalpha,dA0,dA1,dA2,pF0,palpha,pA0,pA1,pA2]=shuffGPP(z1,z2,gpp1(good1),gpp2(good2),a01,a02,100);

%print results on output file
fprintf(fileID,'probabilities :  F_0 alpha_0 A_0 A_1 A_2  \n');
fprintf(fileID,[repmat('%g ',1,5),'\n'],pF0,palpha,pA0,pA1,pA2);

fclose(fileID);