%
%[lamb_all,lamb,lamb_err,unit,eps,b,Mmax,err,BIAS,SD]=TruncGR_Ob(t,M,iop,Mmin,Mmax,Nsynth)
%
% ESTIMATES THE MEAN ACTIVITY RATE WITHIN THE WHOLE SAMPLE AND WITHIN THE 
% COMPLETE PART OF THE SAMPLE, THE ROUND-OFF ERROR OF MAGNITUDE, 
% THE GUTENBERG-RICHTER B-VALUE AND THE UPPER BOUND OF MAGNITUDE
% DISTRIBUTION USING THE UPPER-BOUNDED G-R LED MAGNITUDE DISTRIBUTION MODEL
%
% !! THIS FUNCTION MUST BE EXECUTED AT START-UP OF THE UPPER-BOUNDED 
%       GUTENBERG-RICHETR HAZARD ESTIMATION MODE !!
%
% AUTHOR: S. Lasocki ver 2 01/2015 within IS-EPOS project.
%
% DESCRIPTION: The assumption on the upper-bounded Gutenberg-Richter 
% relation leads to the upper truncated exponential distribution to model 
% magnitude distribution from and above the catalog completness level 
% Mmin. The shape parameter of this distribution and consequently the G-R
% b-value is estimated by maximum likelihood method (Aki-Utsu procedure).   
% The upper limit of the distribution Mmax is evaluated using 
% the Kijko-Sellevol generic formula. If convergence is not reached the
% Whitlock @ Robson simplified formula is used: 
% Mmaxest= 2(max obs M) - (second max obs M)).
% The mean activity rate, lamb, is the number of events >=Mmin into the
% length of the period in which they occurred. Upon the value of the input 
% parameter, iop, the used unit of time can be either day ot month or year.    
% The round-off interval length - eps is the least non-zero difference 
% between sample data or 0.1 if the least difference is greater than 0.1.
%
% REFERENCES:
%Kijko, A., and M.A. Sellevoll (1989) Bull. Seismol. Soc. Am. 79, 3,645-654 
%Lasocki, S., Urban, P. (2011) Acta Geophys 59, 659-673, 
% doi: 10.2478/s11600-010-0049-y
%
% INPUT:
%   t - vector of earthquake occurrence times
%   M - vector of magnitudes from a user selected catalog
%   iop - determines the used unit of time. iop=0 - 'day', iop=1 - 'month', 
%       iop=2 - 'year'
%   Mmin - catalog completeness level. Must be determined externally.
%   Can take any value from [min(M), max(M)].
%   Mmax - upper limit of Magnitude Distribution. Can be set by user, or
%   estimate within the program - it then should be set as Mmax=[].
%
% 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 15 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 the round-off interval of magnitudes.
%   b - Gutenberg-Richter b-value
%   Mmax - upper limit of magnitude distribution
%   err - error parameter on Mmax estimation, err=0 - convergence, err=1 - 
%   no converegence of Kijko-Sellevol estimator, Robinson @ Whitlock
%   method used.
%   BIAS - Mmax estimation Bias (Lasocki and Urban, 2011)
%   SD   - Mmax standard deviation (Lasocki ands Urban, 2011)
%
% 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,b,Mmax,err,BIAS,SD]=TruncGR_Ob(t,M,iop,Mmin,Mmax,Nsynth)
if nargin==5;Nsynth=[];end                                   % K08DEC2019
if isempty(t) || numel(t)<3 || isempty(M(M>=Mmin)) %K03OCT
    t=[1 2];M=[1 2];   end   %K30SEP

n=length(M);
lamb_err=0;
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

% SL 03MAR2015 ----------------------------------
[NM,unit]=time_diff(t(1),t(n),iop);    
lamb_all=n/NM;
[NM,unit]=time_diff(t1,t2,iop);        
lamb=nn/NM;
% SL 03MAR2015 ----------------------------------

if nn<15
    eps=0;b=0;Mmax=0;err=0;
    lamb_err=1;
    BIAS=NaN;SD=NaN;                                               %%%     K 08NOV2019  
    return;
end

eps=magn_accur(M);
xx=M(M>=Mmin); %K21OCT2014
% x=sort(M,'descend');
% for i=1:n
%      if x(i)<Mmin; break; end
%      xx(i)=x(i);                                                      %
% end

clear x;
nn=length(xx);
    
Max_obs=max(xx);
beta0=0;
Mmax1=Max_obs; 
            if isempty(Mmax)==0                                      %%%   K 28JUL2015
                	err_Mmax=0;                                      %%%   K 04DEC2019 
                    fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax);  %%%   K 28JUL2015
                	x0 = 1;              %[0.05,4.0];                %%%   K 28JUL2015 - See exception line 155
                    beta = fzero(fun,x0);                            %%%   K 28JUL2015
                    err=0;                                           %%%   K 28JUL2015
            else                                                     %%%   K 28JUL2015  -  line 150
                    err_Mmax=1;                                      %%%   K 04DEC2019 
for i=1:50,                                                           
	fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax1);
	x0 =1;           %[0.05,4.0];                               %%%   K29JUL2015 - See exception line 155
    beta = fzero(fun,x0);
  	Mmax=Max_obs+moja_calka('f_podc',Mmin,Max_obs,1e-5,nn,beta,Mmin-eps/2,Mmax1);
    if ((abs(Mmax-Mmax1)<0.01)&&(abs(beta-beta0)<0.0001))
        err=0;
        break;
    end
    Mmax1=Mmax;
    beta0=beta;
end
if i==50; 
    err=1.0;
    Mmax=2*xx(1)-xx(2);
fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax);
	x0 =1;
    beta = fzero(fun,x0);
end
            end                                                         %%%   K 28JUL2015
b=beta/log(10);
clear xx

% Exception for v-value
if b<0.05 || b>6.0; error('Unacceptable b-value, abort and select different dataset');end
beta;


% Estimation of Mmax Bias                                               %%%   K 04DEC2019 
% (Lasocki and Urban, 2011, doi:10.2478/s11600-010-0049-y)
if isempty(Nsynth)==0 && err_Mmax==1              % set number of synthetic datasets, e.g. 10000
[BIAS,SD]=Mmax_Bias_GR(t,M,Mmin,Mmax,b,err,Nsynth);
elseif isempty(Nsynth)==0 && err_Mmax==0;              
warning('Mmax must be empty for BIAS calculation');BIAS=[];SD=[];
else; BIAS=0;SD=0;
end

end




function [NM,unit]=time_diff(t1,t2,iop)               % SL 03MAR2015 

% TIME DIFFERENCE BETWEEEN t1,t2 EXPRESSED IN DAY, MONTH OR YEAR UNIT
%
% t1 - start time (in MATLAB numerical format)
% t2 - end time (in MATLAB numerical format)  t2>=t1
%   iop - determines the used unit of time. iop=0 - 'day', iop=1 - 'month', 
%       iop=2 - 'year'
%
%   NM - number of time units from t1 to t2
%   unit - string with name of time unit used ('year' or 'month' or 'day').

if iop==0
    NM=(t2-t1);
    unit='day';
elseif iop==1
    V1=datevec(t1);
    V2=datevec(t2);
    NM=V2(3)/eomday(V2(1),V2(2))+V2(2)+12-V1(2)-V1(3)/eomday(V1(1),V1(2))...
        +(V2(1)-V1(1)-1)*12;
    unit='month';
else
    V1=datevec(t1);
    V2=datevec(t2);
    NM2=V2(3);
    if V2(2)>1
        for k=1:V2(2)-1
            NM2=NM2+eomday(V2(1),k);
        end
    end
    day2=365; if eomday(V2(1),2)==29; day2=366; end;
    NM2=NM2/day2;
    NM1=V1(3);
    if V1(2)>1
        for k=1:V1(2)-1
            NM1=NM1+eomday(V1(1),k);
        end
    end
    day1=365; if eomday(V1(1),2)==29; day1=366; end;
    NM1=(day1-NM1)/day1;
    NM=NM2+NM1+V2(1)-V1(1)-1;
    unit='year';
end
end

function [zero]=bet_est(b,ms,Mmin,Mmax)

%First derivative of the log likelihood function of the upper-bounded 
% exponential distribution (truncated GR model)
% b - parameter of the distribution 'beta'
% ms - mean of the observed magnitudes
% Mmin - catalog completeness level
% Mmax - upper limit of the distribution

M_max_min=Mmax-Mmin;
e_m=exp(-b*M_max_min);
zero=1/b-ms+Mmin-M_max_min*e_m/(1-e_m);
end


function [calka,ier]=moja_calka(funfc,a,b,eps,varargin)

% Integration by means of 16th poit Gauss method. Adopted from CERNLIBRARY

% funfc - string with the name of function to be integrated
% a,b - integration limits
% eps - accurracy
% varargin - other parameters of function to be integrated
% calka - integral
% ier=0 - convergence, ier=1 - no conbergence

persistent W X CONST
W=[0.101228536290376 0.222381034453374 0.313706645877887 ...
0.362683783378362 0.027152459411754 0.062253523938648 ...
0.095158511682493 0.124628971255534 0.149595988816577 ...
0.169156519395003 0.182603415044924 0.189450610455069];
X=[0.960289856497536 0.796666477413627 0.525532409916329 ...
0.183434642495650 0.989400934991650 0.944575023073233 ...
0.865631202387832 0.755404408355003 0.617876244402644 ...
0.458016777657227 0.281603550779259 0.095012509837637];
CONST=1E-12;
delta=CONST*abs(a-b);
calka=0.;
aa=a;
y=b-aa;
ier=0;
while abs(y)>delta,
   bb=aa+y;
   c1=0.5*(aa+bb);
   c2=c1-aa;
   s8=0.;
   s16=0.;
   for i=1:4,
      u=X(i)*c2;
      s8=s8+W(i)*(feval(funfc,c1+u,varargin{:})+feval(funfc,c1-u,varargin{:}));
   end
   for i=5:12,
      u=X(i)*c2;
      s16=s16+W(i)*(feval(funfc,c1+u,varargin{:})+feval(funfc,c1-u,varargin{:}));
   end
   s8=s8*c2;
   s16=s16*c2;
   if abs(s16-s8)>eps*(1+abs(s16))
      y=0.5*y;
      calka=0.;
      ier=1;
   else
      calka=calka+s16;
      aa=bb;
      y=b-aa;
      ier=0;
   end
end
end

function [y]=f_podc(z,n,beta,Mmin,Mmax)

% Integrated function for Mmax estimation. Truncated GR model
% z - column vector of independent variable
% n - the size of 'z'
% beta - the distribution parameter
% Mmin - the catalog completeness level
% Mmax - the upper limit of the distribution

y=Cdfgr(z,beta,Mmin,Mmax).^n;
end

function [y]=Cdfgr(t,beta,Mmin,Mmax)

%CDF of the truncated upper-bounded exponential distribution (truncated G-R
% model
% Mmin - catalog completeness level
% Mmax - upper limit of the distribution
% beta - the distribution parameter
% t - vector of magnitudes (independent variable)
% y - CDF vector

mian=(1-exp(-beta*(Mmax-Mmin)));
y=(1-exp(-beta*(t-Mmin)))/mian;
idx=find(y>1);
y(idx)=ones(size(idx));
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

% --------------------- Mmax BIAS estimation routine ----------------------  K 08NOV2019

function [BIAS,SD]=Mmax_Bias_GR(t,m,Mc,Mmax1,b,err,synth)

if err~=0
warning('process did not converge!!'),pause   
end

MAXm=max(m);beta=b*log(10);N=numel(m(m>=Mc));DeltaM=MAXm-Mc;

for j=1:synth         %set number of synthetic datasets, default is 10000 
% % CDF:
 M=Mc:0.0001:MAXm;upt=1-exp(-beta*(M-Mc));
 dwt=1-exp(-beta*(MAXm-Mc));F=upt./dwt;                                             % j
% linear interpolation to assign magnitude values from a uniform distribution sample
 iM=rand(1,N);M1=interp1q(F',M',iM');
 br(j)=1/(log(10)*(mean(M1)-min(M1)));DM=range(M1);
 Mmax=max(M1);

 % Iteration Process to estimate b and Mmax
b1=1;best=[1.0 10.0];i=1; 
while min(abs(diff(best)))>0.00001
w=exp(b1*(Mmax-Mc));E1=expint(N/(w-1));E2=expint(N*w/(w-1));
%E=expint(w);
Mme=Mmax+(E1-E2)/(b1*exp(-N/(w-1)))+(Mc)*exp(-N);    %Mme=round(Mme/EPS)*EPS;            

if isnan(Mme)
KM=sort(unique(M1),'descend');
Mme=2*KM(1)-KM(2);
end

fun=@(bb) 1/bb+(Mme-Mc)/(1-exp(bb*(Mme-Mc)))-mean(M1)+Mc;  %consider th5 last 0.05 term
b1=fzero(fun,1);best(i)=b1;i=i+1;

if i==50
   warning('process did not converge!!');break 
end

end

be(j)=b1/log(10);
Me(j)=Mme;dm(j)=DM;Mm(j)=Mmax;

end

BIAS=mean(MAXm-Me)
SD=std(MAXm-Me);

%b-mean(be)    %check b-value difference
%histogram(be)

% MAXm: maximum magnitude in the real catalog 
% Mmax: maximum magnitudes observed in the synthetic catalogs (rounded)
% Me: maximum magnitude estimates for the synthetic catalogs 
% Mmax1: maximum magnitude estimated by GRT
end



