% PROGRAM     :  SHAPE [Seismic HAzard Parameters Evaluation] 
% VERSION     :  V_1.2 [Interactive Standalone Version] 
% LAST UPDATED:  April 2020 
% COMPATIBLE  :  Matlab version 2017b or later
% REQUIREMENTS:  "Statistics and Machine Learning" and "Image Processing" Matlab Toolbox installed
% TOOLBOX     :  "Hazard Analysis Toolbox" developed within SERA Project
% USER GUIDE  :  "READ_ME_SHAPE_ver1.pdf"
% CITE AS     :  Leptokaropoulos K. and S. Lasocki (2020), Seismol. Res. Lett., doi: 10.1785/0220190319
% CONTACT     :  kleptoka@igf.edu.pl (Dr. Kostas Leptokaropoulos)
% --------------------------------------------------------------------------------------------------------------------
% Time-and-Technology Dependent Seismic Hazard Assessment (SHA)  
% --------------------------------------------------------------------------------------------------------------------
% INPUT:
% !!! ----------------------------  INPUT DATA REQUIREMENTS ----------------------------- !!!
% the program works with ASCII input data files (e.g. *.txt). The files needed are:
% >   File with the parameters of seismic data [mandatory]
% >   File with the parameters of production data [optional]
% >   File specifying time windows for SHA analysis [optional]
% >   File with the fields description of the corresponding parameters in the seismic data file
% >   File with the fields description of the corresponding parameters in the production data files
% FOR DETAILS on data requirements please refer to the document:
%                            "READ_ME_SHAPE_ver1.pdf"
% --------------------------------------------------------------------------------------------------------------------
% OVERVIEW: THE PROGRAM takes as input a Seismic and optionally, a Production data parameter file to provide Seismic 
% Hazard Assessment for specified time-windows after filtering appropriately the data following User's specifications. 
% --------------------------------------------------------------------------------------------------------------------
% AUTHORS: K. Leptokaropoulos and S. Lasocki
% Based on magnitude distribution/stationary hazard parameters estimation functions originally developed by S. Lasocki
% Last Updated: 04/2020, within SERA PROJECT, EU Horizon 2020 R&I programme under grant agreement No.730900
% CURRENT VERSION: v1.2        ****    [INTERACTIVE STANDALONE VERSION!!]
%% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% PLEASE refer to the accompanying document:
%                                           "READ_ME_SHAPE_ver1.pdf" 
% for description of the Application and its requirements.
%% -----------------------------------------------------------------------------------------------------------------------
% DESCRIPTION: The Application performs time-dependent Seismic Hazard Analysis (SHA), 
% taking into account the activity rate and the magnitude distribution of seismicity 
% for selected time windows. The hazard parameters estimated are: 
%    1) The Mean Return Period (MRP) of a given magnitude, M, which is defined as the 
%       average elapsed time between the occurrence of consecutive events of M and 
%    2) The Exceedance Probability (EPR) of a given magnitude, M, within a given time 
%       period of length, T, which is defined as the probability of an earthquake of 
%       M to occur during T.  
% These hazard parameters are estimated for different time windows which are constructed 
% upon User’s particular specifications. 4 different magnitude distribution models can 
% be chosen. The input files must be in ASCII format (e.g. *.txt). A brief description 
% of the preparation process is given here:
% (please see in the script for more comments and details for each STEP)
% STEP_1. MODE SELECTION: Select '1' for Seismic Data, '2' for both Seismic & Production Data
% STEP_2. DATA SELECTION: Depending on the MODE selection (STEP 1), the user is requested 
%         here to select data and data field files
% STEP_3. MAGNITUDE SCALE SELECTION 
% STEP_4. DATA FILTERING: Before proceeding to Seismic Hazard Analysis the User is given 
%         the option to filter the data according to
%         - Time
%         - Epicentral Location
%         - Events Depths
%         - Magnitude
% STEP_5. TIME WINDOWS GENERATION: The data in here divided according to time windows 
%         defined by the User, by means of 4 different modes
%         - Time
%         - Events
%         - Graphical 
%         - Read from File 
% STEP_6. SHA PARAMETERS SELECTION: Selection of Magnitude Distribution Model, Time Unit, 
%         Magnitude(for exceedance probability and mean return period estimate) and Time 
%         Period (for exceedance probability estimate)
% STEP_7. OUTPUTS: Generate and save output files, and figures
% ---------------------------------------------------------------------------------------------------------------------
%% INPUT:   All input data are sufficiently explained in the script as well as 
%               while running the code (interaction with the user). NOTE that
%               all input files (seismic catalog, production data, time windows) 
%               must be in ASCII format (i.e. *.txt). 
%               Please refer to the APPLICATION DOCUMENTATION for further 
%               instructions and input data requirement specifications:  "READ_ME_SHAPE_ver1.pdf"
% ----------------------------------------------------------------------------------------------------------------------
%% OUTPUT:
%                 <> Output Report with summary of the Results as well as data and parameters used 
%                 <> Output Figure with the results in *.mat and *.jpeg formats
%                 <> Output Matlab Structure with input parameter values and output results, having
%                    as many cells as the number of time windows generated. 
%                                 Structure fields are:
%                    - Time      : vector with origin times of the events included in each time window
%                    - M         : vector with events magnitudes
%                    - Mmin      : Completeness magnitude
%                    - eps       : Magnitude round-off interval
%                    - lambd     : mean activity rate
%                    - lambd_err : events number sufficiency (0-all parametes estiamated, 1-all parameters set as NaNs)
%                    - unit      : Time Unit
%                    - method    : Magnitude Distribution Model
%                    - b         : b-value of GR law    
%                                  [applies only when "method" is set to 'GRU' or 'GRT']
%                    - h         : Kernel smoothing factor 
%                                  [applies only when "method" is set to 'NPU' or 'NPT']
%                    - xx        : Background sample for kernel magnitude estimate 
%                                  [applies only when "method" is set to 'NPU' or 'NPT']
%                    - ambd      : weigthing factors for the adaptive kernel 
%                                  [applies only when "method" is set to 'NPU' or 'NPT']
%                    - ierr      : h convergence indicator (0-converges,1-multiple zeros, 2-no zeros) 
%                                  [applies only when "method" is set to 'NPU' or 'NPT']
%                    - Mmax      : Upper limit of magnitude distribution (truncated) 
%                                  [applies only when "method" is set to 'GRT' or 'NPT']
%                    - err       : Mmax convergence indicator (0-converge, 1-no converge) 
%                                  [applies only when "method" is set to 'GRT' or 'NPT']
%                    - PDF       : 2 columns, first representing magnitudes and second the Probability Density Function of those magnitudes
%                    - PDF       : 2 columns, first representing magnitudes and second the Cumulative Distribution Function of those magnitudes
% -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- 
% REFERENCES:
% Kijko A, Lasocki S, Graham G (2001), Pure Appl. Geophys. 158:1655–1676
% Kijko A, Sellevoll MA (1989), Bull Seismol. Soc. Am. 79:645–654
% Lasocki S (2017), Chapter 11.3 in Rockburst Mechanisms, Oxford, pp 366-380
% Lasocki S, Urban P (2011), Acta Geophys. 59:659–673
% Lasocki S, Orlecka-Sikora B (2008), Tectonophysics 456:28–37
% Leptokaropoulos K, Staszek M, Cielesta S, Urban P, Olszewska D, Lizurek G (2017), Acta Geophys. 65:493-505 
% Leptokaropoulos K, Lasocki S (2020), Seismol. Res. Lett., doi: 10.1785/0220190319.
% ---------------------------------------------------------------------------------------------------------------------
% 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.
% -------------------------------------------------------------------------------------------------------------------


clear;clc
mkdir Outputs_SHA
set(0, 'DefaultUICOntrolFontSize', 12)

% Browse Directury with Functions

%% ***********    STEP_1 and STEP_2: MODE and DATA SELECTION    **************
[Catalog,PROD_Data,Mode1,dstr1,dstr2,s,s1,ss,s2]=Data_Hand_A2M;
%% ***********          STEP_3: SELECT MAGNITUDE SCALE          **************
[Ctime,Cmag,Mtype]=Select_Magnitude_Scale(Catalog);
Mc=min(Cmag);
%% ***********              STEP_4: DATA FILTERING              ************** 

% Ask for additional filtering
opts.Interpreter='tex';opts.Default='Yes';
quest='Do you wish to filter Data?';
answer=questdlg(quest,'Data Filtering','Yes','No',opts);
while strcmp(answer,'Yes')
cd Filtering
% FILTERING FOR MC, Depth, Time and Epicentral coordinates
[out]=Toolbox_SHA_TEST_2
    if out==4
    [Ctime,Cmag,Catalog,Mc]=FiltMc(Ctime,Cmag,Catalog,s1);
    elseif out==3
    [Ctime,Cmag,Catalog,Mc]=FiltDep(Ctime,Cmag,Catalog,s1);
    elseif out==1
    [Ctime,Cmag,Catalog,Mc]=FiltTime(Ctime,Cmag,Catalog,PROD_Data,s1);
    elseif out==2     
    [Ctime,Cmag,Catalog,Mc]=FiltSpace(Ctime,Cmag,Catalog,s1);
    end

quest='Do you wish to filter Data?';
answer=questdlg(quest,'Data Filtering','Yes','No',opts);
cd ../    
end
disp('Proceed to Parameters Setting')



%% % *************         STEP_5: CREATE TIME WINDOWS         ************** 

time_windows=struct;time_windows.Time=[];
 to=Ctime-Ctime(1);tmin=min(to);tmax=max(to);

lista = {'Time','Events','Graphical','File'};                              % select one of the windows creation mode
answer1 = listdlg('PromptString','Please Select a Window Creation Mode: ', ...
	         'SelectionMode','single','ListString',lista,'ListSize',[350 250],'InitialValue',3);
answer1=char(lista(answer1));         
switch answer1
    %% TIME
    case 'Time'
prompt = {'\fontsize{12} Enter window size (DAYS):','\fontsize{12} Enter window step (DAYS):'};
dlgtitle = 'Configure Time Windows';
dims = [1 50]; opts.Interpreter='tex';
definput = {'30','1'};
answer = inputdlg(prompt,dlgtitle,dims,definput,opts);        
window_size=str2double(answer(1));dt=str2double(answer(2));        
        if window_size>tmax;n=1;warning('time window is set larger than data time span');end
        n=ceil((tmax-window_size)/dt);
        for i=1:n
        time_windows(i).Time=Ctime(to>=(i-1)*dt & to<(i-1)*dt+window_size);
        time_windows(i).M=Cmag(to>=(i-1)*dt & to<(i-1)*dt+window_size);
        time_windows(i).Tstart=Ctime(1)+(i-1)*dt;
        time_windows(i).Tend=Ctime(1)+(i-1)*dt+window_size;
        end
     %% EVENTS
    case 'Events'
prompt = {'\fontsize{12} Enter window size (EVENTS):','\fontsize{12} Enter window step (DAYS):'};
dlgtitle = 'Configure Time Windows';
dims = [1 50]; opts.Interpreter='tex';
definput = {'100','1'};
answer = inputdlg(prompt,dlgtitle,dims,definput,opts);        
window_size=str2double(answer(1));dt=str2double(answer(2));

 if window_size>numel(to);window_size=numel(to);warning('events window is set larger than given data');end    
   n=ceil(tmax/dt);
        for i=1:n
        To=find(to>=(i-1)*dt);To=To(1);
        if To<=length(Ctime)-window_size+1;
        time_windows(i).Time=Ctime(To:To+window_size-1);
        time_windows(i).M=Cmag(To:To+window_size-1);
        time_windows(i).Tstart=Ctime(1)+(i-1)*dt;
        time_windows(i).Tend=max(time_windows(i).Time);
        end
        end
    %% GRAPHICAL        
    case 'Graphical'
fig1=figure('Position',[300 100 1600 500]);
plot(Ctime,1:numel(Ctime));xlabel(['Time'],'FontSize',16);ylabel(['Events'],'FontSize',16);hold on;
if isempty(PROD_Data)==0;
yyaxis right;plot(PROD_Data(1).val,PROD_Data(s2).val);ylabel(PROD_Data(s2).field,'FontSize',16,'interpreter','none');
end
datetick('x',20);
title(['Please select starting/ending points for time periods    (press enter to submit)'],'FontSize',16)
[x,y]=ginput;close(fig1)

if length(x)>1
time_start_vector=x(1:length(x)-1);time_end_vector=x(2:length(x));
else; time_start_vector=Ctime(1)-1;time_end_vector=Ctime(length(Ctime))+1;
end

if (length(time_start_vector) ~= length(time_end_vector)); error('time start and end vectors must have the same size'); end 
	time_windows=struct;
	time_windows.Time=[];
	time_windows.M=[];
 	
 	for i = 1:length(time_start_vector)
 		startTime = time_start_vector(i);
 		endTime = time_end_vector(i);
 		time_windows(i).Time=Ctime(Ctime > startTime & Ctime <= endTime);
        time_windows(i).M=Cmag(Ctime > startTime & Ctime <= endTime);
        time_windows(i).Tstart = startTime;
        time_windows(i).Tend = endTime;
    end
    %% Read from FILE
    case 'File'                                                     %K 27JUN2019
        cd TIME_WINDOWS                                             %K 27JUN2019
    d=dir;dstr={d.name};                                            %K 27JUN2019
    dlabela='Select Time Windows File:';                            %K 27JUN2019
    [s,ok]=listdlg('PromptString',dlabela,...                       %K 27JUN2019
    'SelectionMode','single','ListString',dstr,...                  %K 27JUN2019
    'ListSize',[350 250]);                                          %K 27JUN2019
    if ok; Twindows=dlmread(dstr{s});end                            %K 27JUN2019
                     
    T1=Twindows(:,1);T2=Twindows(:,2);n=numel(T1);                  %K 27JUN2019
    for i=1:n                                                       %K 27JUN2019
       time_windows(i).Time=Ctime(Ctime>=T1(i) & Ctime<T2(i));      %K 27JUN2019
       time_windows(i).M=Cmag(Ctime>=T1(i) & Ctime<T2(i));          %K 27JUN2019
       time_windows(i).Tstart=T1(i);                                %K 27JUN2019
       time_windows(i).Tend=T2(i);                                  %K 27JUN2019
    end                                                             %K 27JUN2019
    cd ../                                                          %K 27JUN2019
    
end

%% % *************      STEP_6: SELECT PARAMETERS for SHA      ************** 
% Magnitude Distribution Mode
list1 = {'Unbounded Gutenberg-Richter','Upper-bounded Gutenberg-Richter',...
    'Unbounded non-parametric Kernel','Upper-bounded non-parametric Kernel'};
% Time Unit
list2 = {'Day','Month','Year'};
[indx1] = listdlg('PromptString','Select a Method:',...
                           'SelectionMode','single','ListString',list1,'ListSize',[350 250]);
[indx2] = listdlg('PromptString','Select time unit:',...
                           'SelectionMode','single','ListString',list2,'ListSize',[150 250]);

if indx2==1;prompt = {'\fontsize{12} Magnitude:','\fontsize{12} Period Length, (days):'};
elseif indx2==2;prompt = {'\fontsize{12} Magnitude:','\fontsize{12} Period Length, (months):'};
elseif indx2==3;prompt = {'\fontsize{12} Magnitude:','\fontsize{12} Period Length, (years):'};end

if indx1==1;method='GRU';elseif indx1==2;method='GRT';elseif indx1==3;method='NPU';else;method='NPT';end 

% SET Mmax for truncated Magnitude distribution


if method=='GRT' | method=='NPT'

prompt2 = ['\fontsize{12} M_M_a_x (>',num2str(max(Cmag)),' or estimated from catalog):'];
dlgtitle2 = 'Set Mmax';
dims2 = [1 59]; opts2.Interpreter='tex';
definput2 = {'estimated'};
answer4 = inputdlg(prompt2,dlgtitle2,dims2,definput2,opts2);

if strcmp(answer4,'estimated');Mmax=[];
%questM='Do you wish to include Mmax Bias?';
%answer=questdlg(questM,'Bias Calculation?','Yes','No',opts);

dlgtitleM='Mmax Bias estimation';
promptM={'\fontsize{12} Enter number of synthetic samples (N>1):'}; 
dimsM=[1 57];opts.Interpreter='tex';
definputM={num2str(max(1000))};
answer3=inputdlg(promptM,dlgtitleM,dimsM,definputM,opts);
Nsynth=str2double(answer3);
cd SSH/
if strcmp(method,'GRT')==1 & isempty(Mmax)==1;
[p1,p2,p3,p4,p5,pb,Mmax,err,BIAS,SD]=TruncGR_O(Ctime,Cmag,0,Mc,Mmax,Nsynth);    
Mmax=Mmax+BIAS;
elseif strcmp(method,'NPT')==1 & isempty(Mmax)==1;
[p1,p2,p3,p4,p5,p6,p7,p8,p9,Mmax,err,BIAS,SD]=Nonpar_tr_O(Ctime,Cmag,0,Mc,Mmax,Nsynth);
Mmax=Mmax+BIAS;
end
cd ../
disp([newline 'Corrected Mmax=',num2str(Mmax)]);
else Mmax=str2double(answer4{1});
end

else Mmax=[];
end

% SHA parameters (Magnitude and Time Period)
dlgtitle = 'Parameters for Hazard Analysis';
dims = [1 55]; opts.Interpreter='tex';
definput = {num2str(max(Cmag)),'1'};
answer3 = inputdlg(prompt,dlgtitle,dims,definput,opts);

iop=indx2-1;time_period=str2double(answer3(2));mag=str2double(answer3(1));


%% *************       STEP_7: GENERATE AND SAVE OUTPUTS       ************** 

%% RUN MAGDIST
[HP] = TDHMagDistWrapper(method, time_windows, Mc, iop,Mmax)

%% Harzard Parameters Estimate
[MRPer,ExPr]=TDHRetPeriodExcProbWrapper(method,mag,time_period,Mc,HP)

% Return from Directory with Functions

%% PLOT and Save Results

Zplo

Zsave_output


%% -*-*-*-*-*-*-*-*-*-*-*-*- F U N C T I O N S -*-*-*-*-*-*-*-*-*-*-*-*-
%% ------------------------- DATA HANDLING FUNCTION   --------------------------

function [Catalog,PROD_data,Mode1,dstr1,dstr2,s,ss1,ss,ss2]=Data_Hand_A2M
clc

disp('please select mode:');
disp(' "1" - For Seismic Data');
disp(' "2" - For both Seismic and Production Data');
Mode1=input('mode:  ');

switch Mode1
    
    case 1
%% SEISMIC DATA
cd CATALOGS
d=dir;dstr1={d.name};
[SData,SFields,s,s1]=DatLoad(dstr1,1);
cd ../

Na=length(SFields);if SFields(Na)~=' ';SFields(Na+1)=' ';end
[cou,c,Datime,Catalog]=Fields_dat(SData,SFields,1);
ss1=1:length(Catalog);% ss1=SetParams(Catalog);%Catalog=Catalog(ss1);
PROD_data=[];ss=[];s2=[];ss2=[];dstr2=[];
    case 2
%% ------------------------------------------------------------------------------------------------------------------------------------
% BOTH -SEISMIC and PRODUCTION DATA

cd CATALOGS                                             %Seismic Data
d=dir;dstr1={d.name};
[SData,SFields,s,s1]=DatLoad(dstr1,1);
cd ../

Na=length(SFields);if SFields(Na)~=' ';SFields(Na+1)=' ';end
[cou,c,Datime,Catalog]=Fields_dat(SData,SFields,1);
ss1=1:length(Catalog);
% ss1=SetParams(Catalog);%Catalog=Catalog(ss1);
 
cd PRODUCTION_DATA                         % Production (non-Seismic) Data
d=dir;dstr2={d.name};
[OData,OFields,ss,s2]=DatLoad(dstr2,2);
cd ../

Na1=length(OFields);if OFields(Na1)~=' ';OFields(Na1+1)=' ';end
[cou1,c1,Datime1,PROD_data]=Fields_dat(OData,OFields,2);
ss2=SetParams(PROD_data);%PROD_data=PROD_data(ss2);

end
 % save('Catalog_ST2_Test','Catalog')
 % save('Catalog_ST2_Test','Catalog')


%% -----------------------------------------    F U N C T I O N S    -----------------------------------------
function [cou,c,Datime,OUT]=Fields_dat(indata,infields,iop)

Datime=datenum(indata(:,1),indata(:,2),indata(:,3),indata(:,4),indata(:,5),indata(:,6));       % Convert time to matlab format

% Define Fields
c=1;
for i=1:length(infields)-1
if strcmp(infields(i),' ')==1;cou(i)=0;
else cou(i)=c;end
if strcmp(infields(i),' ')==0 & strcmp(infields(i+1),' ')==1;c=c+1;end
end

if iop==1
OUT(1).field='Occurrence_Time';OUT(1).val=Datime;
elseif iop==2
OUT(1).field='Production_Time';OUT(1).val=Datime;
end

for i=2:c-1
OUT(i).field=infields(cou==i);    
OUT(i).val=indata(:,i+5);
end

%Set Field Type for Magnitude Recognition
for i=1:size(OUT,2)
    if strcmp(OUT(i).field,'ML') || strcmp(OUT(i).field,'Mw') || strcmp(OUT(i).field,'M') ...
   || strcmp(OUT(i).field,'Ms') || strcmp(OUT(i).field,'mb') || strcmp(OUT(i).field,'Md')
         OUT(i).fieldType='Magnitude';
    else
    OUT(i).fieldType=[];
    end
end

end

%% -----------------------------------------------------------------------------------------------------------

function [Data,Fields,s,f]=DatLoad(dstr,da)        
% Select Seismic Catalog
if da==1;dlabel1='Select Seismic Catalog:';dlabel2='Select Catalog FIELDS file:';
elseif da==2;dlabel1='Select Production Data file:';dlabel2='Select Production FIELDS file:';end

[s,ok]=listdlg('PromptString',dlabel1,...
    'SelectionMode','single',...
    'ListString',dstr,'ListSize',[350 250]);
if ok; Data=load(dstr{s});end

[f,ok]=listdlg('PromptString',dlabel2,...
    'SelectionMode','single',...
    'ListString',dstr,'ListSize',[350 250]);
if ok; Fields=fileread(dstr{f});end

end


%%---------------------------------------------------------------------------------------------------------------
function ss1=SetParams(Data)
%Select Parameters from Seismic Catalog
[s1,ok]=listdlg('PromptString','Select field:',...
    'ListString',{Data.field},'SelectionMode','single','ListSize',[350 250]);
for i=1:length(s1)
    ss1=s1(isnan(s1)==0);
end

end
% ----------------------------------------------------------------------------------------------------------------


end
%%----------------------------------------------------------------------------------------------------------------
function [Ctime,Cmag,Mtype]=Select_Magnitude_Scale(Catalog)

cou=1;
for i=1:length(Catalog)
    if strcmp(Catalog(i).fieldType,'Magnitude')==1
        C(cou).field=Catalog(i).field;cou=cou+1;
    end
end

% Check for no magnitude
if cou==1;error('MyComponent:incorrectType',...
'No magnitude column detected!! Please check: Magnitude fields must be noted as one of the following:\nM  ,  Mw  ,  ML  ,  Ms  ,  mb  ,  Md\n or select the entire sample for analysis');end

%Select Parameters from Seismic Catalog -
[ss1,ok]=listdlg('PromptString','Please Select M scale:',...
    'ListString',{C.field}, 'SelectionMode','single');

 id=findfield(Catalog,C(ss1).field);
 id_time=findfield(Catalog,'Occurrence_Time');

 Mtype=Catalog(id).field;
  
Cmag=Catalog(id).val;Ctime=Catalog(id_time).val;

end
  
%% --------------------------------------------------------------------------------------
% finds a field defined by a certain (string) name
function [id] = findfield( catalog,field )
id=0;
j=1;
while j <= size(catalog,2) && id==0
        if (strcmp(catalog(j).field,field)==1)
        id=j;
        end
        j=j+1;
end
end
  
%% -----------!!!!!!!!!!! HAZARD PARAM<ETERS ESTIMATE FUNCTIONS !!!!!!!!!!!-----------

function [HP] = TDHMagDistWrapper(method, time_win_data, mmin, iop,Mmax)
cd SSH
for j=1:size(time_win_data,2);if isempty(time_win_data(j).Time);M1(j)=NaN;else;M1(j)=max(time_win_data(j).M);end;end
Mmaxcat=max(M1);Mmax=max([Mmaxcat Mmax])

	for i=1:size(time_win_data,2)
		mags_vec = time_win_data(i).M;
		time_vec = time_win_data(i).Time;
		HP(i).mmin = mmin;
		HP(i).iop = iop;
		HP(i).method = method;
		switch method
			case 'GRU'
				try
					[HP(i).lamb_all, HP(i).lamb, HP(i).lamb_err, HP(i).unit, HP(i).eps, HP(i).b]=UnlimitGR(time_vec, mags_vec, iop, mmin);
                    HP(i).Mmax=NaN;eps=HP(i).eps;
                    [m, PDF_GRU, CDF_GRU]=dist_GRU(mmin,Mmax+eps,eps,mmin,eps,HP(i).b);                                  %K29JAN2020 - magnitude PDF/CDF 
                    HP(i).PDF=[m PDF_GRU];HP(i).CDF=[m CDF_GRU];                                                         %K29JAN2020 - magnitude PDF/CDF 
				catch err
	                HP(i).lamb_all=NaN; HP(i).lamb=NaN; HP(i).lamb_err=2; HP(i).unit=''; HP(i).eps=NaN; HP(i).b=NaN;HP(i).Mmax=NaN;HP(i).PDF=NaN;HP(i).CDF=NaN;
	                warning('%s: %s', err.identifier, err.message);
				end
			case 'GRT'
				try
					[HP(i).lamb_all, HP(i).lamb, HP(i).lamb_err, HP(i).unit, HP(i).eps, HP(i).b, HP(i).Mmax, HP(i).err]=TruncGR_O(time_vec, mags_vec, iop, mmin,Mmax);
                    eps=HP(i).eps;[m, PDF_GRT, CDF_GRT]=dist_GRT(mmin,Mmax+eps,eps,mmin,eps,HP(i).b,Mmax);               %K29JAN2020 - magnitude PDF/CDF 
                    HP(i).PDF=[m PDF_GRT];HP(i).CDF=[m CDF_GRT];                                                         %K29JAN2020 - magnitude PDF/CDF 
				catch err
	                HP(i).lamb_all=NaN; HP(i).lamb=NaN; HP(i).lamb_err=2; HP(i).unit=''; HP(i).eps=NaN; HP(i).b=NaN; HP(i).Mmax=NaN; HP(i).err=NaN;HP(i).PDF=NaN;HP(i).CDF=NaN;
	                warning('%s: %s', err.identifier, err.message);
				end
			case 'NPU'
				try
					[HP(i).lamb_all, HP(i).lamb, HP(i).lamb_err, HP(i).unit, HP(i).eps, HP(i).ierr, HP(i).h, HP(i).xx, HP(i).ambd]=Nonpar_O(time_vec, mags_vec, iop, mmin);
                    HP(i).Mmax=NaN;eps=HP(i).eps;
                    [m, PDF_NPU, CDF_NPU]=dist_NPU(mmin,Mmax+eps,eps,mmin,eps,HP(i).h,HP(i).xx,HP(i).ambd);              %K29JAN2020 - magnitude PDF/CDF
                    HP(i).PDF=[m PDF_NPU];HP(i).CDF=[m CDF_NPU];                                                         %K29JAN2020 - magnitude PDF/CDF 
               catch err
	               HP(i).lamb_all=NaN; HP(i).lamb=NaN; HP(i).lamb_err=2; HP(i).unit=''; HP(i).eps=NaN; HP(i).ierr=NaN; HP(i).h=NaN; HP(i).xx=[]; HP(i).ambd=[];HP(i).Mmax=NaN;HP(i).PDF=NaN;HP(i).CDF=NaN;
	               warning('%s: %s', err.identifier, err.message);
				end
			case 'NPT'
				try
					[HP(i).lamb_all, HP(i).lamb, HP(i).lamb_err, HP(i).unit, HP(i).eps, HP(i).ierr, HP(i).h, HP(i).xx, HP(i).ambd, HP(i).Mmax, HP(i).err]=Nonpar_tr_O(time_vec, mags_vec, iop, mmin,Mmax);
                    eps=HP(i).eps;
                    [m, PDF_NPT, CDF_NPT]=dist_NPT(mmin,Mmax+eps,eps,mmin,eps,HP(i).h,HP(i).xx,HP(i).ambd,Mmax);         %K29JAN2020 - magnitude PDF/CDF
                    HP(i).PDF=[m PDF_NPT];HP(i).CDF=[m CDF_NPT];                                                         %K29JAN2020 - magnitude PDF/CDF 
                catch err
	                HP(i).lamb_all=NaN; HP(i).lamb=NaN; HP(i).lamb_err=2; HP(i).unit=''; HP(i).eps=NaN; HP(i).ierr=NaN; HP(i).h=NaN; HP(i).xx=[]; HP(i).ambd=[]; HP(i).Mmax=NaN; HP(i).err=NaN;HP(i).PDF=NaN;HP(i).CDF=NaN;
	                warning('%s: %s', err.identifier, err.message);
				end
        end
  
   %   K12NOV2015
  % Calculate lamb and lamb_all in case of 0, 1, or 2 events. 
  % It may be generalized in all cases. Rate is now calculated
  % by division of the event number by the duration of the set 
  % - not the time difference between the first and last events.
   if numel(mags_vec(mags_vec>=mmin))<3  
    HP(i).lamb_all= numel(mags_vec)/(time_win_data(i).Tend-time_win_data(i).Tstart);    
    HP(i).lamb= numel(mags_vec(mags_vec>=mmin))/(time_win_data(i).Tend-time_win_data(i).Tstart);    
       switch iop
           case 0
            %OK   
           case 1
    HP(i).lamb_all=HP(i).lamb_all*30;
    HP(i).lamb=HP(i).lamb*30;
           case 2
    HP(i).lamb_all=HP(i).lamb_all*365;           
    HP(i).lamb=HP(i).lamb*365;
       end
   end   
 %   K12NOV2015

    end

cd ../
    
end

%% 
function [MRPer,ExPr]=TDHRetPeriodExcProbWrapper(meth,mag,time_period,Mmin,HP)
        nn=size(HP,2);
        if nn==0; 'All datasets are empty'
                MRPer=[]; ExPr=[];
        else
                for i=1:nn
                    Mmax=HP(i).Mmax;
                        if HP(i).lamb_err==0;
							try
                        		[MRPer(i),ExPr(i)]=SingleRetPeriodExcProbWrapper(meth,mag,time_period,Mmin,HP(i),Mmax);
							catch err
				                MRPer(i)=NaN; ExPr(i)=NaN;
				                warning('%s: %s', err.identifier, err.message);cd ../
							end
                        else 
                        	MRPer(i)=NaN; ExPr(i)=NaN;
                        if MRPer(i)==inf;MRPer(i)=NaN;end   %K 21OCT2016
                        end
                end
        end
end

%% Function - SingleRetPeriodExcProbWrapper

function [MRPer,ExPr]=SingleRetPeriodExcProbWrapper(meth,mag,time_period,Mmin,HP,Mmax)
cd SSH
	%Md=3;Mu=3;dM=3;Mmin=0.6;
	Md = mag;
	Mu = mag;
	dM = mag;
	
	if isnan(Mmax) && isfield(HP,'Mmax')
		Mmax = HP.Mmax;
	end
	
	m = []; rper = []; prob = [];
	
	switch meth
	    case 'GRU'
	    	[m,rper]=Ret_periodGRU(Md,Mu,dM,Mmin,HP.lamb,HP.eps,HP.b);
	    	[m,prob]=ExcProbGRU(0,Md,Mu,dM,time_period,Mmin,HP.lamb,HP.eps,HP.b);
	    case 'GRT'
	    	[m,rper]=Ret_periodGRT(Md,Mu,dM,Mmin,HP.lamb,HP.eps,HP.b,Mmax);
	    	[m,prob]=ExcProbGRT(0,Md,Mu,dM,time_period,Mmin,HP.lamb,HP.eps,HP.b,Mmax);
	    case 'NPU'    
	    	[m,rper]=Ret_periodNPU(Md,Mu,dM,Mmin,HP.lamb,HP.eps,HP.h,HP.xx,HP.ambd);
	    	[m,prob]=ExcProbNPU(0,Md,Mu,dM,time_period,Mmin,HP.lamb,HP.eps,HP.h,HP.xx,HP.ambd); 
	    case 'NPT'  
	    	[m,rper]=Ret_periodNPT(Md,Mu,dM,Mmin,HP.lamb,HP.eps,HP.h,HP.xx,HP.ambd,Mmax);
	    	[m,prob]=ExcProbNPT(0,Md,Mu,dM,time_period,Mmin,HP.lamb,HP.eps,HP.h,HP.xx,HP.ambd,Mmax);
	end

	if isempty(rper)
		MRPer = NaN;
		ExPr = NaN;
	else
		MRPer = rper(1);
		ExPr = prob(1);
	end
cd ../
end

  


