classdef featureMatrices < handle

% featureMatrices object for ECGwrapper
% -------------------------------------
% 
% Description:
% 
% 
% Arguments: (specified as ECGwrapper('arg_name1', arg_val1, ... , 'arg_nameN', arg_valN) )
% 
%     a. ECG specification:
% 
%         a.1 Recording filename where ECG is stored in a valid format.
% 
%           + recording_name: ECG recording to be classified.
%           + recording_format : Valid ECG format. (MIT, ISHNE, AHA, HES, MAT)
% 
%     b. Operating modes
% 
%       
%     c. Modifiers
% 
% 
% Output:
% 
% Examples:
% 
% 
% Limits and Known bugs:
%   Probably a lot :( ... but dont panic! send me feedback if you need
%   help.
% 
% Author: Mariano Llamedo Soria (llamedom at {electron.frba.utn.edu.ar; unizar.es}
% Version: 0.1 beta
% Birthdate  : 20/3/2012
% Last update: 20/3/2012
       
    properties(GetAccess = private, Constant)
        sampling_rate = 360; %Hz
        fnPayload = {'featMat_clust' 'featMat_ldc' 'ECG' 'QRS_locations' 'RR_intervals'};
        fnFMfile = {'featMat_clust' 'featMat_ldc'};
        fnECGfile = {'ECG' 'QRS_locations' 'RR_intervals'};
        SmallValue = 0.0001;
    end
    
    properties ( Access = private )
        delayLP
        autovec
        q_filters
        LP
        LPsize
        scale_idx
        scales
        global_struct
    end
    
    properties (SetAccess = private, GetAccess = public)  
        %read-only 
        name
        class_labeling
        true_labels
        %auxiliary indexes
        file_count = 1;
        File_list = [];
        File_idx = [];
        QRS_idx = [];
        
    end

    properties
        %read-write
        recording_name
        recording_format
        progress_hdl
        tmp_path
        bCache
        cant_pids
        this_pid
        cant_iteration
        this_iteration
    end
    

    methods 
        
        function obj = featureMatrices(class_labeling)

            %% Constants and definitions

            obj.name = 'feat_matrix_construction';
            obj.class_labeling = class_labeling;

        end
        
        function Start(obj)
        

            %% Calculate accesories signals
            
            %load low-pass preprocessing filter.
            load( ['LP Fs ' num2str(obj.sampling_rate) 'Hz - Fc 35Hz - 80db stop.mat' ]);
            obj.LPsize = length(LP.numerator);
            obj.delayLP = round((obj.LPsize-1)/2);
            obj.LP = LP;
            clear LP
            
            % Wavelet transform filter bank design
            obj.scales = 1:6;
            CantScales = length(obj.scales);
            MaxScales = max(obj.scales);
            obj.scale_idx= nan(MaxScales,1);

            filters_cache_filename = ['wt_filters_' num2str(MaxScales) ' scales_' num2str(obj.sampling_rate) ' Hz.mat' ];
            if( exist(filters_cache_filename, 'file') )
                aux = load( filters_cache_filename );
                obj.q_filters = aux.q_filters;
                clear aux
            else
                obj.q_filters = qs_filter_design(MaxScales, obj.sampling_rate);
            end
            
            for ii = 1:CantScales
                %indice para saber que escala corresponde a cada columna de MyWavelets
                obj.scale_idx(obj.scales(ii)) = ii;
            end

            % Update point
            obj.progress_hdl.checkpoint('PCA projection');

%             obj.autovec = PCA_proj_basis(obj.recording_name, obj.recording_format, obj.q_filters);
            
            obj.autovec = [-0.696609733483207 -0.403719836632163 -0.593081084444745;-0.487887339191286 -0.339521516938205 0.804171053814316;-0.526023595928265 0.849550135686930 0.0395441965529324;];
            
            % data structure to keep data through all iterations.
            obj.global_struct = [];
            
        end
        
        function payload = Process(obj, ECG, this_header, this_iter_ECG_relative_start_end_idx, this_iter_annotations)

            obj.progress_hdl.checkpoint('Label parsing');
            
            [QRS_locations, obj.true_labels ] = Annotation_process(this_iter_annotations, obj.recording_format, obj.class_labeling);
            
            % RR interval sequence
            RR_intervals = diff(QRS_locations, 1);
            RR_intervals = colvec([ RR_intervals(1); RR_intervals ]);

            dRR_intervals = diff(RR_intervals, 1);
            dRR_intervals = colvec([ dRR_intervals(2); dRR_intervals(2); dRR_intervals(2:end) ]);

            % interval resampling.
            RR_intervals = RR_intervals * obj.sampling_rate / this_header.freq;
            dRR_intervals = dRR_intervals * obj.sampling_rate / this_header.freq;
            
            %% Preprocessing
            
            % Update point
            obj.progress_hdl.checkpoint('Resampling');
            
            cant_QRS_locations = length(QRS_locations);
            this_iter_QRS_seq_idx = 1:cant_QRS_locations;

            %resample to obj.sampling_rate Hz
            ECG = resample(ECG, obj.sampling_rate, this_header.freq);
            QRS_locations = colvec(round(QRS_locations * obj.sampling_rate / this_header.freq));
        %         this_iter_true_labels = true_labels(this_iter_QRS_seq_idx);
            this_iter_ECG_resampled_size = size(ECG,1);

            if( this_header.nsig > 2 ) 
                %multilead approach.
                % project ECG and wtECG to obtain PCA components
                ECG = ECG * obj.autovec(:,1:2);
        %             wtECG = cellfun(@(a)( a * obj.autovec(:,1:2) ), mat2cell(wtECG, this_iter_ECG_resampled_size, obj.ECG_header.nsig, ones(1,CantScales) ), 'UniformOutput', false);
        %             wtECG = cell2mat(wtECG);
            end

            % Update point
            obj.progress_hdl.checkpoint('Filtering');

            %Low pass filtering @ 35Hz => ECG recovery
            ECG = filter(obj.LP, ECG);
            %delay compensation.
            ECG = [ zeros(obj.delayLP, size(ECG,2) ) ;ECG(obj.LPsize:end,:) ; zeros(obj.delayLP, size(ECG,2))];

            %Quito la linea de base.
        %     ECG = BaselineWanderRemovalMedian( ECG, obj.sampling_rate);
            ECG = BaselineWanderRemovalSplines( ECG, QRS_locations, obj.sampling_rate);

            % Update point
            obj.progress_hdl.checkpoint('Wavelet transform calculation');

            wtECG = qs_wt(ECG, obj.scales, obj.sampling_rate, obj.q_filters);

            %% Features Calculation

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            this_iter_seq_idx = 1:length(QRS_locations);

            % log(RRactual)
            %%%%%%%%%%%%%%%
            featMat_ldc = log(colvec(RR_intervals(this_iter_QRS_seq_idx))/obj.sampling_rate);

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(RRpost)
            %%%%%%%%%%%%%%%

            featMat_ldc = [ featMat_ldc log(colvec(RR_intervals(min(cant_QRS_locations, this_iter_QRS_seq_idx+1)))/obj.sampling_rate)];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(RRmean1)
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)( find(  QRS_locations >= (QRS_locations(a) - 1*60*obj.sampling_rate) & ... 
                                            QRS_locations <= QRS_locations(a) )), ...
                               this_iter_seq_idx, 'UniformOutput', false);

            aux_featval = cellfun(@(a)(mean(RR_intervals(a))/obj.sampling_rate), aux_idx);
            featMat_ldc = [ featMat_ldc colvec(log(aux_featval)) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(RRmean20)
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)( find(  QRS_locations >= (QRS_locations(a) - 20*60*obj.sampling_rate) & ... 
                                            QRS_locations <= QRS_locations(a) )), ...
                               this_iter_QRS_seq_idx, 'UniformOutput', false);

            aux_featval = cellfun(@(a)(mean(RR_intervals(a))/obj.sampling_rate), aux_idx);
            featMat_ldc = [ featMat_ldc colvec(log(aux_featval)) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % AutoCorr1PlanoWTZeroCross
            %%%%%%%%%%%%%%%
            % AutoCorr1PlanoWTModMaxPos
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)( max(1, QRS_locations(a) - round(0.08*obj.sampling_rate)): ...
                                     min(this_iter_ECG_resampled_size, QRS_locations(a) + round(0.08*obj.sampling_rate))) , ...
                               this_iter_seq_idx, 'UniformOutput', false);

            aux_featval = cellfun(@(a)(squeeze(wtECG(a,:,obj.scale_idx(4)))), aux_idx, 'UniformOutput', false);

            %calculate PCA matrix in this slices for feature
            %AutoCorr1PlanoWTModMaxPos
            autovec_slices = cellfun(@(a)(autovec_calculation(a)), aux_featval, 'UniformOutput', false );

            aux_idx = arrayfun(@(a)( max(1, QRS_locations(a) - round(0.13*obj.sampling_rate)): ...
                                     min(this_iter_ECG_resampled_size, QRS_locations(a) + round(0.2*obj.sampling_rate))) , ...
                               this_iter_seq_idx, 'UniformOutput', false);

            aux_featval = cellfun(@(a,b)(squeeze(wtECG(a,:,obj.scale_idx(4))) * b), aux_idx, autovec_slices, 'UniformOutput', false);
            [ aux_mp aux_zc ] = cellfun(@(a)(CalcModMaxPos(a(:,1))), aux_featval );

            featMat_ldc = [ featMat_ldc colvec(aux_zc) colvec(aux_mp) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % AutoCorr2PlanoWTZeroCross
            %%%%%%%%%%%%%%%
            % AutoCorr2PlanoWTModMaxPos
            %%%%%%%%%%%%%%%

            [ aux_mp aux_zc ] = cellfun(@(a)(CalcModMaxPos(a(:,2))), aux_featval );

            featMat_ldc = [ featMat_ldc colvec(aux_zc) colvec(aux_mp) ];


            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            %clear some auxiliar variables.
            clear aux*

            %calculate clustering features.


            % log(RRactual)
            %%%%%%%%%%%%%%%
            featMat_clust = featMat_ldc(:,1);

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(RRant)
            %%%%%%%%%%%%%%%
            featMat_clust = [ featMat_clust log(colvec(RR_intervals(max(1, this_iter_QRS_seq_idx-1)))/obj.sampling_rate)];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(Prematuridad_loc)
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)(max(1,a-1):min(cant_QRS_locations,a+1)), this_iter_QRS_seq_idx, 'UniformOutput', false);
            %tengo que contemplar los casos extremos
            if( length(aux_idx{1}) < 3 )
                aux_idx{1} = [1 aux_idx{1}];
            end
            if( length(aux_idx{end}) < 3 )
                aux_idx{end} = [aux_idx{end} cant_QRS_locations];
            end

            aux_featval = cellfun(@(a)( RR_intervals(a(2))/sum(RR_intervals(a)) ), aux_idx);
            featMat_clust = [ featMat_clust colvec(log(aux_featval)) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(AbsRRlocalVar)
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)(max(1,a-1):min(cant_QRS_locations,a+1)), this_iter_QRS_seq_idx, 'UniformOutput', false);
            if( length(aux_idx{1}) < 3 )
                aux_idx{1} = [1 aux_idx{1}];
            end
            if( length(aux_idx{end}) < 3 )
                aux_idx{end} = [aux_idx{end} cant_QRS_locations];
            end

            aux_featval = cellfun(@(a)(obj.SmallValue+sum(abs(dRR_intervals(a)))/obj.sampling_rate), aux_idx);

            featMat_clust = [ featMat_clust colvec(log(aux_featval)) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(RRmean20)
            %%%%%%%%%%%%%%%

            featMat_clust = [ featMat_clust featMat_ldc(:,4) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % log(QRS_max_scale_proj12)
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)( max(1, QRS_locations(a) - round(0.08*obj.sampling_rate)): ...
                                     min(this_iter_ECG_resampled_size, QRS_locations(a) + round(0.08*obj.sampling_rate))) , ...
                               this_iter_seq_idx, 'UniformOutput', false);

            aux_featval = cellfun(@(a)(wtECG(a,:,:)), aux_idx, 'UniformOutput', false);

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            aux_featval = cellfun(@(a)(squeeze(sum(abs(a(:,1,:))))'), aux_featval, 'UniformOutput', false );
            aux_featval = cellfun(@(a)(obj.scales * colvec(a) ./ sum(a)), aux_featval );

        %     iAreaAbs = squeeze(sum(abs(wtAux)))';
        %     MaxProjArea = obj.scales * iAreaAbs ./ sum(iAreaAbs);

            featMat_clust = [ featMat_clust colvec(log(aux_featval)) ];

            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % AutoCorr1PlanoWTModMaxPos
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)( max(1, QRS_locations(a) - round(0.13*obj.sampling_rate)): ...
                                     min(this_iter_ECG_resampled_size, QRS_locations(a) + round(0.2*obj.sampling_rate))) , ...
                               this_iter_seq_idx, 'UniformOutput', false);

            aux_featval = cellfun(@(a,b)(squeeze(wtECG(a,:,obj.scale_idx(4))) * b), aux_idx, autovec_slices, 'UniformOutput', false);
            aux_featval = cellfun(@(a)(CalcModMaxPos(a(:,1))), aux_featval);

            featMat_clust = [ featMat_clust colvec(aux_featval) ];


            % Update point
            obj.progress_hdl.checkpoint('Features calculation');

            % CrossCorr_QRST_ModMaxVal13
            %%%%%%%%%%%%%%%

            aux_idx = arrayfun(@(a)( max(1, QRS_locations(a) - round(0.13*obj.sampling_rate)): ...
                                     min( [ this_iter_ECG_resampled_size, ...
                                            QRS_locations(a) + round(0.4*obj.sampling_rate), ...
                                            QRS_locations(a) + RR_intervals(a) - round(0.13*obj.sampling_rate) ] ) ), ...
                               this_iter_seq_idx, 'UniformOutput', false);


            if( this_header.nsig > 2 ) 
                % wtECG already projected in obj.autovec.
                aux_featval = cellfun(@(a)( CalcModMaxVal( squeeze(wtECG(a,:,obj.scale_idx(3))) ) ), aux_idx);
            else
                aux_featval = cellfun(@(a)( CalcModMaxVal( squeeze(wtECG(a,:,obj.scale_idx(3))) * obj.autovec ) ), aux_idx);
            end

            featMat_clust = [ featMat_clust colvec(aux_featval) ];

            %save ECG for user interface and cluster analysis
            for ii = 1:length(obj.fnPayload)
                payload.(obj.fnPayload{ii}) = eval(obj.fnPayload{ii});
            end
    
        end
        
        function allPayloads = Concatenate(obj, allPayloads, payload )
            
            if( isempty(allPayloads) )
                for ii = 1:length(obj.fnFMfile)
                     allPayloads.(obj.fnFMfile{ii}) = payload.(obj.fnFMfile{ii});
                end
            else
                for ii = 1:length(obj.fnFMfile)
                     allPayloads.(obj.fnFMfile{ii}) = [allPayloads.(obj.fnFMfile{ii}); payload.(obj.fnFMfile{ii})];
                end
            end
            
            [~, rec_filename] = fileparts(obj.recording_name);

            for ii = 1:length(obj.fnECGfile)
                 aux_payload.(obj.fnECGfile{ii}) = payload.(obj.fnECGfile{ii});
            end
            
            save([obj.tmp_path 'tmpfile_' rec_filename '_ECG_cantpids_' num2str(obj.cant_pids) '_' num2str(obj.this_pid) '_' num2str(obj.this_iteration) '_' num2str(obj.cant_iteration)  '.mat'], '-struct', 'aux_payload' );

            aux_ECG_filename = [obj.tmp_path 'tmpfile_' rec_filename '_ECG_cantpids_' num2str(obj.cant_pids) '_thispid_' num2str(obj.this_pid) '_iteration_' num2str(obj.this_iteration) '_of_' num2str(obj.cant_iteration)  '.mat'];
            %do this for cluster supercomputers reasons
            movefile(   [obj.tmp_path 'tmpfile_' rec_filename '_ECG_cantpids_' num2str(obj.cant_pids) '_' num2str(obj.this_pid) '_' num2str(obj.this_iteration) '_' num2str(obj.cant_iteration)  '.mat'], ...
                aux_ECG_filename, 'f' );
            
            cant_QRS = length(payload.QRS_locations);
            obj.QRS_idx = [obj.QRS_idx; uint16(colvec(1:cant_QRS))];
            obj.File_idx = [obj.File_idx; repmat(uint16(obj.file_count), cant_QRS, 1)];
            obj.File_list = [obj.File_list; cellstr(aux_ECG_filename)];
            obj.file_count = obj.file_count + 1;

        end
        
        
        
    end
   
end
   
function autovec = autovec_calculation(wtECGslice)

    mean_wtECGslice = mean(wtECGslice);    
    wtECGslice_cov = cov( bsxfun( @minus, wtECGslice, mean_wtECGslice ));
    [autovec autoval] = eig(wtECGslice_cov); 
    [~, autoval_idx] = sort(diag(autoval), 'descend');
    autovec = autovec(:,autoval_idx);
    
end

function [ ModMaxPos ZeroCross ] = CalcModMaxPos( wtSignal )

    ModMaxPos = nan;
    ZeroCross = nan;

    lwtSignal = size(wtSignal, 1);

    if( all( wtSignal == 0 ) )
        return
    end

    ac1 = conv(wtSignal(:,1), flipud(wtSignal(:,1)));

    ac1 = ac1(lwtSignal:end);
    ac1 = 1/ac1(1)*ac1;

    interp_idx = (1:0.1:lwtSignal)';

    ac1_interp = spline( 1:lwtSignal, ac1, interp_idx );
    iZeroCrossPos_ac1 = myzerocros(ac1);

    if( isempty(iZeroCrossPos_ac1))
        return
    else
        if( iZeroCrossPos_ac1(1) < lwtSignal )
            if( sign(ac1(iZeroCrossPos_ac1(1))) ~= sign(ac1(iZeroCrossPos_ac1(1)+1)) )
                ZeroCross = iZeroCrossPos_ac1(1) - ac1(iZeroCrossPos_ac1(1)) / ( ac1(iZeroCrossPos_ac1(1)+1) - ac1(iZeroCrossPos_ac1(1)) );
            else
                ZeroCross = (iZeroCrossPos_ac1(1)-1) - ac1(iZeroCrossPos_ac1(1)-1) / ( ac1(iZeroCrossPos_ac1(1)) - ac1(iZeroCrossPos_ac1(1)-1) );
            end
        else
            return
        end   
    end

    ModMaxPos_ac1 = modmax( ac1_interp, 2, 0, 0 );

    if( ~isempty(ModMaxPos_ac1) )
        valid_ac1_max_idx = find(interp_idx(ModMaxPos_ac1) > ZeroCross);
        if( ~isempty(valid_ac1_max_idx) )
            ModMaxPos = interp_idx(ModMaxPos_ac1(valid_ac1_max_idx(1)));
        else
            return
        end
    else
        return
    end

end

function ModMaxVal = CalcModMaxVal( wtSignal )

    % std_wtSignal = std(wtSignal);
    % wtSignal = bsxfun( @times, wtSignal, 1./std_wtSignal);

    cc = conv(wtSignal(:,1), flipud(wtSignal(:,2)));

    [~, max_pos] = max(abs(cc));

    ModMaxVal = cc(max_pos);

end