%   [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd]=Nonpar(t,M,iop,Mmin)
%
% BASED ON MAGNITUDE SAMPLE DATA M DETERMINES THE ROUND-OFF INTERVAL LENGTH
% OF THE MAGNITUDE DATA - eps, THE SMOOTHING FACTOR - h, CONSTRUCTS 
% THE BACKGROUND SAMPLE - xx AND CALCULATES THE WEIGHTING FACTORS - ambd 
% FOR A USE OF THE NONPARAMETRIC ADAPTATIVE KERNEL ESTIMATORS OF MAGNITUDE 
% DISTRIBUTION.
%
% !! THIS FUNCTION MUST BE EXECUTED AT START-UP OF THE UNBOUNDED 
%       NON-PARAMETRIC HAZARD ESTIMATION MODE !!
%
% AUTHOR: S. Lasocki 06/2014 within IS-EPOS project.
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative 
% to estimating the magnitude distribution functions. The smoothing factor 
% h, is estimated using the least-squares cross-validation for the Gaussian
% kernel function. The final form of the kernel is the adaptive kernel.
% In order to avoid repetitions, which cannot appear in a sample when the
% kernel estimators are used, the magnitude sample data are randomized
% within the magnitude round-off interval. The round-off interval length -
% eps is the least non-zero difference between sample data or 0.1 is the
% least difference if greater than 0.1. The randomization is done
% assuming exponential distribution of m in [m0-eps/2, m0+eps/2], where m0
% is the sample data point and eps is the length of roud-off inteval. The 
% shape parameter of the exponential distribution is estimated from the whole
% data sample assuming the exponential distribution. The background sample 
% - xx comprises the randomized values of magnitude doubled symmetrically 
% with respect to the value Mmin-eps/2: length(xx)=2*length(M). Weigthing 
% factors row vector for the adaptive kernel is of the same size as xx. 
% See: the references below for a more comprehensive description.
% 
% This is a beta version of the program. Further developments are foreseen.
%
% REFERENCES:
%Silverman B.W. (1986) Density Estimation for Statistics and Data Analysis, 
%   Chapman and Hall, London 
%Kijko A., Lasocki S., Graham G. (2001) Pure appl. geophys. 158, 1655-1665
%Lasocki S., Orlecka-Sikora B. (2008) Tectonophysics 456, 28-37
%
% INPUT:
%   t - vector of earthquake occurrence times
%   M - vector of earthquake magnitudes (sample data)
%   iop - determines the used unit of time. iop=0 - 'day', iop=1 - 'month', 
%       iop=2 - 'year'
%   Mmin - lower bound of the distribution - catalog completeness level
%
% OUTPUT
%   lamb_all - mean activity rate for all events
%   lamb - mean activity rate for events >= Mmin
%   lamb_err - error paramter on the number of events >=Mmin. lamb_err=0
%       for 50 or more events >=Mmin and the parameter estimation is
%       continued, lamb_err=1 otherwise, all output paramters except 
%       lamb_all and lamb are set to zero and the function execution is 
%       terminated.  
%   unit - string with name of time unit used ('year' or 'month' or 'day').
%   eps - length of round-off interval of magnitudes.  
%   ierr - h-convergence indicator. ierr=0 if the estimation procedure of 
%   the optimal smoothing factor has converged (the zero of the h functional 
%   has been found, ierr=1 when multiple zeros of h functional were 
%   encountered - the largest h is accepted, ierr = 2 when h functional did  
%   not zeroe - the approximate h value is taken.
%   h - kernel smoothing factor.
%   xx - the background sample for the nonparametric estimators of magnitude 
%   distribution
%   ambd - the weigthing factors for the adaptive kernel
%
% LICENSE
%     This file is a part of the IS-EPOS e-PLATFORM.
%
%     This is free software: you can redistribute it and/or modify it under 
%     the terms of the GNU General Public License as published by the Free 
%     Software Foundation, either version 3 of the License, or 
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 

function [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd]=...
            Nonpar(t,M,iop,Mmin)

lamb_err=0;
n=length(M);
t1=t(1);
for i=1:n
    if M(i)>=Mmin; break; end
    t1=t(i+1);
end
    t2=t(n);
for i=n:1
    if M(i)>=Mmin; break; end
    t2=t(i-1);
end
nn=0;
for i=1:n
    if M(i)>=Mmin
        nn=nn+1;
    end
end
    
if iop==0
    lamb_all=n/round(t(n)-t(1));
    lamb=nn/round(t2-t1);
    unit='day';
elseif iop==1
lamb_all=30*n/(t(n)-t(1));      % K20OCT2014
lamb=30*nn/(t2-t1);               % K20OCT2014
    unit='month';
else
lamb_all=365*n/(t(n)-t(1));       % K20OCT2014
lamb=365*nn/(t2-t1);                % K20OCT2014
    unit='year';
end

if nn<50
    eps=0;ierr=0;h=0;
    lamb_err=1;
    return;
end

eps=magn_accur(M);
n=0;
for i=1:length(M)
    if M(i)>=Mmin;
        n=n+1;
        x(n)=M(i);
    end
end
x=sort(x)';
beta=1/(mean(x)-Mmin+eps/2);
[xx]=korekta(x,Mmin,eps,beta);
xx=sort(xx);
clear x;
xx = podwajanie(xx,Mmin-eps/2);
[h,ierr]=hopt(xx);
[ambd]=scaling(xx,h);
end

function [m_corr]=korekta(m,Mmin,eps,beta)

% RANDOMIZATION OF MAGNITUDE WITHIN THE ACCURACY INTERVAL
%
% m - input vector of magnitudes 
% Mmin - catalog completeness level
% eps - accuracy of magnitude
% beta - the parameter of the unbounded exponential distribution
%
% m_corr - vector of randomized magnitudes
%
F1=1-exp(-beta*(m-Mmin-0.5*eps));
F2=1-exp(-beta*(m-Mmin+0.5*eps));
u=rand(size(m));
w=u.*(F2-F1)+F1;
m_corr=Mmin-log(1-w)./beta;
end

function x2 = podwajanie(x,x0)

% DOUBLES THE SAMPLE

% If the sample x(i) is is truncated from the left hand side and belongs 
% to the interval [x0,inf) or it is truncated from the right hand side and
% belongs to the interval (-inf,x0]
%   then the doubled sample is [-x(i)+2x0,x(i)]
% x - is the column data vector
% x2 - is the column vector of data doubled and sorted in the ascending
% order

x2=[-x+2*x0
    x];
x2=sort(x2);
end

function [h,ierr]=hopt(x)

%Estimation of the optimal smoothing factor by means of the least squares
%method
% x - column data vector
% The result is an optimal smoothing factor
% ierr=0 - convergence, ierr=1 - multiple h, ierr=2 - approximate h is used 
% The function calls the procedure FZERO for the function 'funct'
% NEW VERSION 2 - without a square matrix. Also equipped with extra zeros
% search

% MODIFIED JUNE 2014

ierr=0;
n=length(x);
x=sort(x);
interval=[0.000001 2*std(x)/n^0.2];
x1=funct(interval(1),x);
x2=funct(interval(2),x);
if x1*x2<0
   [hh(1),fval,exitflag]=fzero(@funct,interval,[],x);

% Extra zeros search
    jj=1;
    for kk=2:7
        interval(1)=1.1*hh(jj);
        interval(2)=interval(1)+(kk-1)*hh(jj);
        x1=funct(interval(1),x);
        x2=funct(interval(2),x);
        if x1*x2<0
            jj=jj+1;
            [hh(jj),fval,exitflag]=fzero(@funct,interval,[],x);
        end
    end
    if jj>1;ierr=1;end
    h=max(hh);

   if exitflag==1;return;end

end
h=0.891836*(mean(x)-x(1))/(n^0.2);
ierr=2;
end

function [fct]=funct(t,x)
p2=1.41421356;
n=length(x);
yy=zeros(size(x));
for i=1:n,
   xij=(x-x(i)).^2/t^2;
   y=exp(-xij/4).*((xij/2-1)/p2)-2*exp(-xij/2).*(xij-1);
   yy(i)=sum(y);
end;
fct=sum(yy)-2*n;
clear xij y yy;
end


function [ambd]=scaling(x,h)

% EVALUATES A VECTOR OF SCALING FACTORS FOR THE NONPARAMETRIC ADAPTATIVE
% ESTIMATION

% x - the n dimensional column vector of data values sorted in the ascending
% order
% h - the optimal smoothing factor
% ambd - the resultant n dimensional row vector of local scaling factors 

n=length(x);
c=sqrt(2*pi);
gau=zeros(1,n);
for i=1:n,
   gau(i)=sum(exp(-0.5*((x(i)-x)/h).^2))/c/n/h;
end
g=exp(mean(log(gau)));
ambd=sqrt(g./gau);
end

function [eps]=magn_accur(M)
x=sort(M);
d=x(2:length(x))-x(1:length(x)-1);
eps=min(d(d>0));
if eps>0.1; eps=0.1;end
end