%This script creates an ensemble of Artificial Neural Networks (ANNs) 
%to interpolate delta15N of nitrate observations globally using ocean 
%datasets from the World Ocean Atlas 2009 

%Author: Aaron Bagnell
%Email: abagnell@ucsb.edu
%Release: 1
%Release Date: May 2019

%Load d15NO3 ANN Input Data
load PR19_d15NO3_ANN_anclim_inputs.mat

%Find Grid Cells Where There is d15NO3 observations 
ii = find(d15NO3_obs > -9.999 & temperature > -9.999 & nitrate > -9.999 & chl_a > -9.999);

%Select Input/Output Feautures
inputs = [temperature(ii) salinity(ii) nitrate(ii) oxygen(ii) phosphate(ii)...
 log10(chl_a(ii))];

outputs = d15NO3_obs(ii);

%Create Matrix of Features to Generalize ANN Solution
ii2 = find(temperature> -9.999 & nitrate > -9.999 & chl_a > -9.999);
sitnew = [temperature(ii2) salinity(ii2) nitrate(ii2) oxygen(ii2) phosphate(ii2)...
 log10(chl_a(ii2))];

%Create Ensemble of ANNs
i = 0;
while i < 25 %Number of Ensemble Members
varfinal = 0*d15NO3_obs(:)+NaN;

     x = inputs';
     t = outputs';

     %Create Training and Validation Sets
     nl = size(inputs,1);
     ip = randperm(nl);
     x = x(:,ip);
     t = t(:,ip);

     train_pct = 0.9; %90 Percent of Data Is Used for Training, 10 Percent for Validation
     valid_pct = 1-train_pct;
     itrain = 1:floor(train_pct*nl);
     itest = itrain(end) + 1:nl;

     % Create a Fitting Network
     hiddenLayerSize = [25];
     trainFcn = 'trainbr'; %Backpropogation Algoritm is Levenberg-Marquardt
                           %with Bayesian Regularization
     net = fitnet(hiddenLayerSize,trainFcn);

     %Specify Transfer (Activation) Function for Each Layer
     for j = 1:length(hiddenLayerSize)
          net.layers{j}.transferFcn = 'tansig';
     end
     net.layers{length(hiddenLayerSize)+1}.transferFcn = 'tansig';

     % Setup Division of Data for Training, Validation, Testing
     net.divideFcn = 'divideind';
     net.divideParam.trainInd = itrain;
     net.divideParam.valInd   = itest;
     net.divideParam.testInd  = []; 

     %Set Number of Allowed Failed Training Iterations and Maximum Number of Iterations
     net.trainParam.max_fail = 5;
     net.trainParam.epochs = 1000;
    
     % Train the Network
     [net,tr] = train(net,x,t);

     varfinal(ii2) = net(sitnew');

     mdl = fitlm(outputs(ip(itest)),varfinal(ii(ip(itest))));
         
     %Allow Models Through that Meet Pass Criteria
     if mdl.Rsquared.ordinary >= 0.81
           
          i = i+1;
            
          varfinal_it(:,i) = varfinal;

          fprintf(sprintf('R^2 = %d & RMSE = %d for iteration %d\n',mdl.Rsquared.ordinary,mdl.RMSE,i))
            
     end

end

%Save Results from d15NO3 Ensemble
load PR19_d15NO3_ANN_anclim_inputs.mat lat lon z

d15NO3 = reshape(mean(varfinal_it,2),size(lat));
d15NO3(isnan(d15NO3)==1) = -9.999;
d15NO3_stdev = reshape(std(varfinal_it,[],2),size(lat));
d15NO3_stdev(isnan(d15NO3_stdev)==1) = -9.999;

save('PR19_d15NO3_ANN_anclim_outputs.mat','d15NO3','d15NO3_stdev','lat','lon','z')
