% FUNCTION:         MM_MB
% VERSION:          [Wrapper Standalone Version]  V2.8
% COMPATIBLE with Matlab version 2017b or later
% REQUIREMENTS: "Statistics and Machine Learning" Matlab Toolbox installed
% TOOLBOX: "Magnitude Complexity Toolbox" within SERA Project
% DOCUMENT: "READ_ME_App_2B_v2_Description_MM_MB.docx"
% --------------------------------------------------------------------------------------------------------
% Test performed for the Distribution of a given dataset (time-series)
% checking whether PDF demonstrates multi-modes/ multi-bumps 
% --------------------------------------------------------------------------------------------------------
% OVERVIEW: This Application is a Matlab function which performs testing 
% of hypotheses of 1) multimodality and 2) existence of multi-bumps in a given 
% magnitude distribution. Please check also the accompanying auxiliary scripts 
% 'MM_MB_wrapper' and 'MM_MB_Plot' for a specific application (scenario)
% and plotting results.
% --------------------------------------------------------------------------------------------------------
% AUTHORS: K. Leptokaropoulos, S. Lasocki and P. Urban, 
% based on functions originally developed by S. Lasocki (Lasocki and Papadimitriou, 2006)
% last updated: 03/2019, within SERA PROJECT, EU Horizon 2020 R&I
% programme under grant agreement No.730900
% CURRENT VERSION: v2.8        ****   [Wrapper Standalone Version]
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% PLEASE refer to the accompanying document:
%                                                             "READ_ME_App_2B_v2_Description_MM_MB.docx" 
% for description of the Application and its requirements.
% --------------------------------------------------------------------------------------------------------
% DESCRIPTION:
% The function studies the magniutude distribution complexity by means of the 
% Mulitmodality Test (Silverman, 1986; Efron and Tibshirani, 1993). Two null 
% hypotheses (H0s) are tested:
% H01 - multimodality: The magnitude PDF is unimodal
% H02 - multi-bump: The magnitude PDF has one bump to the right of the 
% mode. 
% A mode is a local maximum of probability density and a bump is an interval 
% [a,b] such that the probability density is concave over [a,b] but not over any 
% larger interval (Silverman, 1986). The importance of modes and bumps relies 
% on the fact that multiple occurrences of these features  in a PDF indicate, for 
% most standard densities, a mixing of components (e.g. Cox, 1966).
% --------------------------------------------------------------------------------------------------------
% INPUT:  The function takes as input any matlab vector (input parameter "M").
%           The input data can be uploaded by the use of "ADTestMag_wrapper"
%           script, from an ASCII file (e.g. *.txt). Such file should contain a vector
%           (raw or column) of the Data that the User wishes to process. The User
%           is afterwards requested to enter values for some additional parameters.
%                Input Parameters Overview:
%                --- M: Time-Series (e.g. Magnitude) vector, read i.e. from an ASCII 
%                     file.
%               --- Mc: Corresponds to the input data completeness threshold.
%               --- m: number of points to divide the data sample, default 100, 
%                        recommended 100-1000. 
%      % INPUTS for MODE and BUMP HUNTING 
%               - n_boot: Number of bootstrap iterations for both MM and MB 
%                               (default: 100)
%               - delta_h: Smoothing parameter step for successive trials in defining 
%                                the critical h for MM and MB process (h-critical accuracy)
%                               (default: 0.001) 
%               - h : initial value of the smoothing factor to apply in defining the 
%                       critical h for the MM testing process (default:0.01)
%               - MMmeth: Method for multimodality Testing. Possible values:
%                         'Efron',   'Silverman'
%               - MBmeth: Method for multibumps Testing. Possible values:
%                         'Silverman',   'Efron'
%                NOTE: The initial value for MB process is taken equal to  the 
%                               critical value defined in the MM process minus its
%                               accuracy (i.e. hcrit_Modes-delta_h, see OUTPUT below)
% ---------------------------------------------------------------------------------------------------------
% OUTPUT:
%                 - Output Report with data and parameters used 
%                                      ('Output_MM_MB.txt file')
%                 - Output Parameters:
%         % DOUBLES:
%                        * n - Number of observations used
%                        * bval- b-value of the G-R law
%                        * Rmodes - The estimated significance of null hypothesis (H01)
%                        * hcrit_modes - estimated critical smoothness parameter for MM test
%                        * Rbumps - The estimated significance of null hypothesis (H02)
%                        * hcrit_bumps - estimated critical smoothness parameter for MB test
%                        * zer1 - point where 1st derivative is zero (extremum)
%                        * zer2 - point(s) where 2nd derivative is zero (inflection points)
%                                     [NOTE: zer2 can be a vector]
%        % VECTORS: 
%                        * gau: Data vector PDF estimated for h critical from Multimodality 
%                        * gau_b: Data vector PDF estimated for h critical from Multibumps 
%                        * poch: 1st derivative of Data vector PDF for h critical from Multibumps 
%                        * poch2: 2nd derivative of Data vector PDF for h critical from Multibumps
%                        * x1,x2: minimum and maximum data points after randomization 
%                                      (parameters needed for plotting)
%                 - Output Figures: Please refer to the Auxiliary Function:      "MM_MB_Plot"
%                                             
% ---------------------------------------------------------------------------------------------------------
% REFERENCES:
% -- Cox, D. R., (1966), Notes on the analysis of mixed frequency distributions. 
%     Br. J. Math. Stat. Psychol. 19, 39-47, doi.org/10.1111/j.2044-8317.1966.tb00353.x.
% -- Efron, B., and Tibshirani R.J. (1993) An Introduction to the Bootstrap,
%     CRC Press, Boca Raton, Fla. 
% -- Lasocki S. and E. E,  Papadimitriou (2006), "Magnitude distribution 
%    complexity revealed in seismicity from Greece", J. Geophys. Res., 
% 111, B11309, doi:10.1029/2005JB003794.
% -- Silverman, B. W., (1986), Density estimation for statistics and data analysis, 
%     CRC press, Boa Raton, Fla. 
% ---------------------------------------------------------------------------------------------------------
% LICENSE 
%     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 [n,bval,Rmodes,hcrit_modes,Rbumps,hcrit_bumps,gau,gau_b,poch,poch2,zer1,zer2,x1,x2]...
    =MM_MB_V2_8(M,Mc,m,n_boot,h,delta_h,MMmeth,MBmeth)
tic 
mkdir Outputs_MM_MB
% define round-off interval and select data above Mc 
sm=sort(M);dm=sm(2:length(sm))-sm(1:length(sm)-1);
dm=dm(dm>0);eps=min(dm);
M=M(M>=Mc);

% calculate beta
Rp=round(-log10(eps));                                   %% check the RP parameter
 [beta]=beta_AK_UT_Mbin (Mc,M,Rp);    

 % Randomize magnitude
 m_corr(:,1) = korekta(M,Mc,eps,beta);

 % SETTING PARAMETERS
n=length(m_corr); 
x1=min(m_corr)-0.001;           %a bit lower than Mmin
x2=max(m_corr)+0.001;         %a bit higher than Mmax



% MULTIMODALITY TESTING
[hcrit_modes,gau]=critical_smoothing(m_corr,n,delta_h,x1,x2,m,h);

% Run Multimodality testing
switch MMmeth
    case 'Efron' 
    [Rmodes] = test_multimodality_e (m_corr,n,x1,x2,m,hcrit_modes,n_boot);   
    case 'Silverman' 
    [Rmodes] = test_multimodality (m_corr,n,x1,x2,m,hcrit_modes,n_boot);   
end

% ------------------------------------------------------------------------
% [zer1]=zera_1st(m_corr,n,x1,x2,m,hcrit-delta_h);  
% ------------------------------------------------------------------------

% BUMP HUNTING
[hcrit_bumps,gau_b,poch,poch2]=critical_smoothing_bumps(m_corr,n,delta_h,x1,x2,m,hcrit_modes-delta_h);

% Run Bump hunting
switch MBmeth
    case 'Silverman' 
    [Rbumps] = bump_hunt (m_corr,n,x1,x2,m,hcrit_bumps,n_boot);   
    case 'Efron'  
    [Rbumps] = bump_hunt_e (m_corr,n,x1,x2,m,hcrit_bumps,n_boot);   
end

bval=beta/log(10);
% Display and Save the Results:
 re=table(n,bval,Rmodes,hcrit_modes,Rbumps,hcrit_bumps)

[zer1]=zera_1st(m_corr,n,x1,x2,m,hcrit_bumps-delta_h);  % use the same h in both cases
[zer2]=zera_2nd(m_corr,n,x1,x2,m,hcrit_bumps-delta_h);


toc

% Save Outputs
cd Outputs_MM_MB\
SaveOuts(eps,Mc,m,n_boot,n_boot,n,bval,Rmodes,hcrit_modes,Rbumps,hcrit_bumps)
cd ../
end

% -----  plotting magnitude density and its derivatives ------
% ****************************************************

%% ******************************************************************
% *************************** FUNCTIONS ***************************
% ****-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-****
%% --------------------------------------------------------------------------------------
%% --------------------------------------------------------------------------------------
%
function [beta]=beta_AK_UT_Mbin(Mmin,m,Rp)
%
%   m - magnitude vector
%   Mmin - completeness magitude threshold
%   beta - beta value. b(G-R)=beta/log(10)
%   Rp - Rounding precision, (1 - one decimal, 2 - two decimals, etc)
beta=1/(mean(m)-(Mmin-0.5*10^(-Rp)));
end

%% --------------------------------------------------------------------------------------
% Magnitude randomization
%
function [m_corr]=korekta(m,Mmin,eps,beta)
%
%   m - magnitude vector
%   Mmin - completeness magitude threshold
%   beta - beta value. b(G-R)=beta/log(10)
% EPS - magnitude round-off interval
%
% m_corr - randomized magnitude vector
%
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
%% --------------------------------------------------------------------------------------
% Finding critical smoothing parameter value for which the non parametric 
% PDF demonstrates one mode in the selected interval
%
function [hcrit,gau]=critical_smoothing(mm,n,delta_h,x1,x2,m,h)
%
x=sort(mm);
p=2.0;
c=sqrt(2*pi);
%x1=input('Lower limit ');
%x2=input('Upper limit ');
%m=input('No. of points ');
y=linspace(x1,x2,m);
%h=input('Initial smoothing factor ');
while p>1 
   for i=1:m,
        poch(i)=-sum((y(i)-x).*exp(-0.5*((y(i)-x)./h).^2));
   end
   p=0;
   for i=2:m,
       if poch(i-1)*poch(i)<0
           p=p+1;
       end
   end
   h=h+delta_h;
end
hcrit=h-delta_h;
for i=1:m,
   gau(i)=sum(exp(-0.5*((y(i)-x)./hcrit).^2))/c/n/hcrit;
end
end
%% --------------------------------------------------------------------------------------
% "SMOOTHED BOOTSTRAP"
%
function [R]=test_multimodality(m,n,x1,x2,n_point,h_crit,mm)
%
x=sort(m);
c=sqrt(2*pi);
%x1=input('Lower limit ');
%x2=input('Upper limit ');
%m=input('No. of points ');
y=linspace(x1,x2,n_point);
%h_crit=input('Critical smoothing factor ');
%mm=input('No of trials ');
no_wielom=0;
for j=1:mm,
	ind=unidrnd(n,n,1);
	eps=normrnd(0,1,n,1);
	z=x(ind)+h_crit.*eps;
	z=sort(z);
    poch(1)=-sum((y(1)-z).*exp(-0.5*((y(1)-z)./h_crit).^2));
    p=0;
    for i=2:n_point,
        poch(i)=-sum((y(i)-z).*exp(-0.5*((y(i)-z)./h_crit).^2));
        if poch(i-1)*poch(i)<0
           p=p+1;
        end
    end
	if p>1
        no_wielom=no_wielom+1;
	end
end
R=no_wielom/mm;
end
%% --------------------------------------------------------------------------------------
% "SMOOTHED BOOTSTRAP"
% EFRON
%
function [R]=test_multimodality_e(m,n,x1,x2,n_point,h_crit,mm)
%
x=sort(m);
s2=sqrt(1+h_crit^2/var(x));
c=sqrt(2*pi);
%x1=input('Lower limit ');
%x2=input('Upper limit ');
%m=input('No. of points ');
y=linspace(x1,x2,n_point);
%h_crit=input('Critical smoothing factor ');
%mm=input('No of trials ');
no_wielom=0;
for j=1:mm,
	ind=unidrnd(n,n,1);
	eps=normrnd(0,1,n,1);
    me=mean(x(ind));
	z=me+(x(ind)-me+h_crit.*eps)./s2;
	z=sort(z);
    poch(1)=-sum((y(1)-z).*exp(-0.5*((y(1)-z)./h_crit).^2));
    p=0;
    for i=2:n_point,
        poch(i)=-sum((y(i)-z).*exp(-0.5*((y(i)-z)./h_crit).^2));
        if poch(i-1)*poch(i)<0
           p=p+1;
        end
    end
	if p>1
        no_wielom=no_wielom+1;
	end
end
R=no_wielom/mm;
end
%% --------------------------------------------------------------------------------------
% Finding critical smoothing parameter value for which the non parametric 
% PDF demonstrates one bump in the selected interval
%
function [hcrit,gau,poch,poch2]=critical_smoothing_bumps(mm,n,delta_h,x1,x2,m,h)
%
x=sort(mm);
p=2.0;
c=sqrt(2*pi);
%x1=input('Lower limit ');
%x2=input('Upper limit ');
%m=input('No. of points ');
y=linspace(x1,x2,m);
%h=input('Initial smoothing factor ');
while p>1 
    for i=1:m,
%        poch(i)=-sum((y(i)-x).*exp(-0.5*((y(i)-x)./h).^2));
        poch2(i)=sum((((y(i)-x)./h).^2-1).*exp(-0.5*((y(i)-x)./h).^2));
    end
   p=0;
for i=2:m,
   if poch2(i-1)*poch2(i)<0
      p=p+1;
   end
end
h=h+delta_h;
end
hcrit=h-delta_h;
for i=1:m,
    gau(i)=sum(exp(-0.5*((y(i)-x)/hcrit).^2))/c/n/hcrit;
    poch(i)=-sum((y(i)-x).*exp(-0.5*((y(i)-x)./hcrit).^2));
    poch2(i)=sum((((y(i)-x)./hcrit).^2-1).*exp(-0.5*((y(i)-x)./hcrit).^2));
end
end
%% --------------------------------------------------------------------------------------
% "SMOOTHED BOOTSTRAP"
%
function [R]=bump_hunt(m,n,x1,x2,n_point,h_crit,mm)
%
x=sort(m);
c=sqrt(2*pi);
%x1=input('Lower limit ');
%x2=input('Upper limit ');
%m=input('No. of points ');
y=linspace(x1,x2,n_point);
%h_crit=input('Critical smoothing factor ');
%mm=input('No of trials ');
no_bump=0;
for j=1:mm,
	ind=unidrnd(n,n,1);
	eps=normrnd(0,1,n,1);
	z=x(ind)+h_crit.*eps;
	z=sort(z);
%    poch(1)=-sum((y(1)-z).*exp(-0.5*((y(1)-z)./h_crit).^2));
    poch2(1)=sum((((y(1)-z)./h_crit).^2-1).*exp(-0.5*((y(1)-z)./h_crit).^2));
    p=0;
    for i=2:n_point,
 %       poch(i)=-sum((y(i)-z).*exp(-0.5*((y(i)-z)./h_crit).^2));
        poch2(i)=sum((((y(i)-z)./h_crit).^2-1).*exp(-0.5*((y(i)-z)./h_crit).^2));
        if poch2(i-1)*poch2(i)<0
           p=p+1;
        end
    end
	if p>1
        no_bump=no_bump+1;
	end
end
R=no_bump/mm;
end
%% --------------------------------------------------------------------------------------
% "SMOOTHED BOOTSTRAP"
% EFRON
%
function [R]=bump_hunt_e(m,n,x1,x2,n_point,h_crit,mm)
%
x=sort(m);
s2=sqrt(1+h_crit^2/var(x));
c=sqrt(2*pi);
%x1=input('Lower limit ');
%x2=input('Upper limit ');
%m=input('No. of points ');
y=linspace(x1,x2,n_point);
%h_crit=input('Critical smoothing factor ');
%mm=input('No of trials ');
no_bump=0;
for j=1:mm,
	ind=unidrnd(n,n,1);
	eps=normrnd(0,1,n,1);
    me=mean(x(ind));
	z=me+(x(ind)-me+h_crit.*eps)./s2;
	z=sort(z);
%    poch(1)=-sum((y(1)-z).*exp(-0.5*((y(1)-z)./h_crit).^2));
    poch2(1)=sum((((y(1)-z)./h_crit).^2-1).*exp(-0.5*((y(1)-z)./h_crit).^2));
    p=0;
    for i=2:n_point,
 %       poch(i)=-sum((y(i)-z).*exp(-0.5*((y(i)-z)./h_crit).^2));
        poch2(i)=sum((((y(i)-z)./h_crit).^2-1).*exp(-0.5*((y(i)-z)./h_crit).^2));
        if poch2(i-1)*poch2(i)<0
           p=p+1;
        end
    end
	if p>1
        no_bump=no_bump+1;
	end
end
R=no_bump/mm;
end

%% --------------------------------------------------------------------------------------
% Determination of PDF first derivative zeros estimates in the selected interval
%
function [zer1]=zera_1st(mm,n,x1,x2,m,h)
%
x=sort(mm);
c=sqrt(2*pi);
y=linspace(x1,x2,m);
    for i=1:m,
         poch(i)=-sum((y(i)-x).*exp(-0.5*((y(i)-x)./h).^2));
     end
     p=0;
     for i=2:m,
         if poch(i-1)*poch(i)<0
             p=p+1;
             zer1(p)=(y(i)+y(i-1))/2;                     %K     23JAN2019
         end
     end

end
%% --------------------------------------------------------------------------------------
% Determination of PDF second derivative zeros estimates in the selected interval
%

function [zer2]=zera_2nd(mm,n,x1,x2,m,h)
%
x=sort(mm);
c=sqrt(2*pi);
y=linspace(x1,x2,m);
    for i=1:m,
         poch2(i)=sum((((y(i)-x)./h).^2-1).*exp(-0.5*((y(i)-x)./h).^2));
     end
     p1=0;
     for i=2:m,
         if poch2(i-1)*poch2(i)<0
             p1=p1+1;
             zer2(p1)=(y(i)+y(i-1))/2;                     %K     23JAN2019
         end
     end
  
end
%% --------------------------------------------------------------------------------------------------------
% --------------------------------------- SAVE OUTPUTS in the report file ---------------------------------------
% Save Outputs
function SaveOuts(EPS,Mmin,m,n_boot1,n_boot2,n,bval,Rmodes,hcrit_modes,Rbumps,hcrit_bumps)
% ---- Save *.txt file with Parameters Report ----
%cd Outputs/
fid=fopen('REPORT_Multimodality.txt','w');
fprintf(fid,['RESULTS from MULTIMODALITY/MULTIBUMP TESTING (created on ', datestr(now),')\n']);
fprintf(fid,'------------------------------------------------------------------------------\n');
fprintf(fid,['<Data Round-off interval                              >: ', num2str(EPS),'\n']);
fprintf(fid,['<Number of points to divide the sample                >: ', num2str(m),'\n']);
fprintf(fid,['<Completeness Threshold                               >: ', num2str(Mmin),'\n']);
fprintf(fid,['<Number of events used                                >: ', num2str(n),'\n']);
fprintf(fid,['<Gutenberg-Richter b-value                            >: ', num2str(bval,'%5.3f'),'\n']);
fprintf(fid,['<Number of bootstrap iterations (multimodality)       >: ', num2str(n_boot1),'\n']);
fprintf(fid,['<Number of bootstrap iterations (multibumps)          >: ', num2str(n_boot2),'\n']);
fprintf(fid,'------------------------------------------------------------------------------\n');
fprintf(fid,['<Critical Smoothing Parameter,h (for modes)           >: ', num2str(hcrit_modes,'%6.4f'),'\n']);
fprintf(fid,['<Critical Smoothing Parameter,h (for bumps)           >: ', num2str(hcrit_bumps,'%6.4f'),'\n']);
fprintf(fid,['<p-value of H01 that input data PDF is unimodal       >: ', num2str(Rmodes,'%4.3f'),'\n']);
fprintf(fid,['<p-value of H02 that input data PDF has no more than \n']);
fprintf(fid,['         one bump to the right of the global maximum  >: ', num2str(Rbumps,'%4.3f'),'\n']);
  


fclose(fid);

end