% FUNCTION:        T2ED
% VERSION:          [Wrapper Standalone Version]  V2.8
% COMPATIBLE with Matlab version 2017b or later
% REQUIREMENTS: "Statistics and Machine Learning" Matlab Toolbox installed
% TOOLBOX: "Clustering/Transformation to ED Toolbox" within SERA Project
% DOCUMENT: "READ_ME_App_1A_v2_Description_T2ED.docx"
% ----------------------------------------------------------------------------------------------------------------
% Transformation of multiple parameters into Equavalent Dimensions
% ---------------------------------------------------------------------------------------------------------------
% OVERVIEW: The function takes as input seismic or/and operational data 
% parameters and transform them into their equivalent dimensions following 
% the methodology introduced by Lasocki (2014).Please check also the scripts
% 'ED_ToolBox_wrapper' and 'plottesth' for specified applications (scenarios)
% and plotting results.
% -----------------------------------------------------------------------------------------------------------------
% AUTHORS: K. Leptokaropoulos, S. Lasocki and P. Urban 
% based on code 'EQ_DIM' originally developed by S. Lasocki (Lasocki, 2014)
% Last Updated: 07/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 documentation:
%                                                                    "READ_ME_App_1A_v2_Description_T2ED.docx" 
% for description of the Application and its requirements.
% -----------------------------------------------------------------------------------------------------------------
% DESCRIPTION: This function constitutes the main part of the "Clustering/
% Transformation to ED" Toolbox which is performed in a Standalone Mode
% taking as input data in ASCII format (See Input Data Requirements Specification
% in the Application Documentation). Transformation to Equivalent Dimensions 
% of Seismic and /or Operation data parameters is ultimately achieved. This 
% process is adequately described in Lasocki, [2014] and also in the comments 
% in the 'EQ_DIM' function. A brief description of the preparation process is 
% given here (found in detail in 'READ_ME_App_1A_v2_Description_T2ED'). 
% These  following steps are followed within 'ED_ToolBox_Wrapper.m' script, where
% the User can define/select the input data & parameters (function arguments):
%   STEP 1- Data Uploading                                                 
%   STEP 2- Data Handling and Conversion                          ['Data_Hand_A2M.m' script]           
%   STEP 3- Time and Magnitude scale Columns Importing                       (Seismic Data)
%   STEP 4- Mc filtering                     
%   STEP 5- Date Column Importing                                                (Production Data)
%   STEP 6- Data Selection                     (define vectors for analysis and reshape data)
%   STEP 7- Parameters values Selection                                    (for ED_DIM function)
%   STEP 8- Selection of Mode for Background/Testing Data              
%   STEP 9- Defining Background/Testing Data               
%   STEP 10- RUN 'T2ED' PROGRAM (this script) after selecting a Scenario
%   STEP 11- DISPLAY (optional) AND SAVE OUTPUTS:                                        
%  ----------------------------------------------------------------------------------------------------------------
% INPUT: [Documentation is available in ['App_1A_v3_Description_T2ED.doc'].                
%                - SEIS_data     : Seismic Catalog                                       (ASCII format)
%                - PROD_data    : Production Data                                       (ASCII format)
%                - Ctime           : Time Column                                            (from Catallog)
%                - PRODtime      : Time (Date) Column                  (from Production Data File)
%                - randmode      : Sample Randomization Mode                        (See EQ_DIM)
%                - sample_mult  : Sample Multiplication Mode                           (See EQ_DIM)
%                - DT               : Time lag         (i.e. response of Seismic to Operational data) 
%                - option          : Mode for picking Background and Testing data:
%                                             - 1, for Events          [only when S. Catalog is selected]
%                                             - 2, for Time 
%                                             - 3, for Using the entire datasets
%                - BG                 : Background Data picking
%                - TS                  : Testing Data picking
%                - Mc                 : Completeness Magnitude                  
% -------------------------------------------------------------------------------------------------------------------
% OUTPUT:
%                 * Output Report with data and parameters used 
%                 * Tdata: Output Structure, with input (Original) as well as output 
%                               (Transformed) data. The structure has as many cells as the 
%                               number of selected by User parameters.The structure fields 
%                               are vectors unless stated otherwise. These fields are:
%                    - xt               : The transformed Testing data (parameter values in 
%                                            the Equivalent Dimensions, [0 1])
%                    - xB              : The transformed Background data (parameter values 
%                                            in the Equivalent Dimensions, [0 1])
%                    - ierr            : [scalar] h-convergence indicator  (See EQ_DIM for details)
%                    - h                : [scalar] kernel smoothing factor
%              (1) - xx              : background sample considered for transformation
%              (1) - ambd        : weigthing factors for the adaptive kernel
%                    - field          : [string] indicating the corresponding field 
%                                          (i.e. the transformed parameter)
%              - index_Testing: Indicator of Testing Data from the original Dataset 
%                                          (see ‘origval_all' field) which were transformed  
%       - index_Background: Indicator of Background Data from the original
%                                          Dataset (see ‘origval_all' field) which were transformed 
%                    - all              : transformed parameters vector with size of the original 
%                                          parameter vector (‘origval_all’, including NaN's)
%                    - origval      : vector with the original parameters that were transformed
%                                           (without NaN's)
%                    - origval_all: original input vector with the event parameters that were 
%                                           transformed (including NaN's). (Note that when both 
%                                           seismic and Production (non-seismic) data are used, 
%                                           the Production data are interpolated as event parameters.
%                    - Source      : [string] Source of the data, either 'Catalog', or 'Production'
% (1) - These parameters may be doubled or tripled compared to xB, according to the 
%         selected “sampl_mult” value
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -                  
%                 * Output Figures: histograms with original and transformed background &
%                                                 full samples as well as adaptive kernel weighting factors.   
%                                                 See 'plottesth' for example - option can be On/Off
% --------------------------------------------------------------------------------------------------------------------
% REFERENCES:
% Lasocki S. (2014), Transformation to equivalent dimension - a new methodology to
% study earthquake clustering, Geophys. J. Int., 197,1224-1235,doi:10.1093/gji/ggu062
% --------------------------------------------------------------------------------------------------------------------
% 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 [Tdata]=T2ED_V2_8(SEIS_data,PROD_data,Ctime,PRODtime,randmode,sample_mult,DT,option,BG,TS,Mc)


if isempty(SEIS_data) && option=='1';error('Events Mode is not applicable for Production Data, please choose option="2", or use a Seismic Catalog');end
    if isempty(BG)==0;P1=BG(1);P2=BG(2);Pt1=TS(1);Pt2=TS(2);else
    P1=[];P2=BG;Pt1=BG;Pt2=BG;end
if option=='1';
    t1=Ctime(BG(1));t2=Ctime(BG(2));tt1=Ctime(TS(1));tt2=Ctime(TS(2));
end 


% Create new Structure with necessary fields - Catalog
if isempty(SEIS_data);Data=[];Mc='not applicable';else
[Data]=homogenizeD(SEIS_data);end

% Create new Structure with necessary fields - Production data
if isempty(PROD_data);Data2=[];else
[Data2]=homogenizeD(PROD_data);end



% ********   STEP_2: DATA SELECTION & FILTERING   ********** 

% load Seismic Data (Catalog)  

% load Production Data (non-seismic)

   % Select Mc and filter parameters for M>=Mc 
   %[Ctime,Cmag,Data,Mc]=FiltMc(Catalog,Cvector);

   % Treat Production data
   %[Data2]=NonS_Data(d,Fields,FieldUnit,PRODvector)

   %Select Parameters for T2ED
 
   
   %MergeData (Seismic/ Non-Seismic)
   if isempty(Data)==0 && isempty(Data2)==0
   [Data3,Data2]=MergeSNSdata(PRODtime,Ctime,DT,Data,Data2);
   elseif isempty(Data)==0 && isempty(Data2)==1
   Data3=Data;DT='not applicable';
   elseif isempty(Data)==1 && isempty(Data2)==0
   Data3=Data2;DT='not applicable';end
   
  %% _______________________________________________________%% 


% ********    STEP_3: SELECT BACKGROUND SAMPLE   ********
% disp('Select Starting and Ending Point for Background And Testing Data');
% if Mode1==2;disp(' - "1" for typing observations');else;disp(' - "1" for typing events');end
% disp(' - "2" for typing time');
% disp(' - "3" for selecting from plot');
% disp(' - "0" for using the entire data as background sample');
% BGD=input(' ');

% [option,P1,P2,t1,t2]=BcGround_Sample(BGD,Mode1,Ctime,PRODtime,Data2,s2,Description,1);
% [option,Pt1,Pt2,tt1,tt2]=BcGround_Sample(BGD,Mode1,Ctime,PRODtime,Data2,s2,Description,2);

% DISPLAY parameters report in the screen
           % parameters_report(Mode1,dstr,s,ss,s1,s2,Catalog,FieldDescription,randmode,sample_mult,Mc,DT,t1,t2,tt1,tt2);

% **********     STEP_4: RUN EQ_DIM PROGRAM    **********
% Run T2ED in loop
if isempty(Data);D=Data2;t=PRODtime;else;D=Data3;t=Ctime;end
nn1=length(D);
for i=1:nn1
x=D(i).val;    
[xt,xBG,ierr,h,xx,ambd,IX1,IX2]=T2ED_Wrapper_MODE1(option,P1,P2,Pt1,Pt2,x,randmode,sample_mult,D(i),t(D(i).index),D,i);
Tdata(i).xt=xt;
Tdata(i).xBG=xBG;
Tdata(i).ierr=ierr;
Tdata(i).h=h;
Tdata(i).xx=xx;
Tdata(i).ambd=ambd;
Tdata(i).field=(['Transformed_',D(i).field]);
Tdata(i).index_Testing=IX1;
Tdata(i).index_Background=IX2;
Tdata(i).all=nan(size(D(i).index));Tdata(i).all(IX1)=xt;Tdata(i).all(IX2)=xBG;
Tdata(i).origval=x;
Tdata(i).origval_all=D(i).all;
Tdata(i).Source=D(i).source;
%clearvars xt ierr h xx ambdb
end

% *********   STEP_5: DISPLAY AND SAVE OUTPUTS   ********
% Save Outputs
if option=='1'
SaveOuts(Data3,randmode,sample_mult,Mc,DT,t1,t2,tt1,tt2,Tdata)
else SaveOuts(Data3,randmode,sample_mult,Mc,DT,P1,P2,Pt1,Pt2,Tdata)
end

%% ASCII OUTPUT FILE
cd Outputs_ED\
fid0=fopen('Data_Original.txt','w');
fid1=fopen('Data_Transformed.txt','w');
fid2=fopen('Data_fields.txt','w');
for i=1:size(Tdata,2)
D_O(:,i)=Tdata(i).origval_all;
D_T(:,i)=Tdata(i).all;
D_F{i}=Tdata(i).field;
fprintf(fid2,'%s    ',D_F{i});
end
fprintf(fid0,[repmat('%36.16f      ',[1,size(D_O,2)]),'\n'],D_O');
fprintf(fid1,[repmat('%12.10f    ',[1,size(D_T,2)]),'\n'],D_T');
fclose(fid0);fclose(fid1);fclose(fid2);cd ../

disp(' ')
disp('Please Find the Results in Folder: "Outputs_ED" ')


% *************************** FUNCTIONS ***************************
% ****-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-****

%% HOMOGENIZE Catalog-Production files
function [Data]=homogenizeD(INdata)
cou=1;
    for i=1:size(INdata,2)
        x=INdata(i).val;
        index=isnan(x)==0;
        x=x(isnan(x)==0);
        Data(cou).val=x;
        Data(cou).field=INdata(i).field;
        
        Data(cou).source=INdata(i).source;
        Data(cou).index=index;
        Data(cou).all=nan(size(Data(cou).index)); Data(cou).all(index)=x;
        cou=cou+1;
    end
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

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

%% --------------------------------------------------------------------------------------
% Merge Seismic. Non-Seismic Data
function [Data,Data2]=MergeSNSdata(PRODtime,Ctime,DT,Data,Data2);
% TRANSFORM TIME to adjust to EQ Occurrence 
s2=1:size(Data2,2);  

TT=PRODtime+DT;
T=RandomizeT(TT);
cou=1;
for i=1:length(s2)
    Data2(i).all=interp1(T,Data2(i).all,Ctime);
    x=Data2(i).all;
    index=isnan(x)==0;
    x=x(isnan(x)==0);
    Data2(cou).index=index;
    Data2(cou).Tval=x;
    Data2(cou).Tval_all=nan(size(Data2(cou).index)); Data2(cou).Tval_all(index)=x;
    cou=cou+1;
end


% Merge Seismic and Non-Seismic Parameters
N1=size(Data,2);N2=size(Data2,2);
for i=1:N2
    Data(N1+i).field=Data2(i).field;
    Data(N1+i).index=Data2(i).index;
    
    Data(N1+i).all=Data2(i).Tval_all;            %watch out for the inversed parameters!!!!
    Data(N1+i).val=Data2(i).Tval;                %watch out for the inversed parameters!!!!
    Data(N1+i).source=Data2(i).source;
end

end

%% --------------------------------------------------------------------------------------
% integer time points randomization to float
% For interpolation to occur, randomization of grid vector is necessary
% Uniform distribution is used for randomization within eps/10 accuracy
function [T]=RandomizeT(T1)

dT1=T1(2:length(T1))-T1(1:length(T1)-1);dT1=sort(dT1);
DT1=dT1(dT1>0);eps=min(DT1)/10;
ra=rand(length(T1),1);
T=T1-0.5*eps+eps.*ra;

end

%% --------------------------------------------------------------------------------------
function [xt,xBG,ierr,h,xx,ambd,IX1,IX2]=T2ED_Wrapper_MODE1(option,P1,P2,Pt1,Pt2,x,randmode,sample_mult,Data,t,D,i)

% STEP 1 & STEP 2 - MODE SELECTION & DATA PREPARATION

% STEP 3 - PARAMETERS SELECTION 
% A) BACKGROUND DATA OPTIONS:
% 'Points' - Option='1'
% 'Time' - 'Option='2'

%B)  -------------------------------- Parameters Value List --------------------------------  
% B1)  randmode - string determining the mode of data randomization (3rd input): 
%           'no'   - no randomization, 
%           'exp'  - exponential randomization, 
%           'norm' - Gaussian randomization,
%           'uni'  - uniform randomization
% B2) sample_mult - string determining the way of sample multiplication (4th input): 
%           'no'   - no multiplication, 
%           'left' - dubbling from the left hand side, 
%           'right - doubbling from the right hand side, 
%           'both' - trippling from both sides.

dix=find(Data.index==1);

switch option
    
    case  '1'
 n1=(dix<P1);n1=sum(n1)+1;    
 n2=(dix<=P2);n2=sum(n2);
 nt1=(dix<Pt1);nt1=sum(nt1)+1;    
 nt2=(dix<=Pt2);nt2=sum(nt2);

%if n2>length(x);n2=length(x);end   

    case '2'
        
 T1=P1;
 T2=P2; 
 p1=T1-t;p2=T2-t;
 n1=find(p1<0);n2=find(p2>0);
 n1=min(n1);n2=max(n2);

 Tt1=Pt1;
 Tt2=Pt2; 
 pt1=Tt1-t;pt2=Tt2-t;
 nt1=find(pt1<0);nt2=find(pt2>0);
 nt1=min(nt1);nt2=max(nt2);

 case '3'
    n1=1;n2=length(dix);nt1=n1;nt2=n2;
    
end

IX1=dix(nt1:nt2);
IX2=dix(n1:n2);

%%%% ERROR HANDLING

%if n2>length(x); error(['Please select n2 smaller than ',num2str(length(x)),' events']);
%elseif isempty(x(n1:n2)); error('No events detected! Please select other period');end

if isempty(x(n1:n2)); warning('No events detected! Please select other period');
    xt=nan;xBG=nan;ierr=nan;h=nan;xx=nan;ambd=nan;
elseif isempty(x(nt1:nt2)); warning('No events detected! Please select other period');
   xt=nan;xBG=nan;ierr=nan;h=nan;xx=nan;ambd=nan;
else 
% STEP 4 - RUN T2ED PROGRAM
[xt,xBG,ierr,h,xx,ambd]=EQ_DIMv_2_0(x,n1,n2,nt1,nt2,randmode,sample_mult);

%dealing with sample multiplication
if strcmp(sample_mult,'left')==1;N=length(xBG);xBG=xBG(N/2+1:N);
elseif strcmp(sample_mult,'right')==1;N=length(xBG);xBG=xBG(N/2+1:N);
elseif strcmp(sample_mult,'both')==1;N=length(xBG);xBG=xBG(N/3+1:2*N/3);end    

%round values out of bounds
 xBG(xBG>1)=1;xt(xt>1)=1;
 xBG(xBG<0)=0;xt(xt<0)=0;

%[yt,ierry,hy,xxy,ambdy]=EQ_DIM(M,n1,n2,randmode,sample_mult);

% STEP 5 - DISPLAY OUTPUT WARNING 
% Output  Warning: for h-convergence indicator.
if ierr==0   %warning('Process has converged')                     %16APR2019 - no need to display comment in this case
elseif ierr==1; warning('Multiple h were encountered - largest h is accepted')
elseif ierr==2; warning('Process did not converge - approximate h value is taken')
end
end

% STEP 6 - STORE OUTPUT PARAMETERS 
% see doc for details:
% xt, ierrx, hx,xxx,ambdx

% STEP 7 - OUTPUT VISUALIZATIONS 
% see doc for details:

% STEP 8 - REPEAT PROCESS FOR THE OTHER PARAMETERS (VECTORS) 
% see doc for details:
% optional plot for check 

end

%% --------------------------------------------------------------------------------------
% Save Outputs
function SaveOuts(Data,randmode,sample_mult,Mc,DT,t1,t2,tt1,tt2,Tdata)
% ---- Save *.txt file with Parameters Report ----
mkdir Outputs_ED/
cd Outputs_ED/
fid=fopen('OUTPUT_REPORT.txt','w');
fprintf(fid,['SELECTED DATA & PARAMETERS REPORT (created on ', datestr(now),')\n']);
fprintf(fid,'-----------------------------------------------------------------\n');

fprintf(fid,['<Parameters Selected for Transformation to ED >:\n']);
s2=size(Data,2);

for i=1:s2
    fprintf(fid,['                                              -  ',Data(i).field,'    (',Data(i).source,')\n']);
end
fprintf(fid,'-----------------------------------------------------------------\n');
    if Mc>-20
fprintf(fid,['<Completeness Magnitude                       >: ',num2str(Mc,'%3.1f'),'\n']);
    else 
fprintf(fid,['<Completeness Magnitude                       >: ',num2str('All magnitudes included'),'\n']);
    end
fprintf(fid,['<Time Delay (Operations-Seismicity [days])    >: ', num2str(DT),'\n']);
fprintf(fid,['<Randomization Mode                           >: ', char(randmode),'\n']);
fprintf(fid,['<Sample Multiplication Mode                   >: ', char(sample_mult),'\n']);
if isempty(t1)
fprintf(fid,['<Starting Poing for background data           >: All Data Selected\n']);
fprintf(fid,['<Ending Poing for background data             >: All Data Selected\n']);
fprintf(fid,['<Starting Poing for testing data              >: All Data Selected\n']);
fprintf(fid,['<Ending Poing for testing data                >: All Data Selected\n']);
else
fprintf(fid,['<Starting Poing for background data           >: ', datestr(t1),'\n']);
fprintf(fid,['<Ending Poing for background data             >: ', datestr(t2),'\n']);
fprintf(fid,['<Starting Poing for testing data              >: ', datestr(tt1),'\n']);
fprintf(fid,['<Ending Poing for testing data                >: ', datestr(tt2),'\n']);
end
fclose(fid);

% Save Output Structure
save('Tdata.mat','Tdata')
cd ../

end





end