classdef ECGtask_classification_features_calc < ECGtask

% ECGtask for ECGwrapper (for Matlab)
% ---------------------------------
% 
% Description:
% 
% Abstract class for defining ECGtask interface
% 
% 
% Author: Mariano Llamedo Soria (llamedom at {electron.frba.utn.edu.ar; unizar.es}
% Version: 0.1 beta
% Birthdate  : 18/2/2013
% Last update: 18/2/2013
       
    properties(GetAccess = public, Constant)
        name = 'ECGclassificationTask';
        target_units = 'uV';
        doPayload = true;
        % if user = memory;
        % memory_constant is the fraction respect to user.MaxPossibleArrayBytes
        % which determines the maximum input data size.

    end
    
    properties(GetAccess = public, SetAccess = private)
        true_labels
        memory_constant = 0.15;
        started
    end
    
    properties( Access = 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
        q_filters
        LP
        LPsize
        scale_idx
        scales    
        series
        cant_QRS_locations 
    end
    
    properties
        progress_handle 
        recording_format
        class_labeling
        autovec
        tmp_path
                
    end
    
    methods
           
        function obj = ECGclassificationTask(obj)
            % not implemented
        end
        
        function Start(obj, ECG_header, ECG_annotations)
            
            % the dimension of the input block should be affected for the
            % internal sampling rate, since a resampling is performed
            % internally.
            obj.memory_constant = obj.memory_constant * ECG_header.freq / 360;
            
            %% 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;
            
            % Wavelet transform filter bank design
            obj.scales = 3:6;
            CantScales = length(obj.scales);
            MaxScales = max(obj.scales);
            obj.scale_idx= nan(MaxScales,1);

            filters_cache_filename = ['wt_filters_' num2str(obj.scales) ' scales_' num2str(obj.sampling_rate) ' Hz.mat' ];
            if( exist(filters_cache_filename, 'file') )
                aux = load( filters_cache_filename );
                obj.q_filters = aux.q_filters;
            else
                obj.q_filters = qs_filter_design(obj.scales, 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

            if( isfield( ECG_annotations, 'anntyp'))
                obj.true_labels = ECG_annotations.anntyp;
            end
            QRS_locations = ECG_annotations.time;
            
            % 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.
            obj.series.RR_intervals = round(RR_intervals * obj.sampling_rate / ECG_header.freq);
            obj.series.dRR_intervals = round(dRR_intervals * obj.sampling_rate / ECG_header.freq);
            obj.series.QRS_locations = QRS_locations;

            obj.cant_QRS_locations = length(obj.series.RR_intervals);
            
        end
        
        function payload = Process(obj, ECG, ECG_sample_start_end_idx, ECG_header, ECG_annotations, ECG_annotations_start_end_idx )
            
            payload = [];
            
            %refered to the whole ECG recording
            this_iter_QRS_idx = ECG_annotations_start_end_idx(1):ECG_annotations_start_end_idx(2);
            RR_intervals = obj.series.RR_intervals(this_iter_QRS_idx);
            
            %refered to this ECG strip
            QRS_locations = ECG_annotations.time;
            
            %% Preprocessing
            
            % Update point
            obj.progress_handle.checkpoint('Resampling');
            
            %resample to obj.sampling_rate Hz
            ECG = resample(ECG, obj.sampling_rate, ECG_header.freq);
            QRS_locations = colvec(round(QRS_locations * obj.sampling_rate / ECG_header.freq));
        %         this_iter_true_labels = true_labels(this_iter_QRS_idx);
            this_iter_ECG_resampled_size = size(ECG,1);

            if( ECG_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_handle.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_handle.checkpoint('Wavelet transform calculation');

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

            %% Features Calculation

            % Update point
            obj.progress_handle.checkpoint('Features calculation');
            
            this_iter_cant_QRS = length(QRS_locations);
            this_iter_seq_idx = 1:this_iter_cant_QRS;

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

            obj.progress_handle.checkpoint('Features calculation');

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

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

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

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

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

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

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

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

            aux_idx = arrayfun(@(a)( findStartEnd(  obj.series.RR_intervals >= (obj.series.RR_intervals(a) - 20*60*obj.sampling_rate) & ... 
                                                    obj.series.RR_intervals <= obj.series.RR_intervals(a) )), ...
                               this_iter_QRS_idx, 'UniformOutput', false);

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

            % Update point
            obj.progress_handle.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);

            %calculate PCA matrix in this slices for feature
            %AutoCorr1PlanoWTModMaxPos

            aux_idx2 = 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(b(1):b(2),:,obj.scale_idx(4))) * autovec_calculation(squeeze(wtECG(a(1):a(2),:,obj.scale_idx(4)))) ), aux_idx, aux_idx2, '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_handle.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_handle.checkpoint('Features calculation');

            %calculate clustering features.


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

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

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

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

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

            aux_idx = arrayfun(@(a)(max(1,a-1):min(obj.cant_QRS_locations,a+1)), this_iter_QRS_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} obj.cant_QRS_locations];
            end

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

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

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

            aux_idx = arrayfun(@(a)(max(1,a-1):min(obj.cant_QRS_locations,a+1)), this_iter_QRS_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} obj.cant_QRS_locations];
            end

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

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

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

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

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

            % Update point
            obj.progress_handle.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)( Calc_QRS_max_scale_proj( wtECG(a(1):a(2),:,:), obj.scales) ), aux_idx);

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

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

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

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

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


            % 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( ECG_header.nsig > 2 ) 
                % wtECG already projected in autovec.
                aux_featval = cellfun(@(a)( CalcModMaxVal( squeeze(wtECG(a(1):a(2),:,obj.scale_idx(3))) ) ), aux_idx);
            else
                aux_featval = cellfun(@(a)( CalcModMaxVal( squeeze(wtECG(a(1):a(2),:,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 payload = Finish(obj, payload, ECG_header)
            % not implemented
        end
        
        function plA = Concatenate(obj, plA, plB)
            
            if( isempty(plA) )
                for ii = 1:length(obj.fnFMfile)
                     plA.(obj.fnFMfile{ii}) = plB.(obj.fnFMfile{ii});
                end
            else
                for ii = 1:length(obj.fnFMfile)
                     plA.(obj.fnFMfile{ii}) = [plA.(obj.fnFMfile{ii}); plB.(obj.fnFMfile{ii})];
                end
            end
                        
        end

    end
    
end


%%%%%%%%%%%%%%%%%%
% Other functions
%%%%%%%%%%%%%%%%%%

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

function start_end_aux = findStartEnd( bAux )

    start_aux = find(  bAux, 1, 'first' );
    end_aux = find(  bAux, 1, 'last' );
    start_end_aux = [start_aux end_aux];
                                        
end

function aux_slice = Calc_QRS_max_scale_proj(aux_slice, scales)
    aux_slice = squeeze(sum(abs(aux_slice(:,1,:))))';
    aux_slice = scales * colvec(aux_slice) ./ sum(aux_slice);
end