% nnitrain19.m
%
% to train the Neural Network Inversion (nni) for
% ozone number density profile.
% Consolidation of nnitrain12.m / ref=p A3.153
% target is logProfile
%
% DiF / v1.0 / 6-Oct-2023
%=============================================================================
%
% computation of the standard deviation of the radiance profiles is added
% and saved in o3Aux (dry run without performing the NN training)
%
% DiF / v1.1 / 26-Oct-2023
%=============================================================================
% also computation of percentiles 2.5 & 97.5 % of the radiance profiles is added
% and saved in o3Aux (dry run without performing the NN training)
%
% DiF / v1.2 / 24-Nov-2023
%=============================================================================

iset=1; % to setup data
iOut=1; % to plot and to generate function
iterMax=200000;
chan=[1 2 3 4 5 6 7 8];nch=length(chan);

if iset==1

    load('ozSet12.mat','vmr'); % true ozone values
    z=[0:1:60]';nz=length(z);
    inorm=41;
    nc=25600;
%     npcChan=[5 7 9 12 13 14 15 16]; % [0.1 % on log(R)]
    npcChan=[8 9 9 14 18 19 20 21]; % [0.3 % on R]
    X=[];Y=[];Ai=[];
    
    YY=zeros(nc,nz,nch); % log-radiance training set
    YYY=zeros(nc,nz,nch); % normalized log-radiances
    YMoz=zeros(nz,nch); % channel mean log-radiance profiles
    dYMoz=zeros(nz,nch); % channel std of log-radiance profiles
    YMoz025=zeros(nz,nch); % channel 2.5 percentile of log-radiance profiles
    YMoz975=zeros(nz,nch); % channel 97.5 percentile of log-radiance profiles
    Uoz=zeros(nz,nz,nch); % channel PC basis 
    for ich=chan
        [X,Y,LOK]=buildBigXY21(ich,1:1:nz,nc,1,1); 
        YY(:,:,ich)=Y;
        YYY(:,:,ich)=(Y(:,:)-Y(:,inorm));
    end
    l=find(all(LOK ,2));
    nl=length(l);
    
    for ich=chan
        np=npcChan(ich);
        Ydum=squeeze(YYY(l,:,ich));
        Ym=mean(Ydum);
        dY=Ydum-Ym;
        C=dY'*dY;
        [dum,MU,U]=svd(C,0);
        A=dY*U;
%         Yp=Ym+A(:,1:np)*U(:,1:np)';
        YMoz(:,ich)=Ym;
        dYMoz(:,ich)=std(Ydum);
        Yperc=prctile(Ydum,[2.5 16 50 84 97.5]);
        YMoz025(:,ich)=Yperc(1,:);
        YMoz975(:,ich)=Yperc(5,:);
        Uoz(:,:,ich)=U;
        ApCh=A(:,1:np);
        Ai=[Ai ApCh];
    end
end


input=[Ai X(l,33:34)];
output=vmr(l,1:nz);

rng(0, "threefry");
gpurng(0, "threefry");

[dum nXi]=size(output);

x = input';
t = output';



trainFcn = 'trainscg';  % Scaled conj grad backpropagation. (for GPU !!!)
%     %trainFcn = 'trainlm';  % L-M
% 
% hid1 =60;
% hid2 =60;
% hid3=60;
hid1 =61;
hid2 =61;



net = fitnet([hid1,hid2],trainFcn); % two hiddenlayers NN
% net = fitnet([hid1,hid2,hid3],trainFcn); % three hiddenlayers NN

% 
%     % Choose Input and Output Pre/Post-Processing Functions
%     % For a list of all processing functions type: help nnprocess
net.input.processFcns = {'removeconstantrows','mapminmax'};
net.output.processFcns = {'removeconstantrows','mapminmax'};
% 
%     % Setup Division of Data for Training, Validation, Testing
%     % For a list of all data division functions type: help nndivision
net.divideFcn = 'divideblock';  % Divide data 
net.divideMode = 'sample';  % Divide up every sample
net.divideParam.trainRatio = 70/100;
net.divideParam.valRatio = 15/100; % validation
net.divideParam.testRatio = 15/100;
% 
%     % Choose a Performance Function
%     % For a list of all performance functions type: help nnperformance
net.performFcn = 'mse';  % Mean Squared Error
% 
%     % Choose Plot Functions
%     % For a list of all plot functions type: help nnplot
net.plotFcns = {'plotperform','plottrainstate','ploterrhist', ...
    'plotregression', 'plotfit'};
% 
net.trainParam.epochs=iterMax;
net.trainParam.max_Fail=iterMax;
net.trainParam.sigma=5e-5; % def= 5e-5
net.trainParam.lambda=5e-7; % def= 5e-7
net.trainParam.showWindow=0;
net.trainParam.showCommandLine=1;
tic
%     % Train the Network
[net,tr] = train(net,x,t,'useGPU','yes');
et=toc;
% 
%     % Test the Network
y = net(x);
e = gsubtract(t,y);
performance = perform(net,t,y);
%     % Recalculate Training, Validation and Test Performance
trainTargets = t .* tr.trainMask{1};
valTargets = t .* tr.valMask{1};
testTargets = t .* tr.testMask{1};
trP = perform(net,trainTargets,y);
valPerformance = perform(net,valTargets,y);
teP = perform(net,testTargets,y);
%     
    
%     performance 
itar=find(isfinite(testTargets(1,:)));
ytar=testTargets(1:nXi,itar);
mat=corrcoef(y(1:nXi,itar),ytar);
reg=mat(1,2);
resul(1:4)=[et trP teP reg]

% 
clf1;plotperform(tr)
clf2;plotregression(t,y)

if iOut==1 
    zl=z;
    yTrue=vmr(l,1:nz);
    yPred=y';
    dy=yPred-yTrue;
    
    % icr=1:100;for ic=icr,figure(5);clf;plot(yTrue(ic,:),zl,'b',yPred(ic,:),zl,'r');ylim([5 35]);grid minor;pause,end
    figure(3);clf;set(gcf,'pos',[1208         329        1277        1003])
    yTrueM=mean(yTrue);
    yTrueS=std(yTrue);
    yPredM=mean(yPred);
    yPredS=std(yPred);
    plot(yTrueM,zl,'ob',yPredM,zl,'r','lineW',2);grid minor
    hold on
    plot(yTrueM+yTrueS,zl,'ob',yTrueM-yTrueS,zl,'ob','lineW',1);
    plot(yPredM+yPredS,zl,'r',yPredM-yPredS,zl,'r','lineW',1);
    xlabel('O_3 vmr');ylabel('z [km]');
    title(['nearL  true(blue) and predicted (red) with std ('  int2str(length(l)) ' cases)'])
    print -dpng nniOz3a
    
    valper=[0.1586 0.5 0.8414];
    clear dsz
    for iz=1:nz
        dum=wpercentile(dy(:,iz),[],valper);
        ds2=0.5*(dum(3)-dum(1));
        dsz(iz)=ds2;
    end
    
    figure(4);
    title(['ds=[p(0.84)-p(0.16)]/2 mean(ds)=' num2str(mean(dsz),3)])
    plot(dsz,zl);grid on;xlabel('ds [ppm]');ylabel('z [km]')
    title('Ozone median retrieval error')
    print -dpng nniOz3b
    
    save ozAux3 z inorm YYY npcChan Uoz YMoz dYMoz YMoz025 YMoz975
    genFunction(net,'nniOz3');  

end


