clear all;
close all;
clc;
format long;
A = importdata('Outputs/ExcCDF-RawData.txt');
x1=A(:,1)';
x2=A(:,2)';
y=A(:,3)';

%%--------------3D plot of  depth-angle-CDF of excess -------------------------
figure(1)
scatter3(x1,x2,y,2,'filled');hold on;
axis([0 89 0 2 -0.0008 0.0008]);
view(40,20);
 xlabel('\theta (\pi/180)');  ylabel('Depth (¦Ìm)'); zlabel('Cumulant I-V excess');
 
 %%  ------------NN Training on CDF of Excess---------------------
inputData  = [x1;x2];                                          % Input: x1,x2 (Depth and inclination angle)
outputData = y;                                                  %Output: Cumulant Excess calculated by IM3D

%----------------NN fitting Setting-------------------
% The number of hidden nodes is set to 20,10,5.
% Traibr method is used for training.
% The transfer functions of hidden layer and output layer are Tansig and purelin, respectively. 

net=newff(inputData,outputData,[20,10,5],{'tansig','tansig','purelin'},'trainbr','learngdm','msereg');
  load('EXC_BP.mat')                                           %Please comment it out if you would like to train a new surrogate model
   net.trainparam.goal = 1e-15;                          %Traning objective function: mean-square error 
   net.trainparam.show = 100;                          
   net.trainparam.epochs = 50000;                     % maximum number of iterations
   net.trainparam.min_grad=1e-15;
   net.trainParam.Ir=0.1;                                      % Leaning rate 
   net.trainParam.max_fail=1000;                        % Validation 

zz_net=(0.0:0.005:2.085);
dz_net=zz_net(2)-zz_net(1);
Ang_net=0:1:89;
[net,b]=train(net,[x1;x2],y)
[xx,yy]=meshgrid(Ang_net,zz_net);
xsim=xx(:)'
ysim=yy(:)'
 simout=sim(net,[xsim;ysim])
 D=simout;

figure(2);                                                            % Comparion of NN-trained CDF and CDF raw data
scatter3(x1,x2,y,2,'filled','r');hold on;
scatter3(xsim,ysim,D,2,'filled','b');     
%
xlabel('\theta');  ylabel('Depth (¦Ìm)');  zlabel('Net excess CDF distributation');
axis([0 89 0 2 -0.0008 0.0008]);

%save EXC_BP net

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Data (depth profiles) Checking: NN-fitted Excess vs. non-fitted data
 %%----------------------Depth profiles of excesses, 2D plot-------------------------
dep=(0.045:0.03:2.085); length(dep);
dz=dep(2)-dep(1);
Ang=0:1:89;  length(Ang);

figure(3);  
Ang0_model=0:1:89;                                        %Choose angle=0,10,20.... for comparion            
for ii=1:1:length(Ang0_model)
    a=Ang0_model(ii) ;
    postheta=find(abs(Ang_net-Ang0_model(ii))<=1e-4);
    for jj=1:length(zz_net) 
        tt=jj+length(zz_net)*(postheta-1);
        test_net(jj)=simout(tt) ;
    end 
    exc_net=gradient(test_net)./gradient(zz_net);
    dep_net=zz_net+0.5*dz;                                 % taking the derivative of CDF with respect to the ion depth
    plot(dep_net,exc_net,'-','Linewidth',0.1) ; hold on;
end  
%

dep=(0.045:0.03:2.085); length(dep);                 %Depth calculated by IM3D
dz=dep(2)-dep(1);
B = importdata('Outputs/ExcPDF-RawData.txt');              % Excess PDF raw data calculated by IM3D
x1r=B(:,1)';
x2r=B(:,2)';
yr=B(:,3)';
for ii=1:1:length(Ang0_model)
    a=Ang0_model(ii);  
    postheta=find(abs(Ang_net-Ang0_model(ii))<=1e-4);
    for jj=1:length(dep) 
        tt=jj+length(dep)*(postheta-1);
        test(jj)=yr(tt)  ;     
    end 
    exc0=test        
  %  plot(dep,exc0,'o'); hold on;
end  
legend('0','10','20','30','40','50','60','70','80','Location','SouthEast')
axis([0 2 -0.005 0.003]);
grid on;
xlabel('Depth (¦Ìm)');  ylabel('V-I excess PDF distributation');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% -------------------------Output NN-fitted Excess----------------------------
 %  (1) Output NN-fitted CDF Excess into file
fid_Ecdf_fitted= fopen('Outputs/ExcCDF-NNFitted.txt','w')
for j=1:length(xsim)
fprintf(fid_Ecdf_fitted,'%20.15f\t',xsim(j),ysim(j),simout(j));
fprintf(fid_Ecdf_fitted,'\n');
end
fclose(fid_Ecdf_fitted);
% (2) PDF obtained by taking the derivative of CDF with respect to the ion depth
fid_Epdf_fitted= fopen('Outputs/ExcPDF-NNFitted.txt','w')
Dpdf=gradient(simout)./gradient(ysim);
ypdf=ysim+0.5*dz; 
for j=1:1:length(xsim)
fprintf(fid_Epdf_fitted,'%20.15f\t',xsim(j),ypdf(j),Dpdf(j));
fprintf(fid_Epdf_fitted,'\n');
end
fclose(fid_Epdf_fitted);  

%  (3) Output NN-fitted total excess to a file for further COST evaluation
[xx,yy]=meshgrid(Ang_net,zz_net);
pdf_mat=reshape( Dpdf, size(xx) );
fid = fopen('Outputs/total-excess.txt','w')
[b1,b2] =size (pdf_mat);

for i=1:b1
    for j=1:b2          
       fprintf(fid,'%30.15f',pdf_mat(i,j)');  
    end  
     fprintf(fid,'\n');  
end  
fclose(fid);  
