function Rsl=RV(Rsl,Prm,Dat)
%**************************************************************************
%Illustration of non-stationary extreme value modelling in MATLAB accompanying  
%Threshold modeling of nonstationary extremes
%Paul Northrop, Philip Jonathan and David Randell
%in 
%Extreme Value Modeling and Risk Analysis: Methods and Applications
%Dipak Dey and Jun Yan (Editors)
%2015
%**************************************************************************
%
%**************************************************************************
%Return value simulation
%Rsl          1 x 1 Rsl structure
%Prm          1 x 1 Prm structure
%Dat          1 x 1 Dat structure
%
%Rsl          1 x 1 Rsl structure
%**************************************************************************

%Define sector boundaries
BinEdg=linspace(Prm.BssDmn(1),Prm.BssDmn(2),Prm.RVnSct+1);
BinEdg(1)=-Inf;
BinEdg(end)=Inf;

%Allocate memory for return value estimtes
Rsl.RV=NaN(Prm.RVnSct,Prm.RVnRls*Prm.BSnRsm);

fprintf('Generating %g Realisations from Extreme value model\n',Prm.RVnRls)

%Run return value diagnostic only if return period factor is unity
if Prm.RVRtrPrdFct==1;
    DgnOn=true;
    p=linspace(0,1,100);
    Rsl.RVDgn=NaN(Prm.RVnSct+1,100,Prm.RVnRls*Prm.BSnRsm);
else
    DgnOn=false;
end;

cI=0;
for iRls=1:Prm.RVnRls
    
    fprintf('.')
    if rem(iRls,20)==0
        fprintf('\n')
    end
    
    for iB=1:Prm.BSnRsm
        
        cI=cI+1;

        %generate number of storms severity from Poisson
        C=poissrnd(Rsl.Pss.MLE(:,iB)*diff(Rsl.X(1:2))/Rsl.Pss.DltX*Prm.RVRtrPrdFct);
        nStr=sum(C);  %number of storms for each realisation
        
        %populate covariate value for each storm
        count=0;
        X=NaN(nStr,1);
        for iStr=1:length(Rsl.X)  
            I=(count+1):(count+C(iStr));
            count=count+C(iStr);
            X(I)=iStr;
        end
        
        %generate storm sevirity from GP distirbutino
        Y=gprnd(Rsl.GP.Shp.MLE(X,iB),Rsl.GP.Scl.MLE(X,iB),Rsl.QR.MLE(X,iB));
        
        %accummulate return value estimates
        [~,IBin]=histc(X,BinEdg); %find covariate bin
        Rsl.RV(:,cI)=accumarray(IBin,Y,[],@max); %find max per bin
        
        if DgnOn; %diagnostic plot if necessary
            t1=cell2mat(accumarray(IBin,Y,[],@(x){quantile(x,p)}));           
            Rsl.RVDgn(:,:,cI)=[t1;quantile(Y,p)];
        end;
        
    end;
end;
fprintf('\n')

%Omni-directional return value estimate
Omni_RV=max(Rsl.RV,[],1);

%Prepare figure labels
BinEdg(1)=Prm.BssDmn(1);
BinEdg(end)=Prm.BssDmn(2);
tLbl=cell(Prm.RVnSct,1);
for i=1:Prm.RVnSct
    tLbl{i}=sprintf('[%g,%g]',BinEdg(i),BinEdg(i+1));
end

%Diagnostic plot if necessary
if DgnOn; 
    IExc=Dat.Y>Rsl.QR.ThrDat;
    XExc=Dat.X(IExc);
    YExc=Dat.Y(IExc);
    [~,IBinDat]=histc(XExc,BinEdg); %find covariate bin
    t1=cell2mat(accumarray(IBinDat,YExc,[],@(x){quantile(x,p)}));
    RVDgnDat=[t1;quantile(YExc,p)];
    
    figure(4)
    clf
    qRV=quantile(Rsl.RVDgn,[0.025,0.5,0.975],3);   
    C=jet(Prm.RVnSct);    
    nPlt=ceil(sqrt(Prm.RVnSct+1));
    subplot(nPlt,nPlt,1)
    plot(qRV(end,:,1),p,'r--','linewidth',2)
    hold on
    plot(qRV(end,:,2),p,'r-','linewidth',2)
    plot(qRV(end,:,3),p,'r--','linewidth',2)    
    plot(RVDgnDat(end,:),p,'k-','linewidth',2)
    title('Omni')
    for iS=1:Prm.RVnSct
        subplot(nPlt,nPlt,iS+1)
        plot(qRV(iS,:,1),p,'r--','linewidth',2)
        hold on
        plot(qRV(iS,:,2),p,'r-','linewidth',2)                
        plot(qRV(iS,:,3),p,'r--','linewidth',2)  
        plot(RVDgnDat(iS,:),p,'k-','linewidth',2)
        title(tLbl(iS))
        if iS+1==nPlt^2-nPlt+1;
            xlabel('Value');
            ylabel('CDF');
        end;
    end    
end

%Return value plot
figure(5)
clf
C=jet(Prm.RVnSct);
I=(1:Prm.RVnRls*Prm.BSnRsm)/Prm.RVnRls*Prm.BSnRsm;
for iS=1:Prm.RVnSct
    plot(sort(Rsl.RV(iS,:)),I,'-','color',C(iS,:),'linewidth',2)
    hold on
end
plot(sort(Omni_RV),I,'k-','linewidth',2)
xlabel('Value');
ylabel('CDF');
legend([tLbl;{'Omni'}])

return;
