%%These codes are provided to support paper 'Quantifying soil gaseous nitrogen losses from nitrification and denitrification based on isotope benchmarking model' submitted to Global Biogeochemical Cycles
%% developed by Wenxin Ba and Haoming Yu
%%This function script is supposed to run in Matlab ver. R2020a
%%Compatibility to other version of Matlab is not tested.

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Part 1 Global mapping of soil δ15N in croplands

%%%%%%%Constructing a Random Forest Model
clearvars;
close all;
data = readtable('E:\A_data\soil\Data_train2.3.csv');
Data_o = data(:, [1,5,6,7,8,9,10]);
Data_train=table2array(Data_o)
rng(222);
ind = randsample(1:2, size(Data_train, 1), true, [0.8, 0.2]);
df = Data_train(ind == 1, :);
df_test = Data_train(ind == 2, :);
tree = df;
%  Divide the training set and testing set
PN = df(:, 2:end);
TN = df(:, 1);
PM = df_test(:, 2:end);
TM = df_test(:, 1);
PN=PN';
[pn, ps_input] = mapminmax(PN, 0, 1);
pn=pn';
PM=PM';
pm = mapminmax('apply', PM, ps_input);
pm=pm';
TN=TN';
[tn, ps_output] = mapminmax(TN, 0, 1);
tn=tn';
%%  Model parameter setting and training model
trees = 600; 
leaf  = 3; 
OOBPrediction = 'on';  
OOBPredictorImportance = 'on';
Method = 'regression';
treenet = TreeBagger(trees, pn, tn, 'OOBPredictorImportance', OOBPredictorImportance,...
      'Method', Method, 'OOBPrediction', OOBPrediction, 'minleaf', leaf);
importance = treenet.OOBPermutedPredictorDeltaError;
pyuce = predict(treenet, pm );
Pyuce = mapminmax('reverse', pyuce, ps_output);
Pyuce =Pyuce';

%%%%%%% Model application
load(strcat('E:\mean_temperature.mat'),'temperature');
load(strcat('E:\mean_perirr.mat'),'perirr');
load(strcat('E:\mean_silt2.mat'),'silt');

I=ones(2160,4320);
yearly_temperatureexp=exp((temperature+I.*273.15)./100);
siltexp=exp(silt./100);
BNFln=log(BNF);
perirrln=log(perirr);
C=[1,5,6,8,9,11,13];
BNFln(isinf(BNFln))=nan;
perirrln(isinf(perirrln))=nan;
siltexp(silt<0)=nan;

data_input2=[reshape(yearly_temperatureexp,1,[]);reshape(siltexp,1,[]);reshape(BNFln,1,[]);reshape(fer_15,1,[]);reshape(perirrln,1,[]);reshape(noy,1,[])];%
nonEmptyRows = find(~any(isnan(data_input2), 1));
data_input20 = data_input2(:,~any(isnan(data_input2), 1));
[data_input0, ps_input] = mapminmax(data_input20 , 0, 1);
prec_data=data_input0';
pyuce = predict(treenet, prec_data);
Pyuce = mapminmax('reverse', pyuce, ps_output);
soil_15N=zeros(2160,4320);
soil_15N(nonEmptyRows)=Pyuce';

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Part 2 Isotope-benchmarking model for cropland ecosystems
U_vol=20;
U_de=13;
U_ni=47.5;
U_netuptake=0.432;
BNF15=-2;
DEP15=0;

N_input=TN_ha_BNF2+fertilizer+yearly_nhx+yearly_noy;
N15_input=(fer_15+BNF15.*TN_ha_BNF2+DEP15.*(yearly_nhx+yearly_noy))./((TN_ha_BNF2+fertilizer+yearly_nhx+yearly_noy));
        a=1.2*I./(I+exp(-0.04*(swvl*100-66.5*I)));
        b=0.76*I./(I+exp(0.07*(swvl*100-38.7*I)))+0.32*I;
        c=0.77*I-0.2*I./(I+exp(0.07*(swvl*100-59.7*I)));
        fnit=(I./((I-a)./a+(I-b)./b+I)).*c;
        fdenit=I-fnit;
U_gas=(fnit.*U_ni+fdenit.*U_de);
f_netuptake=TN_ha_uptake./N_input;
f_vol=NH3_sum./N_input;
f_gas=(soil_15N-N15_input-f_vol.*((U_vol).*I)-f_netuptake.*((U_netuptake).*I))./((U_gas).*I);
N_gas=f_gas.*N_input;

%%%%function
function [f_gas,N_gas] = Denitrification_flux(input_data,input_params)
I=ones(length(input_data(:,1)),1);
soil_15N=input_data(:,1);
TN_ha_BNF2=input_data(:,2);
residue=input_data(:,3);
yearly_nhx=input_data(:,4);
yearly_noy=input_data(:,5);
fertilizer=input_data(:,6);
rfer_15=input_params(:,7);
TN_ha_uptake=input_data(:,7);
NH3_sum=input_data(:,8);
swvl=input_data(:,9);
area_crop_sum=input_data(:,10);
ind_noarea=find(area_crop_sum==0);

fertilizer(isnan(fertilizer))=0;
fertilizer(ind_noarea)=nan;
rfer_15(isnan(rfer_15))=0;
rfer_15(ind_noarea)=nan;
rfer_15(isinf(rfer_15))=nan;
residue(isnan(residue))=0;
residue(ind_noarea)=nan;

U_de=input_params(:,1);
U_ni=input_params(:,2);
U_vol=input_params(:,3);
U_netuptake=input_params(:,4);
BNF15=input_params(:,5);
DEP15=input_params(:,6);


fer_15=rfer_15.*fertilizer;

N_input=TN_ha_BNF2+fertilizer+yearly_nhx+yearly_noy;
N15_input=(fer_15+BNF15.*TN_ha_BNF2+DEP15.*(yearly_nhx+yearly_noy))./((TN_ha_BNF2+fertilizer+yearly_nhx+yearly_noy));
        a=1.2*I./(I+exp(-0.04*(swvl*100-66.5*I)));
        b=0.76*I./(I+exp(0.07*(swvl*100-38.7*I)))+0.32*I;
        c=0.77*I-0.2*I./(I+exp(0.07*(swvl*100-59.7*I)));
        fnit=(I./((I-a)./a+(I-b)./b+I)).*c;
        fdenit=I-fnit;
U_gas=(fnit.*U_ni+fdenit.*U_de);
        max(max(U_gas))
        min(min(U_gas))
f_netuptake=TN_ha_uptake./N_input;
f_vol=NH3_sum./N_input;
f_vol(isinf(f_vol)) = nan;
f_gas=(soil_15N-N15_input-f_vol.*((U_vol).*I)-f_netuptake.*((U_netuptake).*I))./((U_gas).*I);
N_gas=f_gas.*N_input;

end

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Part 3 Uncertainty estimation (Monte Carlo method)
clear;
parameter= xlsread('E:\data_se_name_6.xlsx');
load(strcat('D:\data\raw data\rice_sowing_ave.mat'),'area');
area(isnan(area))=0;
country_ID=imread('D:\data\raw data\country_id_m.tif');
area(country_ID<0)=0;
area_crop_sum=area;
load(strcat('D:\result\residu_China.mat'),'residue');
I=ones(2160,4320);
load(strcat('D:\result\rfer_15_China.mat'),'rfer_15');
load(strcat('D:\result\NH3_sum_China.mat'),'NH3_sum');
load(strcat('D:\data\rice_mean_BNF.mat'),'BNF');
TN_ha_BNF2=BNF;
load(strcat('D:\data\mean_nhx.mat'),'nhx');
load(strcat('D:\data\mean_noy.mat'),'noy');
yearly_nhx=nhx;
yearly_noy=noy;
load(strcat('D:\result\soil_15N_rice_end.mat'),'soil_15N');
load(strcat('D:\result\TN_ha_uptake_China.mat'),'TN_ha_uptake');
load(strcat('D:\data\raw data\mean_WFPS_BD.mat'),'WFPS');
swvl=WFPS;
load(strcat('D:\result\fertilizer_end_China.mat'),'fertilizer');
fertilizer_data_ha=fertilizer;

area_crop_sum(area_crop_sum==0)=nan;
input_data0 = [reshape(soil_15N,[],1),reshape(TN_ha_BNF2,[],1),reshape(residue,[],1),reshape(yearly_nhx,[],1),reshape(yearly_noy,[],1),reshape(fertilizer,[],1),reshape(rfer_15,[],1),reshape(TN_ha_uptake,[],1),reshape(NH3_sum,[],1),reshape(swvl,[],1),reshape(area_crop_sum,[],1)];
nonEmptyRows = find(~any(isnan(input_data0), 2));
input_data1 = input_data0(~any(isnan(input_data0), 2), :);
input_data =input_data1;
Q=ones(length(input_data),1)
U_vol=20;
U_de=13;
U_ni=47.5;
U_netuptake=0.432;
BNF15=-2;
DEP15=0;

input_params = [U_de,U_ni,U_vol,U_netuptake,BNF15,DEP15,input_data(:,7)];
input_data(:,7)=[];
[fgas0,N_gas,fu_ratio]=Denitrification_flux(input_data,input_params);

num_samples = 1000;
% Initialize result array
results = zeros(length(input_data(:,1)),3,3);
    results_inputdata=zeros(length(input_data(:,1)),num_samples);
    results_params=zeros(length(input_data(:,1)),num_samples);
    results_sum=zeros(length(input_data(:,1)),num_samples);
    gas_result=zeros(6,num_samples);
for i = 1:num_samples%Randomly generate 1000 parameters
        
        input_params0=input_params;
        perturbed_data=input_data;
        for j=1:6
        parameter_sample =normrnd(parameter(j,1), parameter(j,5), size(input_data(:,1)));%基于均值，标准差生成正态分布的随机序列
        input_params0(:,j)=parameter_sample;
        end
        
        for j = 1:10
        random_data = unifrnd (-1, 1, size(input_data(:,1)));
        perturbed_data(:,j) = perturbed_data(:,j) + perturbed_data(:,j).* random_data.* 0.1;  
        end

        random_data = unifrnd (-1, 1, size(input_data(:,1)));
        input_params0(:,7)= input_params0(:,7) + input_params0(:,7).* random_data.* 0.1;
        [perturbed_f_gas,N_gas]=Denitrification_flux(input_data, input_params0);
        results_params(:,i) = (fgas- fgas0)./fgas0;%Calculate the relative change
        
        [perturbed_f_gas,N_gas] = Denitrification_flux(perturbed_data, input_params);
        results_inputdata(:,i) = (perturbed_f_gas - fgas0)./fgas0;%Calculate the relative change
        
        [perturbed_f_gas,N_gas] = Denitrification_flux(perturbed_data, input_params0);
        results_sum(:,i) = (perturbed_f_gas - fgas0)./fgas0;%Calculate the relative change

end        

U=ones(1,1000)   
for i = 1:length(input_data(:,1))
        
    [results(i,1,1),results(i,1,2)]=mystd(U,results_params(i, :));
    [results(i,2,1),results(i,2,2)]=mystd(U,results_inputdata(i, :));
    [results(i,3,1),results(i,3,2)]=mystd(U,results_sum(i, :));
    results(i,1,2)= prctile(results_params(i, :), 97.5);
    results(i,1,3)= prctile(results_params(i, :), 2.5);
    
    results(i,2,2)= prctile(results_inputdata(i, :), 97.5);
    results(i,2,3)= prctile(results_inputdata(i, :), 2.5);
    
    results(i,3,2)= prctile(results_sum(i, :), 97.5);
    results(i,3,3)= prctile(results_sum(i, :),2.5);

end
    
    Uncertainty = zeros(2160,4320,3,3);

    for i=1:3
        for j=1:3
            results0= zeros(9331200, 1);
            results0(nonEmptyRows)=results(:,i,j);
            Uncertainty(:,:,i,j)=reshape(results0,2160,[]);
        end
    end
save(strcat('D:\Uncentainty\Sum_Uncertainty_std_rice_data1.mat'),'Uncertainty');
