clear all
close all
clc

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

rowok=1:1:size(dataset,1);

%to test only one plot uncomment the following lines
% site='GL';
% rowok=find(strcmp(site,dataset.site));

%to test only one year uncomment the following lines
% yr=2021;
% rowok=table2array(dataset(:,4))==yr; %years

data=table2array(dataset(rowok,2:end));
gpp=data(:,14);
rs=data(:,6);

%column numbers of candidate additional variables in 'data'
addedvar=[8,11];

%additional variables provided by the user: uncomment one of the follwoing 
% myi=[]; % empty array to build the model from scratch 
myi=[11,4]; % filled with column numbers to fit the model using preselected vatiables  


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

%% first step: GPP vs rs

fprintf(fileID,'first step: \n');
fprintf(fileID,'F_0 alpha_0 pvals varres varexpl AIC h \n');

% exclude rows with NaN values
good=(~isnan(gpp) & ~isnan(rs));
a1g=[-5;-0.02]; %guess

% Michaelis-Menten nonlinear fit 
modelfun='y~(b1*b2*x1/(b1+b2*x1))';
fit1=fitnlm(rs(good),gpp(good),modelfun,a1g);
a1=fit1.Coefficients.Estimate;
pval=fit1.Coefficients.pValue;
res=fit1.Residuals.Raw;
varexpl=var(fit1.Fitted)/var(gpp(good));
AIC1=length(gpp(good))*log(sum(res.^2))+2*length(a1);
h=lillietest(res);

% write fitted values and staistics on output file
fprintf(fileID,[repmat('%g ',1,7),'\n'],a1,pval,varexpl,AIC1,h); 

%% second step: GPP vs additional variable #1

fprintf(fileID,'\n');
fprintf(fileID,'second step \n');
fprintf(fileID,'variable F_0 alpha_0 a_0 a_1 pvals varexpl AIC h \n');

n=1; % flag for variable selection

% loop over candidate additional variables
if(isempty(myi))
for i=addedvar
    
    clear good1 X a2 f2 fit2
    
    good1=~isnan(data(:,i)) & ~isnan(gpp) & ~isnan(rs); 
    
    % compute the starting values for nonlinear fit
    gpp2=gpp(good1)./(a1(1)*a1(2).*rs(good1)./(a1(1)+a1(2).*rs(good1)));
    f2=regress(gpp2,[ones(size(data(good1,i))) data(good1,i)]);
    a2g=[a1; f2(1,1);f2(2,1)];
    
    % nonlinear fit
    X=[rs(good1),data(good1,i)];
    modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2)';
    fit2=fitnlm(X,gpp(good1),modelfun,a2g); 
    a2=fit2.Coefficients.Estimate;
    pval=fit2.Coefficients.pValue;
    res=fit2.Residuals.Raw;
    varexpl=var(fit2.Fitted)/var(gpp(good1));
    AIC=length(gpp(good1))*log(sum(res.^2))+2*length(a2);
    h=lillietest(res);
    
    fprintf(fileID,['%s ',repmat('%g ',1,11),'\n'],cell2mat(dataset(:,i+1).Properties.VariableNames),a2,pval,varexpl,AIC,h);
    
    if(n==1 && AIC<AIC1) %first variable optimizing the AIC compared to step 1
        best=AIC;
        i1=i;
        n=0;
    elseif(n==0 && AIC<best) %variables optimizing the AIC compared to other varibles in the loop
        best=AIC;
        i1=i;
    end
end
end

fprintf(fileID,'\n');
fprintf(fileID,'selected: \n');

if(~isempty(myi))
    clear good1 a2 f2 fit2
    
    x1=data(:,myi(1));
    good1=~isnan(x1) & ~isnan(gpp) & ~isnan(rs);
    
    % compute the starting values for nonlinear fit
    gpp2=gpp(good1)./(a1(1)*a1(2).*rs(good1)./(a1(1)+a1(2).*rs(good1)));
    f2=regress(gpp2,[ones(size(x1(good1))) x1(good1)]);
    a2g=[a1; f2(1,1);f2(2,1)];
    
    % nonlinear fit
    modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2)';
    fit2=fitnlm([rs(good1),x1(good1)],gpp(good1),modelfun,a2g);
    a2=fit2.Coefficients.Estimate;
    pval=fit2.Coefficients.pValue;
    res=fit2.Residuals.Raw;
    varexpl=var(fit2.Fitted)/var(gpp(good1));
    AIC2=length(gpp(good1))*log(sum(res.^2))+2*length(a2);
    h=lillietest(res);
    
    % write fitted values and staistics on output file
    fprintf(fileID,['%s ',repmat('%g ',1,11),'\n'],cell2mat(dataset(:,myi(1)+1).Properties.VariableNames),a2,pval,varexpl,AIC2,h);
    
elseif(isempty(myi) && n==0) % at least one additional varible optimizing the AIC compared to step 1
    
    clear good1 a2 f2 fit2
    addedvar(addedvar==i1)=[]; % delete the selected variable from the array of candidates

    x1=data(:,i1);
    good1=~isnan(x1) & ~isnan(gpp) & ~isnan(rs);
    
    % compute the starting values for nonlinear fit
    gpp2=gpp(good1)./(a1(1)*a1(2).*rs(good1)./(a1(1)+a1(2).*rs(good1)));
    f2=regress(gpp2,[ones(size(x1(good1))) x1(good1)]);
    a2g=[a1; f2(1,1);f2(2,1)];
    
    % nonlinear fit
    modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2)';
    fit2=fitnlm([rs(good1),x1(good1)],gpp(good1),modelfun,a2g);
    a2=fit2.Coefficients.Estimate;
    pval=fit2.Coefficients.pValue;
    res=fit2.Residuals.Raw;
    varexpl=var(fit2.Fitted)/var(gpp(good1));
    AIC2=length(gpp(good1))*log(sum(res.^2))+2*length(a2);
    h=lillietest(res);
    
    % write fitted values and staistics on output file
    fprintf(fileID,['%s ',repmat('%g ',1,11),'\n'],cell2mat(dataset(:,i1+1).Properties.VariableNames),a2,pval,varexpl,AIC2,h);
    
else % no additional varible in the loop optimizing the AIC -> test DOY and stop
    
    clear good1 a2 f2 fit2
    i1=5;
    x1=data(:,i1);
    good1=~isnan(x1) & ~isnan(gpp) & ~isnan(rs);
    
    % compute the starting values for nonlinear fit
    gpp2=gpp(good1)./(a1(1)*a1(2).*rs(good1)./(a1(1)+a1(2).*rs(good1)));
    f2=regress(gpp2,[ones(size(x1(good1))) x1(good1)]);
    a2g=[a1; f2(1,1);f2(2,1)]; %guess
    
    % nonlinear fit
    modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2)';
    fit2=fitnlm([rs(good1),x1(good1)],gpp(good1),modelfun,a2g);
    a2=fit2.Coefficients.Estimate;
    pval=fit2.Coefficients.pValue;
    res=fit2.Residuals.Raw;
    varexpl=var(fit2.Fitted)/var(gpp(good1));
    AIC2=length(gpp(good1))*log(sum(res.^2))+2*length(a2);
    h=lillietest(res);
    
    if(AIC2<AIC1)
        % write fitted values and staistics on output file
        fprintf(fileID,['%s ',repmat('%g ',1,11),'\n'],cell2mat(dataset(:,i1+1).Properties.VariableNames),a2,pval,varexpl,AIC2,h);

        i2=4;
        x2=data(:,i2);
        good2=(~isnan(x2) & ~isnan(x1) & ~isnan(gpp) & ~isnan(rs));
    
        % compute the starting values for nonlinear fit
        gpp2=gpp(good2)./(a2(1)*a2(2).*rs(good2)./(a2(1)+a2(2).*rs(good2)))-a2(4).*x1(good2);
        f3=regress(gpp2,[ones(size(x2(good2))) x2(good2)]);
        a3g=[a2(1); a2(2); f3(1,1); a2(4); f3(2,1)];
    
        % nonlinear fit
        modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2+b5*x3)';
        fit3=fitnlm([rs(good2),x1(good2),x2(good2)],gpp(good2),modelfun,a3g);
        a3=fit3.Coefficients.Estimate;
        pval=fit3.Coefficients.pValue;
        res=fit3.Residuals.Raw;
        varexpl=var(fit3.Fitted)/var(gpp(good2));
        AIC3=length(gpp(good2))*log(sum(res.^2))+2*length(a3);
        h=lillietest(res);
    
        if(AIC3<AIC2)
            fit2=fit3;
            good1=good2;
            % write fitted values and staistics on output file
            fprintf(fileID,['%s ',repmat('%g ',1,13),'\n'],cell2mat(dataset(:,i2+1).Properties.VariableNames),a3,pval,varexpl,AIC3,h);
        end
        
        
    else
        
        clear good1 a2 f2 fit2
        i1=4;
        x1=data(:,i1);
        good1=~isnan(x1) & ~isnan(gpp) & ~isnan(rs);

        % compute the starting values for nonlinear fit
        gpp2=gpp(good1)./(a1(1)*a1(2).*rs(good1)./(a1(1)+a1(2).*rs(good1)));
        f2=regress(gpp2,[ones(size(x1(good1))) x1(good1)]);
        a2g=[a1; f2(1,1);f2(2,1)]; %guess

        % nonlinear fit
        modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2)';
        fit2=fitnlm([rs(good1),x1(good1)],gpp(good1),modelfun,a2g);
        a2=fit2.Coefficients.Estimate;
        pval=fit2.Coefficients.pValue;
        res=fit2.Residuals.Raw;
        varexpl=var(fit2.Fitted)/var(gpp(good1));
        AIC2=length(gpp(good1))*log(sum(res.^2))+2*length(a2);
        h=lillietest(res);
        
        if(AIC2<AIC1)
            % write fitted values and staistics on output file
            fprintf(fileID,['%s ',repmat('%g ',1,11),'\n'],cell2mat(dataset(:,i1+1).Properties.VariableNames),a2,pval,varexpl,AIC2,h);
        end
        
    end
    
    % plot results: measured vs fitted
    figure
    scatter(gpp(good1),fit2.Fitted,50,'r','filled')
    hold on
    plot(min(gpp):max(gpp),min(gpp):max(gpp),'k','LineWidth',1.5)
    set(gca,'FontSize',12)
    xlabel('measured GPP (\mu mol/m^2/s)','FontSize',15)
    ylabel('modeled GPP (\mu mol/m^2/s)','FontSize',15)

    fclose(fileID);
    
    return
end

 %% third step: GPP vs additional variable #2
 
fprintf(fileID,'\n');
fprintf(fileID,'third step \n');
fprintf(fileID,'variable F_0 alpha_0 a_0 a_1 a_2 pvals varexpl AIC h \n');

n=1;

if(isempty(myi))
 for i=addedvar
    
    clear good2 X a3 f3 fit3
    good2=(~isnan(data(:,i)) & ~isnan(x1) & ~isnan(gpp) & ~isnan(rs));
    
    %%guess
    gpp2=gpp(good2)./(a2(1)*a2(2).*rs(good2)./(a2(1)+a2(2).*rs(good2)))-a2(4).*x1(good2);
    f3=regress(gpp2,[ones(size(data(good2,i))) data(good2,i)]);
    a3g=[a2(1); a2(2); f3(1,1); a2(4); f3(2,1)];

    X=[rs(good2),x1(good2),data(good2,i)];
    
    %%fitnonlin
    modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2+b5*x3)';
    fit3=fitnlm(X,gpp(good2),modelfun,a3g);
    a3=fit3.Coefficients.Estimate;
    pval=fit3.Coefficients.pValue;
    res=fit3.Residuals.Raw;
    varexpl=var(fit3.Fitted)/var(gpp(good2));
    AIC=length(gpp(good2))*log(sum(res.^2))+2*length(a3);
    h=lillietest(res);
    
    fprintf(fileID,['%s ',repmat('%g ',1,13),'\n'],cell2mat(dataset(:,i+1).Properties.VariableNames),a3,pval,varexpl,AIC,h);
    
    if(n==1 && AIC<AIC2) 
        best=AIC;
        i2=i;
        n=0;
    elseif(n==0 && AIC<best) 
        best=AIC;
        i2=i;
    end
 end
end

fprintf(fileID,'\n');
fprintf(fileID,'selected: \n');

if(~isempty(myi))
    
    clear good2 a3 f3 fit3
    
    x2=data(:,myi(2));
    good2=(~isnan(x2) & ~isnan(x1) & ~isnan(gpp) & ~isnan(rs));
    
    % compute the starting values for nonlinear fit
    gpp2=gpp(good2)./(a2(1)*a2(2).*rs(good2)./(a2(1)+a2(2).*rs(good2)))-a2(4).*x1(good2);
    f3=regress(gpp2,[ones(size(x2(good2))) x2(good2)]);
    a3g=[a2(1); a2(2); f3(1,1); a2(4); f3(2,1)];
    
    % nonlinear fit
    modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2+b5*x3)';
    fit3=fitnlm([rs(good2),x1(good2),x2(good2)],gpp(good2),modelfun,a3g);
    a3=fit3.Coefficients.Estimate;
    pval=fit3.Coefficients.pValue;
    res=fit3.Residuals.Raw;
    varexpl=var(fit3.Fitted)/var(gpp(good2));
    AIC3=length(gpp(good2))*log(sum(res.^2))+2*length(a3);
    h=lillietest(res);
    
    % write fitted values and staistics on output file
    fprintf(fileID,['%s ',repmat('%g ',1,13),'\n'],cell2mat(dataset(:,myi(2)+1).Properties.VariableNames),a3,pval,varexpl,AIC3,h);
    
elseif(isempty(myi) && n==0)
    
    clear good2 a3 f3 fit3
    addedvar(addedvar==i2)=[]; % delete the selected variable from the array of candidates

    x2=data(:,i2);
    good2=(~isnan(x2) & ~isnan(x1) & ~isnan(gpp) & ~isnan(rs));
    
    % compute the starting values for nonlinear fit
    gpp2=gpp(good2)./(a2(1)*a2(2).*rs(good2)./(a2(1)+a2(2).*rs(good2)))-a2(4).*x1(good2);
    f3=regress(gpp2,[ones(size(x2(good2))) x2(good2)]);
    a3g=[a2(1); a2(2); f3(1,1); a2(4); f3(2,1)];
    
    % nonlinear fit
    modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2+b5*x3)';
    fit3=fitnlm([rs(good2),x1(good2),x2(good2)],gpp(good2),modelfun,a3g);
    a3=fit3.Coefficients.Estimate;
    pval=fit3.Coefficients.pValue;
    res=fit3.Residuals.Raw;
    varexpl=var(fit3.Fitted)/var(gpp(good2));
    AIC3=length(gpp(good2))*log(sum(res.^2))+2*length(a3);
    h=lillietest(res);
    
    % write fitted values and staistics on output file
    fprintf(fileID,['%s ',repmat('%g ',1,13),'\n'],cell2mat(dataset(:,i2+1).Properties.VariableNames),a3,pval,varexpl,AIC3,h);
else  
    
    clear good2 a3 f3 fit3
    i2=5;
    x2=data(:,i2);
    good2=(~isnan(x2) & ~isnan(x1) & ~isnan(gpp) & ~isnan(rs));
    
    % compute the starting values for nonlinear fit
    gpp2=gpp(good2)./(a2(1)*a2(2).*rs(good2)./(a2(1)+a2(2).*rs(good2)))-a2(4).*x1(good2);
    f3=regress(gpp2,[ones(size(x2(good2))) x2(good2)]);
    a3g=[a2(1); a2(2); f3(1,1); a2(4); f3(2,1)];
    
    % nonlinear fit
    modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2+b5*x3)';
    fit3=fitnlm([rs(good2),x1(good2),x2(good2)],gpp(good2),modelfun,a3g);
    a3=fit3.Coefficients.Estimate;
    pval=fit3.Coefficients.pValue;
    res=fit3.Residuals.Raw;
    varexpl=var(fit3.Fitted)/var(gpp(good2));
    AIC3=length(gpp(good2))*log(sum(res.^2))+2*length(a3);
    h=lillietest(res);
    
    if(AIC3<AIC2)
        % write fitted values and staistics on output file
        fprintf(fileID,['%s ',repmat('%g ',1,13),'\n'],cell2mat(dataset(:,i2+1).Properties.VariableNames),a3,pval,varexpl,AIC3,h);
        
        i3=4;
        x3=data(:,i3);
        good3=(~isnan(x3) & ~isnan(x2) & ~isnan(x1) & ~isnan(gpp) & ~isnan(rs));

        % compute the starting values for nonlinear fit
        gpp2=gpp(good3)./(a3(1)*a3(2).*rs(good3)./(a3(1)+a3(2).*rs(good3)))-a3(4).*x1(good3)-a3(5).*x2(good3);
        f4=regress(gpp2,[ones(size(x3(good3))) x3(good3)]);
        a4g=[a3(1); a3(2); f4(1,1); a3(4); a3(5); f4(2,1)];

        % nonlinear fit
        modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2+b5*x3+b6*x4)';
        fit4=fitnlm([rs(good3),x1(good3),x2(good3),x3(good3)],gpp(good3),modelfun,a4g);
        a4=fit4.Coefficients.Estimate;
        pval=fit4.Coefficients.pValue;
        res=fit4.Residuals.Raw;
        varexpl=var(fit4.Fitted)/var(gpp(good3));
        AIC4=length(gpp(good3))*log(sum(res.^2))+2*length(a4);
        h=lillietest(res); 
        
        if(AIC4<AIC3)
            fit3=fit4;
            good2=good3;
            % write fitted values and staistics on output file
            fprintf(fileID,['%s ',repmat('%g ',1,15),'\n'],cell2mat(dataset(:,i3+1).Properties.VariableNames),a4,pval,varexpl,AIC4,h);
        end
    
    
    else
        
        clear good2 a3 f3 fit3
        i2=4;
        x2=data(:,i2);
        good2=(~isnan(x2) & ~isnan(x1) & ~isnan(gpp) & ~isnan(rs));

        % compute the starting values for nonlinear fit
        gpp2=gpp(good2)./(a2(1)*a2(2).*rs(good2)./(a2(1)+a2(2).*rs(good2)))-a2(4).*x1(good2);
        f3=regress(gpp2,[ones(size(x2(good2))) x2(good2)]);
        a3g=[a2(1); a2(2); f3(1,1); a2(4); f3(2,1)];

        % nonlinear fit
        modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2+b5*x3)';
        fit3=fitnlm([rs(good2),x1(good2),x2(good2)],gpp(good2),modelfun,a3g);
        a3=fit3.Coefficients.Estimate;
        pval=fit3.Coefficients.pValue;
        res=fit3.Residuals.Raw;
        varexpl=var(fit3.Fitted)/var(gpp(good2));
        AIC3=length(gpp(good2))*log(sum(res.^2))+2*length(a3);
        h=lillietest(res); 
        
        if(AIC3<AIC2)
            % write fitted values and staistics on output file
            fprintf(fileID,['%s ',repmat('%g ',1,13),'\n'],cell2mat(dataset(:,i2+1).Properties.VariableNames),a3,pval,varexpl,AIC3,h);
        end
        
    end
    
    % plot results: measured vs fitted
    figure
    scatter(gpp(good2),fit3.Fitted,50,'r','filled')
    hold on
    plot(min(gpp):max(gpp),min(gpp):max(gpp),'k','LineWidth',1.5)
    set(gca,'FontSize',12)
    xlabel('measured GPP (\mu mol/m^2/s)','FontSize',15)
    ylabel('modeled GPP (\mu mol/m^2/s)','FontSize',15)
    
    fclose(fileID);
    
    return
end

%% measured vs fitted data

% plot results when DOY has not been tested (see optional): measured vs fitted 
figure
scatter(gpp(good2),fit3.Fitted,50,'r','filled')
hold on
plot(min(gpp):max(gpp),min(gpp):max(gpp),'k','LineWidth',1.5)
set(gca,'FontSize',12)
xlabel('measured GPP (\mu mol/m^2/s)','FontSize',15)
ylabel('modeled GPP (\mu mol/m^2/s)','FontSize',15)

fclose(fileID);
%% optional
% copy and past the third step here to test additional variable #3 
% changes 
% good2-> good3=(~isnan(data(:,i)) & ~isnan(x2) & ~isnan(x1) & ~isnan(gpp) & ~isnan(rs));
% gpp2=gpp(good3)./(a3(1)*a3(2).*rs(good3)./(a3(1)+a3(2).*rs(good3)))-a3(4).*x1(good3)-a3(5).*x1(good3);
% X=[rs(good3),x1(good3),x2(good3),data(good3,i)];
% modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2+b5*x3+b6*x4)';
% a3->a4 f3->f4 fit3->f4
% AIC2->AIC3; AIC3->AIC4

%% recommended: check for parameter stability
% good2=(~isnan(x2) & ~isnan(x1) & ~isnan(gpp) & ~isnan(rs));
% modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2+b5*x3)';
% fit3prime=fitnlm([rs(good2),x1(good2),x2(good2)],gpp(good2),modelfun,a3);
% a3prime=fit3prime.Coefficients.Estimate;
% pval=fit3prime.Coefficients.pValue;
% res=fit3prime.Residuals.Raw;
% varexpl=var(fit3prime.Fitted)/var(gpp(good2));
% AIC3=length(gpp(good2))*log(sum(res.^2))+2*length(a3prime);
% h=lillietest(res);

