%% cSPIKE v 1.0 17.8.2016: Eero Satuvuori
% This is the main part of the cSPIKE. The object offers all the
% functionality of graphical user interface SPIKY as a command line
% version. The main functions are implemented with MEX files and C++ 
% backends. Methods are based on ISI distance, SPIKE distance and SPIKE
% synchronization by Thomas Kreuz et al. The class contains data and
% methods for storing spiketrains as a single set as well as calculating
% the different distance measures. 

%% Class for enclosing spiketrains into one class and handling the set.
classdef SpikeTrainSet < handle
    
    
    properties (Access = private)
        
        % Internal datastores.
        DataArray;
        PooledArray;
        timeStart;
        timeEnd;
        threshold;
        timeShift;
        realISIs;
        
        % Internal accuracy used by the program. This is the maximum
        % accuracy used. Increasing the accuracy may lead to malfunctions.
        accuracy = 10^-12;
    end
    %% PRIVATE METHODS ****************************************************
    methods (Access = private)
        
        %% Checking data availability
        % Function is given a bool value if the information is crucial for
        % the calling function. In case it is an error is thrown. If it is 
        % not, just a bool value of 0 is returned if there is no data yet.
        function bool = HaveData(obj, crucial)
            
            % Returning 1 if there is data and 0 if there is not.
            if isempty(obj.DataArray)
                bool = 0;
                if crucial
                    error('No data available');
                end
            else
                bool = 1;
            end
            
        end
        %% Unique Union
        % Takes a set of intervals as a cell array so that each pair of two
        % values is an interval. The output is a unique union of the
        % intervals.
        function ReturnArray = UniqueUnion(obj, CellArray )
            starts = [];
            ends = [];
            for i = 1:2:length(CellArray)
                starts(end+1) = CellArray{i};
                ends(end+1) = CellArray{i+1};
            end
            
            UniqueArray = 0;
            
            while ~UniqueArray
                moves = 0;
                % Taking each first spike and if it is in an interval, move
                % it to the border of the interval.
                for i = 1:length(starts)
                    for ii = 1:length(starts)
                        
                        if (i ~= ii)
                            
                            if (starts(i) > starts(ii) && starts(i)< ends(ii))
                                
                                starts(i) = starts(ii);
                                moves = moves + 1;
                            end
                        end
                        
                    end
                end
                % Taking each last spike and doing the same.
                for i = 1:length(ends)
                    for ii = 1:length(ends)
                        
                        if (i ~= ii)
                            
                            if (ends(i) > starts(ii) && ends(i) < ends(ii))
                                
                                ends(i) = ends(ii);
                                moves = moves + 1;
                            end
                        end
                        
                    end
                end
                
                % If there were no moves this round, all the spikes are at
                % the edges of their intervals. There may be multiple same
                % intervals.
                if moves  == 0
                    
                    % now the array is unique. No need to continue.
                    UniqueArray = 1;
                    
                    % Taking each interval only once.
                    Ustarts = unique(starts);
                    Uends = unique(ends);
                    
                    % Since the intervals are in time order and unique,
                    % they will be in correct time order after unique
                    % command. The output is formed based on the intervals
                    ReturnArray = cell(0);
                    for i = 1:length(Ustarts)
                        
                        ReturnArray{end+1} = Ustarts(i);
                        ReturnArray{end+1} = Uends(i);
                    end
                end
            end
        end
        
        %% Checking that the DATA given is valid
        % Function checks the data and gives an error
        % if the data is not in a correct form
        function CheckData(obj,Data)
            
            if ~iscell(Data)
                error('Data is not a cell array');
            end
            if length(Data) == 1
                error('One spiketrain');
            end
            for i = 1:length(Data)
                
                if ndims(Data)~=2 || size(Data,2) ~= 1;
                    error('Data has incorrectly aligned cell array');
                end
                
                if ndims(Data{i}) ~=2 && (size(Data{i},1)||size(Data{i},2))
                    error('Data has too many dimensions');
                end
                
                if size(Data{i},1) ~= 1
                    if  ~isempty(Data{i})
                        error('Data has incorrectly aligned vectors');
                    end
                end
                
                if ~isnumeric(Data{i})
                    error('Data is not all numeric');
                end
            end
            
            
        end
        %% Cutting function
        % This function cuts the data so that spikes outside of boundaries
        % are removed.
        function Data = CutData(obj,Data, START, END)
            
            
            for i = 1:length(Data)
                inds = not(abs(sign(sign(START-obj.accuracy - Data{i})...
                       + sign(END+obj.accuracy - Data{i}))));
                Data{i} = Data{i}(inds);
            end
            
        end
        %%
        function Data = RoundData(obj, Data)
            
            dt = obj.accuracy;
            for i = 1:length(Data)
                Data{i} = round(Data{i}/dt)*dt;
            end
            
        end
        %% Check input times
        % Checking that times given are in correct order and within 
        % boundaries
        function CheckTime(obj,Time1,Time2)
            if Time1  > Time2
                error('Time interval boundaries incorrect');
            end
            
            if ~((Time1 >= obj.timeStart) && (Time2 <= obj.timeEnd))
                error('time interval not within data bounds');
            end
            
        end
        
        
        %% Correct the DATA Edges
        % Adds edge correction and makes sure that there are only 1 spike
        % at each time instant. Also rearranges the data into time series.
        
        function Corrected = EdgeCorrection(obj,Data,time1,time2)

            PooledISIs = [];
            Corrected = cell(1,length(Data))';
            for i = 1:length(Data)
                % 0 or 1 spikes cases
                if isempty(Data{i}) || length(Data{i}) == 1
                    if isempty(Data{i})
                        first = -(time2-time1)*1.2;
                        last = (time2-time1)*2.2;
                        Corrected{i} = [first obj.timeStart obj.timeEnd last];
                        Corrected{i} = unique(Corrected{i});
                        thisST = obj.timeEnd-obj.timeStart;
                    else
                        fullST = 0; %flag
                        first = obj.timeStart;
                        last = obj.timeEnd;
                        if Data{i} == first
                            first = -(time2-time1)*1.2;
                            thisST = obj.timeEnd-obj.timeStart;
                            fullST = 1;
                        end
                        if Data{i} == last
                            last = (time2-time1)*2.2;
                            thisST = obj.timeEnd-obj.timeStart;
                            fullST = 1;
                        end
                        
                        if ~fullST
                            thisST = [last-Data{i} Data{i}-first];
                        end
                        
                        Corrected{i} = [first Data{i} last];
                        Corrected{i} = unique(Corrected{i});
                    end
                    
                else
                    %The rest
                    thisST = Data{i}(2:end) - Data{i}(1:end-1);
                    if(Data{i}(1) ~= time1)
                        first = Data{i}(1) - max( Data{i}(2)-Data{i}(1), Data{i}(1) - obj.timeStart );
                        thisST = [thisST Data{i}(1)-first];
                    else
                        % In case the first spike is at the border it is
                        % not treated as auxiliary spike. For calculation
                        % the next spike does not matter so we add an
                        % "auxiliary" spike very far from the dataset to
                        % mark that the border is a real one.
                        first = -(time2-time1)*1.2;
                    end
                    if(Data{i}(end)~= time2)
                        last = Data{i}(end) + max( Data{i}(end)-Data{i}(end-1),obj.timeEnd - Data{i}(end));
                        thisST = [thisST last-Data{i}(end)];
                    else
                        % Same for the other end.
                        last = (time2-time1)*2.2;
                    end
                    
                    
                    
                    Corrected{i} = [first last Data{i}];
                    
                    % At the end making sure that each dataset has only one
                    % spike at each time instant and that they are in time 
                    % order
                    Corrected{i} = unique(Corrected{i});
                end
                PooledISIs = [PooledISIs thisST];
            end
            
            if(isempty (PooledISIs))
                obj.threshold = 0;
            else
                obj.threshold = sqrt(mean(PooledISIs.^2));
            end
            obj.realISIs = PooledISIs;
            
        end
        
        %% Form pooled spiketrain of all spiketrains
        % Pools all spiketrains into one. Uses unique, so that multiple
        % spikes at same times are filtered out.
        function Pooled = Pool(obj , Data)
            
            Pooled = [];
            for i = 1:length(Data)
                Pooled = [Pooled Data{i}];
            end
            Pooled = unique(Pooled);
            
        end
        %% Time shift functions.
        % These functions take care of the difference between "real" times
        % and internal times. The object computes using positive spiketimes
        % and in case the recording begins at negative times, all data is
        % shifted so that the recording begins from 0 internally. For
        % output the times are again shifted back to original.
        function  [time1,time2,Data] = CreateTimeShift(obj,time1,time2,Data)
            
            obj.timeShift = time1;
            time1 = time1-obj.timeShift;
            time2 = time2-obj.timeShift;
            for i = 1:length(Data)
                Data{i}=Data{i}-obj.timeShift;
            end
            
        end
        
        function  [time] = InputTimeShift(obj,time)
            
            time = time-obj.timeShift;
            
        end
        
        
        function  [time] = outputTimeShift(obj,time)
            
            time = time+obj.timeShift;
            
        end
        
        
    end
    %% PUBLIC METHODS *****************************************************
    methods (Access = public)
        %% Constructor
        % Builds a spiketrain set with the spiketrains given and beginning
        % and end times.
        function obj = SpikeTrainSet(spiketrains, beginning, ending )
            obj.timeShift = 0;
            obj.SetData(spiketrains, beginning, ending );
            
        end
        
        %% Setting data
        % Function sets data to the classes own variables. Calls
        % appropriate checking functions to input. Data is only accepted
        % in a single cell array containing spiketrains. If called to a
        % spiketrain set that has already been set, the user is asked for
        % input if it is to be replaced.
        function SetData(obj, Data, time1, time2)
            Data = Data';
            % checking that the data is in correct form.
            obj.CheckData(Data);
            %Rounding data to 12 digits accuracy
            Data = obj.RoundData(Data);
            % Cutting spikes outside of interval
            Data = obj.CutData(Data, time1, time2);
            % If there was no data, add the data
            if obj.HaveData(0)
                
                %presetting answer variable into not n or y
                answer = 'h';
                % If there was data, ask if the data is to be replaced
                while answer ~= 'y' && answer ~= 'Y' && answer ~= 'n' && answer ~= 'N'
                    answer = input('Do you want to replace the data?(y/n)','s');
                end
                %
                if answer == 'y' || answer ~= 'Y'
                    
                    % Swiping old arrays clean
                    obj.DataArray = [];
                    obj.PooledArray= [];
                    
                    % Adding new values to the array
                    obj.DataArray{2} = Data;
                    [time1,time2,Data] = obj.CreateTimeShift(time1,time2,Data);
                    obj.timeStart = time1;
                    obj.timeEnd = time2;
                    obj.DataArray{1} = obj.EdgeCorrection(Data);
                    obj.PooledArray = obj.Pool(Data);
                    
                end
            else
                % Setting the arrays for the first time
                obj.DataArray{2} = Data;
                [time1,time2,Data] = obj.CreateTimeShift(time1,time2,Data);
                obj.timeStart = time1;
                obj.timeEnd = time2;
                obj.DataArray{1} = obj.EdgeCorrection(Data,time1,time2);
                obj.PooledArray = obj.Pool(Data);
                
            end
            
        end
        
        %% ISI METHODS ----------------------------------------------------
        % The methods in this section are applying ISI distance measure to
        % the data. The section also contains Adaptive extension.
        
        %% ISI-distance for interval/point
        % Calculates ISI distance for given interval between time1 and
        % time2. time1 < time2. If no input times are given return value
        % for whole data set interval that is defined. If only one time is
        % given or both are the same, return dissimilarity of the
        % spiketrains at that time instant.
        function  ISId = ISIdistance(obj, time1, time2)
            if nargin == 1
                time1 = outputTimeShift(obj,obj.timeStart);
                time2 = outputTimeShift(obj,obj.timeEnd);
            end
            if nargin == 2
                time2 = time1;
            end
            
            ISId = AdaptiveISIdistance(obj, time1, time2, 0);
            
        end
        %% Adaptive ISI-distance for interval/point
        % Calculates SPIKE distance for given interval between time1 and
        % time2. time1 < time2. If no input times are given return value
        % for whole data set interval that is defined. If only one time is
        % given or both are the same, return dissimilarity of the
        % spiketrains at that time instant.
        function  ISId = AdaptiveISIdistance(obj, time1, time2, threshold)
            if nargin == 1
                time1 = obj.timeStart;
                time2 = obj.timeEnd;
            else
                time1 = obj.InputTimeShift(time1);
                time2 = obj.InputTimeShift(time2);
            end
            if nargin == 2
                time2 = time1;
            end
            if nargin < 4
                threshold = obj.threshold;
            end
            obj.HaveData(1);
            obj.CheckTime(time1,time2);
            if(time1 == time2)
                ISId = mexAdaptiveISIDistance(obj.DataArray{1},threshold,time1);
            else
                ISId = mexAdaptiveISIDistance(obj.DataArray{1},threshold,time1, time2);
            end
        end
        %% ISI-distance/dissimilarity matrix
        % Returns the distance matrix of the spike trains averaged over
        % time interval defined by time 1 and time 2. If not given,
        % averaged over whole interval. If given only one time or the same
        % one twice, the returned matrix is dissimilarity between all
        % pairs.
        function  ISIdM = ISIdistanceMatrix(obj, time1, time2)
            if nargin == 1
                time1 = outputTimeShift(obj,obj.timeStart);
                time2 = outputTimeShift(obj,obj.timeEnd);
            end
            if nargin == 2
                time2 = time1;
            end
            
            ISIdM = AdaptiveISIdistanceMatrix(obj, time1, time2, 0);
        end
        %% Adaptive ISI-distance/dissimilarity matrix
        % Returns the adaptive distance matrix of the spike trains averaged 
        % over time interval defined by time 1 and time 2. If not given,
        % averaged over whole interval. If given only one time or the same
        % one twice, the returned matrix is dissimilarity between all
        % pairs. Threshold can be given as a parameter. If not given the
        % one extracted from the data is used.
        function  ISIdM = AdaptiveISIdistanceMatrix(obj, time1, time2,threshold)
            
            if nargin == 1
                time1 = obj.timeStart;
                time2 = obj.timeEnd;
            else
                time1 = obj.InputTimeShift(time1);
                time2 = obj.InputTimeShift(time2);
            end
            
            if nargin == 2
                time2 = time1;
            end
            if nargin < 4
                threshold = obj.threshold;
            end
            obj.HaveData(1);
            obj.CheckTime(time1,time2);
            
            ISIdM = mexAdaptiveISIDistanceMatrix(obj.DataArray{1},threshold ,time1, time2);
            
        end
        %% Triggered ISI-dissimilarity matrix
        % Returns the averaged dissimilarity matrix of the spike trains
        % averaged over time instants defined in inputs. Each instant is 
        % given as new variables: (1,2,5,6) being average over time 
        % instants 1,2,5 and 6.
        function  ISIdM = TriggeredISImatrix(obj, varargin)
            
            if nargin < 3
                error('Not enough input arguments');
            end
            ISIdM = TriggeredAdaptiveISImatrix(obj,0, varargin);
        end
        
        %% Triggered Adaptive ISI-dissimilarity matrix
        % Returns the averaged dissimilarity matrix of the spike trains
        % averaged over time instants defined in inputs. Each instant is 
        % given as new variables: (1,2,5,6) being average over time 
        % instants 1,2,5 and 6. If a negative threshold value is given the 
        % data exctracted threshold is used.
        function  ISIdM = TriggeredAdaptiveISImatrix(obj,threshold, varargin)
            if threshold < 0
                threshold = obj.threshold;
            end
            if nargin < 4
                error('Not enough input arguments');
            end
            
            obj.HaveData(1);
            time = obj.InputTimeShift(varargin{1});
            obj.CheckTime(time,time);
            
            ISIdM = mexAdaptiveISIDistanceMatrix(obj.DataArray{1},threshold,time, time);
            for i = 2:(nargin-1)
                time = obj.InputTimeShift(varargin{i});
                obj.CheckTime(time,time);
                ISIdM = ISIdM + mexAdaptiveISIDistanceMatrix(obj.DataArray{1},threshold,time, time);
            end
            
            ISIdM = ISIdM./(nargin-2);
        end
        
        %% Averaged ISI-distance matrix
        % Returns the distance matrix of the spike trains averaged over
        % time intervals defined as pairs of inputs. The inputs must be in
        % pairs. If the intervals overlap, the average is computed over
        % unique union of the intervals. Input is given as new pairs of 
        % variables: (1,2,5,6) being average over time intervals [1,2] and 
        % [5,6].
        function  ISIdM = AveragedISIdistanceMatrix(obj, varargin)
            if nargin < 3
                error('Not enough input arguments');
            end
            if mod(nargin,2) == 0;
                error('Not a paired number of time points');
            end
            
            ISIdM = AveragedAdaptiveISIdistanceMatrix(obj,0, varargin);
            
        end
        
        %% Averaged Adaptive ISI-distance matrix
        % Returns the distance matrix of the spike trains averaged over
        % time intervals defined as pairs of inputs. The inputs must be in
        % pairs. If the intervals overlap, the average is computed over
        % unique union of the intervals. Input is given as new pairs of 
        % variables: (1,2,5,6) being average over time intervals [1,2] and 
        % [5,6]. Threshold can be given as a parameter. If a negative
        % threshold value is given the data exctracted threshold is used.
        function  ISIdM = AveragedAdaptiveISIdistanceMatrix(obj,threshold, varargin)
            if threshold < 0
                threshold = obj.threshold;
            end
            if nargin < 4
                error('Not enough input arguments');
            end
            if mod(nargin,2) == 1;
                error('Not a paired number of time points');
            end
            
            obj.HaveData(1);
            
            % Forming unique union of input
            Timeslices = obj.UniqueUnion(varargin);
            
            time1 = Timeslices{1};
            time2 = Timeslices{2};
            time1 = obj.InputTimeShift(time1);
            time2 = obj.InputTimeShift(time2);
            obj.CheckTime(time1,time2);
            TotalTime = time2-time1;
            ISIdM = (time2-time1)*mexISIistanceMatrix(obj.DataArray{1},threshold,time1, time2);
            
            for i = 3:2:length(Timeslices)
                
                time1 = Timeslices{i};
                time2 = Timeslices{i+1};
                time1 = obj.InputTimeShift(time1);
                time2 = obj.InputTimeShift(time2);
                obj.CheckTime(time1,time2);
                ISIdM = ISIdM + (time2-time1)*mexISIistanceMatrix(obj.DataArray{1},threshold,time1, time2);
                TotalTime = TotalTime + (time2-time1);
            end
            
            % averaging over time. Each matrix is weighted
            ISIdM = ISIdM./TotalTime;
        end
        
        %% ISI-distance profile
        % Forms a profile between times 1 and 2. The output is a Matlab
        % class Profile, which contains the profile and can even perform a
        % simple plotting. See object Profile for more information.
        function  ISIp = ISIdistanceProfile(obj, time1, time2)
            if nargin == 1
                time1 = outputTimeShift(obj,obj.timeStart);
                time2 = outputTimeShift(obj,obj.timeEnd);
            end
            if nargin == 2
                time2 = time1;
            end
            
            ISIp = AdaptiveISIdistanceProfile(obj, time1, time2,0);
        end
        %% Adaptive ISI-distance profile
        % Forms a profile between times 1 and 2. The output is a Matlab
        % class Profile, which contains the profile and can even perform a
        % simple plotting. See object Profile for more information.
        function  ISIp = AdaptiveISIdistanceProfile(obj, time1, time2,threshold)
            
            if nargin == 1
                time1 = obj.timeStart;
                time2 = obj.timeEnd;
            else
                time1 = obj.InputTimeShift(time1);
                time2 = obj.InputTimeShift(time2);
            end
            if nargin < 4
                threshold = obj.threshold;
            end
            obj.HaveData(1);
            obj.CheckTime(time1,time2);
            
            internalProfile = mexAdaptiveISIDistanceProfile(obj.DataArray{1},threshold,time1, time2);
            internalProfile(2,:) = obj.outputTimeShift(internalProfile(2,:));
            ISIp = Profile(internalProfile);
            
        end
        %% SPIKE METHODS --------------------------------------------------
        % The methods in this section are applying the SPIKE-distance 
        % measure to the data. The section also contains Adaptive and rate 
        % intedependent extensions. 
        
        %% SPIKE-distance for interval/point
        % Calculates SPIKE distance for given interval between time1 and
        % time2. time1 < time2. If no input times are given return value
        % for whole data set interval that is defined. If only one time is
        % given or both are the same, return dissimilarity of the
        % spiketrains at that time instant.
        function  SPIKEd = SPIKEdistance(obj, time1, time2)
            if nargin == 1
                time1 = outputTimeShift(obj,obj.timeStart);
                time2 = outputTimeShift(obj,obj.timeEnd);
            end
            if nargin == 2
                time2 = time1;
            end
            
            SPIKEd = AdaptiveSPIKEdistance(obj, time1, time2, 0);
        end
        %% Adaptive SPIKE-distance for interval/point
        % Calculates SPIKE distance for given interval between time1 and
        % time2. time1 < time2. If no input times are given return value
        % for whole data set interval that is defined. If only one time is
        % given or both are the same, return dissimilarity of the
        % spiketrains at that time instant.
        function  SPIKEd = AdaptiveSPIKEdistance(obj, time1, time2, threshold)
            
            if nargin == 1
                time1 = obj.timeStart;
                time2 = obj.timeEnd;
            else
                time1 = obj.InputTimeShift(time1);
                time2 = obj.InputTimeShift(time2);
            end
            if nargin == 2
                time2 = time1;
            end
            if nargin < 4
                threshold = obj.threshold;
            end
            
            obj.HaveData(1);
            obj.CheckTime(time1,time2);
            
            if(time1 == time2)
                SPIKEd = mexAdaptiveSPIKEDistance(obj.DataArray{1},threshold,time1);
            else
                SPIKEd = mexAdaptiveSPIKEDistance(obj.DataArray{1},threshold,time1, time2);
            end
        end
        
        %% Adaptive Rate independent SPIKE distance for interval/point
        % Calculates SPIKE distance for given interval between time1 and
        % time2. time1 < time2. If no input times are given return value
        % for whole data set interval that is defined. If only one time is
        % given or both are the same, return dissimilarity of the
        % spiketrains at that time instant. For non-Adaptive version call
        % the function without version or threshold parameters. For calling
        % the adaptive version, set version parameter to 2. This will use
        % the predetermined threshold. If you then want to set threshold
        % yourself give it as fifth parameter.
        function  SPIKEd = AdaptiveRateIndependentSPIKEdistance(obj, time1, time2,version, threshold)
            
            if nargin == 1
                time1 = obj.timeStart;
                time2 = obj.timeEnd;
            else
                time1 = obj.InputTimeShift(time1);
                time2 = obj.InputTimeShift(time2);
            end
            if nargin < 4
                threshold = 0;
            end
            if nargin == 4
                if version == 2
                    threshold = obj.threshold;
                else
                    threshold = 0;
                end
            end
            
            obj.HaveData(1);
            obj.CheckTime(time1,time2);
            
            if(time1 == time2)
                SPIKEd = mexAdaptiveRateIndependentSPIKEDistance(obj.DataArray{1},threshold,time1);
            else
                SPIKEd = mexAdaptiveRateIndependentSPIKEDistance(obj.DataArray{1},threshold,time1, time2);
            end
        end
        %% SPIKE distance/dissimilarity matrix
        % Returns the distance matrix of the spike trains averaged over
        % time interval defined by time 1 and time 2. If not given,
        % averaged over whole interval. If given only one time or the same
        % one twice, the returned matrix is dissimilarity between all
        % pairs.
        function  SPIKEdM = SPIKEdistanceMatrix(obj, time1, time2)
            if nargin == 1
                time1 = outputTimeShift(obj,obj.timeStart);
                time2 = outputTimeShift(obj,obj.timeEnd);
            end
            if nargin == 2
                time2 = time1;
            end
            
            SPIKEdM = AdaptiveSPIKEdistanceMatrix(obj, time1, time2, 0);
        end
        
        %% Adaptive SPIKE distance/dissimilarity matrix
        % Returns the distance matrix of the spike trains averaged over
        % time interval defined by time 1 and time 2. If not given,
        % averaged over whole interval. If given only one time or the same
        % one twice, the returned matrix is dissimilarity between all
        % pairs.
        function  SPIKEdM = AdaptiveSPIKEdistanceMatrix(obj, time1, time2,threshold)
            
            if nargin == 1
                time1 = obj.timeStart;
                time2 = obj.timeEnd;
            else
                time1 = obj.InputTimeShift(time1);
                time2 = obj.InputTimeShift(time2);
            end
            if nargin == 2
                time2 = time1;
            end
            if nargin < 4
                threshold = obj.threshold;
            end
            obj.HaveData(1);
            obj.CheckTime(time1,time2);
            
            SPIKEdM = mexAdaptiveSPIKEDistanceMatrix(obj.DataArray{1},threshold,time1, time2);
            
        end
        
        %% Triggered SPIKE-dissimilarity matrix
        % Returns the averaged dissimilarity matrix of the spike trains
        % averaged over time instants defined in inputs. Each instant is 
        % given as new variables: (1,2,5,6) being average over time 
        % instants 1,2,5 and 6.
        function  SPIKEdM = TriggeredSPIKEmatrix(obj, varargin)
            
            if nargin < 3
                error('Not enough input arguments');
            end
            SPIKEdM = TriggeredAdaptiveSPIKEmatrix(obj,0, varargin);
            
        end
        
        %% Triggered Adaptive SPIKE-dissimilarity matrix
        % Returns the averaged dissimilarity matrix of the spike trains
        % averaged over time instants defined in inputs. Each instant is 
        % given as new variables: (1,2,5,6) being average over time 
        % instants 1,2,5 and 6. If a negative threshold value is given the 
        % data exctracted threshold is used.
        
        function  SPIKEdM = TriggeredAdaptiveSPIKEmatrix(obj,threshold, varargin)
            if threshold < 0
                threshold = obj.threshold;
            end
            if nargin < 4
                error('Not enough input arguments');
            end
            
            obj.HaveData(1);
            time = obj.InputTimeShift(varargin{1});
            obj.CheckTime(time,time);
            
            SPIKEdM = mexAdaptiveSPIKEDistanceMatrix(obj.DataArray{1},threshold,time, time);
            for i = 2:(nargin-1)
                time = obj.InputTimeShift(varargin{i});
                obj.CheckTime(time,time);
                SPIKEdM = SPIKEdM + mexAdaptiveSPIKEDistanceMatrix(obj.DataArray{1},threshold,time, time);
            end
            
            SPIKEdM = SPIKEdM./(nargin-1);
        end
        
        %% Averaged SPIKE-distance matrix
        % Returns the distance matrix of the spike trains averaged over
        % time intervals defined as pairs of inputs. The inputs must be in
        % pairs. If the intervals overlap, the average is computed over
        % unique union of the intervals. Input is given as new pairs of 
        % variables: (1,2,5,6) being average over time intervals [1,2] and 
        % [5,6].
        function  SPIKEdM = AveragedSPIKEdistanceMatrix(obj, varargin)
            
            if nargin < 3
                error('Not enough input arguments');
            end
            if mod(nargin,2) == 0;
                error('Not a paired number of time points');
            end
            SPIKEdM = AveragedAdaptiveSPIKEdistanceMatrix(obj, 0, varargin);
        end
        
        %% Averaged Adaptive SPIKE-distance matrix
        % Returns the distance matrix of the spike trains averaged over
        % time intervals defined as pairs of inputs. The inputs must be in
        % pairs. If the intervals overlap, the average is computed over
        % unique union of the intervals. Input is given as new pairs of 
        % variables: (1,2,5,6) being average over time intervals [1,2] and 
        % [5,6]. Threshold can be given as a parameter. If a negative
        % threshold value is given the data exctracted threshold is used.
        
        
        function  SPIKEdM = AveragedAdaptiveSPIKEdistanceMatrix(obj, threshold, varargin)
            if threshold < 0
                threshold = obj.threshold;
            end
            if nargin < 4
                error('Not enough input arguments');
            end
            if mod(nargin,2) == 1;
                error('Not a paired number of time points');
            end
            
            obj.HaveData(1);
            
            % Forming unique union of input
            Timeslices = obj.UniqueUnion(varargin);
            
            time1 = Timeslices{1};
            time2 = Timeslices{2};
            time1 = obj.InputTimeShift(time1);
            time2 = obj.InputTimeShift(time2);
            obj.CheckTime(time1,time2);
            TotalTime = time2-time1;
            SPIKEdM = (time2-time1)*mexAdaptiveSPIKEDistanceMatrix(obj.DataArray{1},threshold,time1, time2);
            
            for i = 3:2:length(Timeslices)
                
                time1 = Timeslices{i};
                time2 = Timeslices{i+1};
                time1 = obj.InputTimeShift(time1);
                time2 = obj.InputTimeShift(time2);
                obj.CheckTime(time1,time2);
                SPIKEdM = SPIKEdM + (time2-time1)*mexAdaptiveSPIKEDistanceMatrix(obj.DataArray{1},threshold,time1, time2);
                TotalTime = TotalTime + (time2-time1);
            end
            
            % averaging over time. Each matrix is weighted
            SPIKEdM = SPIKEdM./TotalTime;
        end
        
        %% SPIKE-distance profile
        % Forms a profile between times 1 and 2. The output is a Matlab
        % class Profile, which contains the profile and can even perform a
        % simple plotting.
        function  SPIKEp = SPIKEdistanceProfile(obj, time1, time2)
            if nargin == 1
                time1 = outputTimeShift(obj,obj.timeStart);
                time2 = outputTimeShift(obj,obj.timeEnd);
            end
            if nargin == 2
                time2 = time1;
            end
            
            SPIKEp = AdaptiveSPIKEdistanceProfile(obj, time1, time2, 0);
            
        end
        %% Adaptive SPIKE-distance profile
        % Forms a profile between times 1 and 2. The output is a Matlab
        % class Profile, which contains the profile and can even perform a
        % simple plotting.
        function  SPIKEp = AdaptiveSPIKEdistanceProfile(obj, time1, time2,threshold)
            
            if nargin == 1
                time1 = obj.timeStart;
                time2 = obj.timeEnd;
            else
                time1 = obj.InputTimeShift(time1);
                time2 = obj.InputTimeShift(time2);
            end
            if nargin < 4
                threshold = obj.threshold;
            end
            obj.HaveData(1);
            obj.CheckTime(time1,time2);
            
            internalProfile = mexAdaptiveSPIKEDistanceProfile(obj.DataArray{1},threshold,time1, time2);
            internalProfile(2,:) = obj.outputTimeShift(internalProfile(2,:));
            SPIKEp = Profile(internalProfile);
            
        end
        
        
        %% SPIKE-Synchronization METHODS ----------------------------------
        % The methods in this section are applying SPIKE-synchronization
        % measure to the data
        
        %% SPIKE-Synchronization 
        % This function returns the value of the method from the given
        % interval given in the input.
        function  synchro = SPIKEsynchro(obj,time1,time2)
            if nargin == 1
                time1 = outputTimeShift(obj,obj.timeStart);
                time2 = outputTimeShift(obj,obj.timeEnd);
            end
            if nargin == 2
                time2 = time1;
            end
            
            synchro = AdaptiveSPIKEsynchro(obj,time1,time2, 0);
        end
        
        %% Adaptive SPIKE-synchronization 
        % This function returns the value of the method from the given
        % interval given in the input. If no threshold value is given the
        % data exctracted value is used.
        function  synchro = AdaptiveSPIKEsynchro(obj,time1,time2,threshold)
            obj.HaveData(1);
            if nargin == 1
                time1 = obj.timeStart;
                time2 = obj.timeEnd;
            else
                time1 = obj.InputTimeShift(time1);
                time2 = obj.InputTimeShift(time2);
            end
            if nargin < 4
                threshold = obj.threshold;
            end
            obj.HaveData(1);
            obj.CheckTime(time1,time2);
            
            synchro = mexAdaptiveSPIKESynchro(obj.DataArray{1},threshold/2,obj.timeStart,obj.timeEnd,time1,time2);
            
        end
        %% SPIKE-synchronization profile
        % Forms a profile between times 1 and 2. The output is a Matlab
        % class Profile, which contains the profile and can even perform a
        % simple plotting.
        function  [synchro, Sorder, STOrder] = SPIKEsynchroProfile(obj, time1, time2)
            if nargin == 1
                time1 = outputTimeShift(obj,obj.timeStart);
                time2 = outputTimeShift(obj,obj.timeEnd);
            end
            if nargin == 2
                time2 = time1;
            end
            
            [synchro, Sorder, STOrder] = AdaptiveSPIKEsynchroProfile(obj, time1, time2, 0);
            
        end
        %% Adaptive SPIKE-synchronization profile
        % Forms a profile between times 1 and 2. The output is a Matlab
        % class Profile, which contains the profile and can even perform a
        % simple plotting.
        function  [synchro, Sorder, STOrder] = AdaptiveSPIKEsynchroProfile(obj, time1, time2,threshold)
            
            if nargin == 1
                time1 = obj.timeStart;
                time2 = obj.timeEnd;
            else
                time1 = obj.InputTimeShift(time1);
                time2 = obj.InputTimeShift(time2);
            end
            if nargin < 4
                threshold = obj.threshold;
            end
            
            obj.HaveData(1);
            obj.CheckTime(time1,time2);
            % Change later into a real profile
            [synchro, Sorder, STOrder] = mexAdaptiveSPIKESynchroProfile(obj.DataArray{1},threshold/2,obj.timeStart,obj.timeEnd,time1, time2);
            
        end
        
        %% SPIKE-synchronization matrix
        % Returns the distance matrix of the spike trains averaged over
        % time interval defined by time 1 and time 2. If not given,
        % averaged over whole interval. If given only one time or the same
        % one twice, the returned matrix is dissimilarity between all
        % pairs.
        
        function  [SPIKESM,SPIKEOM,normSPIKEOM] = SPIKESynchroMatrix(obj, time1, time2)
            
            if nargin == 1
                time1 = outputTimeShift(obj,obj.timeStart);
                time2 = outputTimeShift(obj,obj.timeEnd);
            end
            if nargin == 2
                time2 = time1;
            end
            
            [SPIKESM,SPIKEOM,normSPIKEOM] = AdaptiveSPIKESynchroMatrix(obj, time1, time2, 0);
            
        end
        
        %% Adaptive SPIKE-synchronization matrix
        % Returns the distance matrix of the spike trains averaged over
        % time interval defined by time 1 and time 2. If not given,
        % averaged over whole interval. If given only one time or the same
        % one twice, the returned matrix is dissimilarity between all
        % pairs. If no threshold value is given thedata exctracted value 
        % is used.
        function  [SPIKESM,SPIKEOM,normSPIKEOM] = AdaptiveSPIKESynchroMatrix(obj, time1, time2, threshold)
            
            if nargin == 1
                time1 = obj.timeStart;
                time2 = obj.timeEnd;
            else
                time1 = obj.InputTimeShift(time1);
                time2 = obj.InputTimeShift(time2);
            end
            if nargin == 2
                time2 = time1;
            end
            if nargin < 4
                threshold = obj.threshold;
            end
            obj.HaveData(1);
            obj.CheckTime(time1,time2);
            
            NroSpiketrains = length(obj.DataArray{1});
            SPIKESM = ones(NroSpiketrains);
            SPIKEOM = zeros(NroSpiketrains);
            normSPIKEOM = zeros(NroSpiketrains);
            for i = 1:NroSpiketrains-1
                for ii = i:NroSpiketrains
                    if i ~= ii
                        pair{1} = obj.DataArray{2}{i};
                        pair{2} = obj.DataArray{2}{ii};
                        spikes = length(obj.DataArray{1}{i})+length(obj.DataArray{1}{ii})-4;
                        if spikes > 0
                            STS = SpikeTrainSet(pair,obj.timeStart,obj.timeEnd);
                            SPIKESM(i,ii) = STS.AdaptiveSPIKEsynchro(time1, time2, threshold);
                            SPIKESM(ii,i) = SPIKESM(i,ii);
                            
                            [~,~,profile] = STS.AdaptiveSPIKEsynchroProfile(time1, time2, threshold);
                            SPIKEOM(i,ii) = sum(profile{1});
                            SPIKEOM(ii,i) = sum(profile{2});
                            
                            normSPIKEOM(i,ii) = SPIKEOM(i,ii)/spikes;
                            normSPIKEOM(ii,i) = SPIKEOM(ii,i)/spikes;
                        end
                    end
                end
            end
            
        end
        
        %% Auxiliary functions --------------------------------------------
        
        %% Get function
        % Gives a copy of the spike data contained in the spike train set
        function  DATA = giveDATA(obj)
            obj.HaveData(1);
            DATA = obj.DataArray;
        end
        %% Threshold value extracted from the data
        % returns the threshold value extracted from this dataset
        function  [THR, Mean] = giveTHR(obj)
            obj.HaveData(1);
            THR = obj.threshold;
            Mean = mean(obj.realISIs);
        end

        %% Start and end of the data set
        % returns an array with start and end times of the recording.
        function  [time1,time2] = giveTIMES(obj)
            obj.HaveData(1);
            time1 = obj.outputTimeShift(obj.timeStart);
            time2 = obj.outputTimeShift(obj.timeEnd);
        end
        
        %% Plotting function
        % Plots the spike trains to current axis
        function  plotSpikeTrainSet(obj,color,widith)
            obj.HaveData(1);
            if nargin == 1
                color = 'k';
            end
            if nargin < 3
                widith = 1;
            end
            
            beginning = obj.outputTimeShift(obj.timeStart);
            ending = obj.outputTimeShift(obj.timeEnd);
            hold on;
            for i = 1:length(obj.DataArray{2})
                
                SpikeTrain = obj.DataArray{2}{i};
                plotSpikeTrain(SpikeTrain,0.8,(length(obj.DataArray{2})+1-i)-0.5,color,widith)
                
            end
            axis([beginning,ending,0,length(obj.DataArray{2})])
            set(gca,'ytick',[]);
        end
        
        
    end
  
end

