%Example script for computing random forests, sensitivities for Holder and
%Gnanadesikan Global Biogeochemical Cycles, (2023). 

predNames={'dfe','mld','nh4','no3','po4','rsds','sios','sos','tos','w50'};
varNamesLong={'Index','Latitude','Longitude','Month','Dissolved iron','Mixed Layer Depth','Ammonium','Nitrate',...
    'Phytoplankton Carbon','Phosphate','Shortwave','Silicate','Salinity','Temperature','Upwelling at 50m'};

load table_modis.asc%Add other files here as needed

j=1;%Full dataset saved as a structure.
final{j}=table_modis;
modelNames{j}='Observations-MODIS';
treenum(j)=200; % For other datasets use 50
ndims{1}=[360 180 12]; %x-y-t dimensions

jRange=1:1;% Alter if adding files

for j=jRange % Find predictor variables available for each model
    resp = [9]; % Response variable (phyc)
    npred=0;
    for k=5:15
     if (sum(final{j}(:,k)~=0)>0&k~=resp)
      npred=npred+1;
      pred{j}(:,npred)=k;
     end
    end 
    varPredNames{j}=varNamesLong(pred{j});
end


% Create training and testing datasets

for j = jRange
    p = 0.8; % Percentage of rows to select for training
    rng(123,'twister'); % rng stands for random number generation; makes the results reproducible (same thing each time)
    splitIdx = randperm(size(final{j},1)); % Generates numbers for random rows of the "final" variable to be split into training and testing datasets
    dataTraining{j} = final{j}(splitIdx(1:round(size(final{j},1).*p)),:); % Moves seleted rows into the dataTraining variable
    dataTesting{j} = final{j}(splitIdx(round(size(final{j},1).*p)+1:end),:); % Moves other selected rows into the dataTesting variable
end

%% Random Forests

tic
clear Mdl_rf
for j = jRange
    rng(123,'twister') % For reproducability
    Mdl_rf{j} = TreeBagger(treenum(j),dataTraining{j}(:,pred{j}),dataTraining{j}(:,resp),'SampleWithReplacement','off','InBagFraction',0.632,'OOBPredictorImportance','on','PredictorNames',varPredNames{j},'PredictorSelection','interaction-curvature','Method','regression','NumPrint',5); % RF model
    fprintf('Random forest model completed for %s \n',modelNames{j})
end

%Compute permutation (predpermute) and replacement sensitivity (predmed)
for j=jRange
     if ~isempty(final{j})
       predictphyc{j}=predict(Mdl_rf{j},final{j}(:,pred{j}));
       prednum=length(pred{j})
         for k=1:prednum
             predictors=final{j}(:,pred{j});
             predtest=dataTesting{j}(:,pred{j});
             predictors(:,k)=prctile(final{j}(:,pred{j}(k)),50);
             predtest(:,k)=dataTesting{j}(randperm(length(dataTesting{j}(:,k))),pred{j}(k));
             predmed{j}(:,pred{j}(k))=predict(Mdl_rf{j},predictors);
             predpermute{j}(:,pred{j}(k))=predict(Mdl_rf{j},predtest);
             [j k]
         end
     end
end

%Regrid prediction/replacement sensitivity to grid
for j=jRange
    predgrid{j}=nan(ndims{j});
    predgrid{j}(final{j}(:,1))=predictphyc{j}(:);
    for k=1:length(pred{j})
       predmed_grid{j,k}=nan(ndims{j});
       predmed_grid{j,k}(final{j}(:,1))=predmed{j}(:,pred{j}(k));
    end
end

%Compute median sensitivity: Store in sensmed
for j=jRange
    for k=1:length(pred{j})
     meddata=ones(199,1)*prctile(final{j}(:,pred{j}),50);
     meddata(:,k)=prctile(final{j}(:,pred{j}(k)),(0.5:0.5:99.5));
     sensmed{j}(:,k)=predict(Mdl_rf{j},meddata);
    end
end
sensaxis=prctile(predcesm2,(0.5:0.5:99.5));