try: 
    import psyco
    psyco.full()
except:
    print "Impossible to import psyco library"
    pass

from importlibs import *



__author__  = "Andrea Maffezzoli"
__revision__ = "$Revision: 5.56 - 20 Jun 2011"
__name__    = "burstdetection"
__all__     = ["BurstDetection_root", "maxinterval", "netbursts"]



class BurstDetection_root(object):
    """
    Root class of maxinterval.
    General informations about the experiment environment, initialized once the
    time that an instance of its subclass maxinterval is created.
    Class attributes are:
    - results_path:        path where store_feat_data() will save data;
    - feat_storing_root_path: path where to bursts features data;
    - pca_results_path:    path where to store latents of PCA applied on ACF, ISI, etc;
    - exp_name:            string of the name of the experiment;
    - step_name:           string of the name of the experimental step.
                           is, those in the list desc_burst_features;
    - burst_features:      the class attribute equal to the global list
    -                      burst_features.
    - exp_steps_dict:      the dict of experimental steps and their timings.
    """
    
    results_path            = None #''
    feat_storing_root_path  = None #''
    pca_results_path        = None #'' 
    exp_name                = None #''
    step_name               = None #''
    burst_features          = None #[]
    exp_steps_dict          = None #{}
    store_results           = False
    
    def __init__(self, *args, **kwargs): pass

    def __call__(self, cls): self.add(cls)

    def add(self, cls):
        """Updating dict bursts_mean_tocluster with the one of cls,
        usually an instance of subclass maxinterval"""
        
        for neuron in cls.neurons_id:
            self.bursts_all_data.setdefault(neuron, []).extend(cls.bursts_all_data[neuron])
            self.bursts_mean_tocluster.setdefault(neuron, []).extend(cls.bursts_mean_tocluster[neuron])
            self.desc_bursts_stat.setdefault(neuron, []).extend(cls.desc_bursts_stat[neuron])
        self._update_info(cls)
    
    def _update_info(self, cls):
        """Updating various maxinterval informations with those
        from cls, usually an instance of subclass maxinterval"""
        
        self.burst_features.extend([z for z in cls.burst_features if z not in self.burst_features])
        self.neurons_id.extend([z for z in cls.neurons_id if z not in self.neurons_id])
        self.neuron_spiking.extend([z for z in cls.neuron_spiking if z not in self.neuron_spiking])
        self.neuron_not_spiking.update(cls.neuron_not_spiking)
    
    def random_burst(self):
        """Randomizing bursts data"""
        
        _it_index = itertools.count()
        data_to_random = sorted([(j[2], _it_index.next(), neuron, j) \
                                for neuron, bursts_data in self.bursts_all_data.iteritems() \
                                for j in bursts_data])
        
        #####################################################################################
        # If using iterator, randomization doesn't happen!!!
        # _it_neurons_in_all_bursts = itertools.imap(operator.itemgetter(2), data_to_random)
        #####################################################################################
        _it_neurons_in_all_bursts = map(operator.itemgetter(2), data_to_random)
        
        t_last_br, i_last_br = data_to_random[0][0], 0
        for br in data_to_random:
            if br[0] - t_last_br >= randomize_window_length:
                shuffled_data = [data_to_random[i] for i in xrange(i_last_br, br[1])]
                numpy.random.shuffle(shuffled_data)
                data_to_random[i_last_br:br[1]] = shuffled_data
                t_last_br, i_last_br = br[0], br[1]
        
        _it_all_bursts = itertools.imap(operator.itemgetter(3), data_to_random) # map(operator.itemgetter(3), data_to_random)
        
        for neuron in self.bursts_all_data: self.bursts_all_data[neuron][:] = []
        for neuron, br in zip(_it_neurons_in_all_bursts, _it_all_bursts):
            self.bursts_all_data[neuron].append(br)
        for neuron in self.bursts_all_data:
            self.bursts_all_data[neuron][:] = sorted(self.bursts_all_data[neuron],
                                                     key = operator.itemgetter(2))
        self.__bursts_stat__()

    def store_feat_data(self,
                        feat_name,
                        step_name,
                        bins_data,
                        data):
        """Storing on .csv file features data of each neuron"""
        
        feat_path = '%s\/%s\/' % (self.feat_storing_root_path, feat_name)
        if not(os.path.exists(feat_path)): os.makedirs(feat_path, 0777)  
        
        feat_filename = '%s%s_%s.csv' % (feat_path, feat_name, step_name)
        data_2_write = map(None, *[['bin:'] + bins_data.tolist()] + \
                                   [[neuron] + list(data.get(neuron, [' '])) \
                                   for neuron in self.neuron_spiking])
        
        csvwriterows_safe(feat_filename, data_2_write)
    
    def binsearch_ts_interval(self, begin_ts, end_ts):
        
        def _inner_func(T): return tuple(T[bisect.bisect_left(T, begin_ts):bisect.bisect(T, end_ts)])
        return _inner_func

    def __repr__(self): return self.__str__()
    
    def __str__(self): return "BurstDetection_root class"



class maxinterval(BurstDetection_root):
    """This class perform burst detection for each neuron or unit, recorded by each channel,
    by replicating the Max Interval Burst Analysis method of Neuroexplorer software 
    (at regard see: http://www.neuroexplorer.com/).
    Once performed the burst detection we calculate various features from timestamps of single unit bursts.
    Here we report such features, together to their acronyms and related brief explanations:
    - index of burst start timestamp;
    - index of burst end timestamp;
    - burst start timestamp;
    - BD: burst time duration;
    - SN: the number of spikes in a burst;
    - BR: the mean burst rate, calculated at end of all bursts detection, 
      because it is the mean rate of all bursts during the entire experimental time interval;  
    - BA: the burst accomodation, or spike frequency accomodation, estimated via the statistical index called
      skewness, that is the statistical standardized third moment, applied on timestamps data of all the spikes
      in a single burst: such index is positive if we have accomodation, that is spikes have bigger rate at the
      begin of the burst, and negative if viceversa the spiking rate is higher towards the end;
    - IBI_low: the number of inter bursts intervals shorter then a given time theshold;
    - IBI_high: the number of inter bursts intervals longer then a given time theshold;
    - RIBI: the ratio IBI_low / (IBI_low + IBI_high), calculated at end of all bursts detection;
    - IBSR: the intra-burst spike rate, that is the mean rate of spikes inside each burst;
    - BK: the statistical index known as kurtosis, that is the statistical standardized fourth moment minus 3, 
      applied on timestamps data of all the spikes in a burst, where a zero value is given by a Gaussian, 
      or a mesokurtic pdf, a positive value by a pdf more peaked and with fatter tails than a Gaussian,
      or a leptokurtic pdf, and a negative values by a distribution less peaked and with tails thinner
      than a Gaussian, or a platykurtic pdf;
    - ACF: the first Principal Component coordinate value of the Autocorrelation functions terms, calculated 
      until a given time delay and timebin length on binned timestamps data, independently to burst detection.
    """
    
    def __init__(self, startreplytime = None, endreplytime = None, all_ts = [], step_name = "", **kwargs):
        
        start = time.time()

        ############################################################################################
        # Various informations about the experiment stored here and once in the root class
        ############################################################################################
        self.start_time, self.end_time = startreplytime, endreplytime
        BurstDetection_root.exp_steps_dict[self.step_name] = (self.start_time, self.end_time)
        self.step_name = step_name
        self.ACF_delay, self.Fano_factor_bin = ACF_delay, Fano_factor_bin
        for k,v in kwargs.iteritems(): setattr(self, k, v)
        self.ACF_bin   = self.ACF_delay / 100.

        if self.results_path:
            self.feat_storing_root_path = self.results_path + 'bursts_data_extracted\/'
            if not(os.path.exists(self.feat_storing_root_path)):
                os.makedirs(self.feat_storing_root_path, 0777)

        ############################################################################################
        # Storing timestamps from transplall matrix, giving attention to have
        # ts strictly < or > in order to not exchange None ts values with
        # startreplytime or endreplytime == 0.0 seconds.
        # The following self.neurons_id is builts with the sam esorting as in
        # transplall matrix, differently to the unsorted tpa_dict.keys().
        # self.neurons_id:              a list of all neurons spiking and not;
        # self.neuron_spiking:          a list only spiking neurons;
        # self.neuron_not_spiking:      a set of only not spking neurons.
        ############################################################################################
        bsearch = self.binsearch_ts_interval(self.start_time , self.end_time) 
        tpa_dict = dict((T[0], bsearch(T[1:])) for T in all_ts)
        self.neurons_id             = [T[0] for T in all_ts]
        self.neuron_spiking         = [T[0] for T in all_ts if len(T) > 1]
        self.neuron_not_spiking     = set(self.neurons_id) - set(self.neuron_spiking)
        self.neurons_no_burst_neither_SN1 = set()

        ############################################################################################
        # self.bursts_all_data:         a dict of all burst data for each neuron;
        # self.bursts_mean_tocluster:   a dict with neurons as keys and as values the 
        #                               statistical indexes of burst features data, 
        #                               for each neuron used in clustering;
        # self.desc_bursts_stat:        a dict with neurons as keys and as values the 
        #                               statistical indexes of descriptive burst features data, 
        #                               for each neuron used in clustering;
        ############################################################################################
        self.bursts_all_data        = {} #dict.fromkeys(self.neuron_spiking)
        self.bursts_mean_tocluster  = {} #dict.fromkeys(self.neuron_spiking)
        self.desc_bursts_stat       = {} #dict.fromkeys(self.neuron_spiking)
        
        if all_ts and startreplytime is not None and endreplytime is not None:
            
            ############################################################################
            # Now finding bursts events from experimental timestamps of each neuron
            ############################################################################
            self.__bursts_scan__(tpa_dict, startreplytime)

            ############################################################################
            # Now calculating features statistical index from bursts data of each neuron
            ############################################################################
            self.__bursts_stat__(tpa_dict, startreplytime, endreplytime)


    def __bursts_scan__(self, tpa_dict, startreplytime):
        """Finding bursts events through a deep inspection of consecutive 
        experimental timestamps of spikes from every neuron"""

        ###################################################################
        # Together to all_features data, we will save also: 'start_index',
        # 'end_index', 'start_timestamp', stored in each burst_dict and so
        # at the begin of each list entry in the finak bursts matrix.
        ###################################################################
        begin_vars = ['start_index', 'end_index', 'start_timestamp']
        self.nr_begin_vars = len(begin_vars)
        self.spikes_nr = dict((neuron, 0) for neuron in self.neurons_id)
        
        all_bursts_ISI  = {}
        all_bursts_IBI  = {}
        all_bursts_ts   = {}
        flag_SN1 = True if 'SN1' in all_features_set else False

        ###print "Performing bursts detection on the interval %.3f-%.3f seconds" %(startreplytime, endreplytime)
        for neuron, s in tpa_dict.iteritems():

            bursts = []
            _bursts_append = bursts.append
            burst_ts, isi = [], []
            _burst_ts_append = burst_ts.append
            burst_dict = dict.fromkeys(all_features)
            nr_ibis_low, nr_ibis_high = 0., 0.
            l, t, end_last_burst = len(s), 0, 0
            self.spikes_nr[neuron] = l
            
            #############################################################
            # t in the while loop is limited to l - 1 = len(s) - 2:
            # by this way we could lose bursts begininning on the before 
            # to last timestamps and finishing on the last one
            #############################################################
            while t < l - 1:
                intra_t = t + 1
                merge_bursts_flag = False
                if s[intra_t] - s[t] > max_t_start:
                    
                    if flag_SN1:
                        burst_dict['start_index']       = t         # start index
                        burst_dict['end_index']         = t         # end index
                        burst_dict['start_timestamp']   = s[t]      # start timestamp
                        burst_dict['SN1']               = 1
                        
                        ##########################################################
                        # Setting IBI_low and IBI_high values for each burst.
                        # Setting ibi_before or the IBI value preceding the burst
                        # Setting variables ibi_low_val and ibi_high_val, 
                        # based on the value fixed for IBI_threshold
                        ##########################################################
                        if end_last_burst > 0: ibi_before = s[t] - s[bursts[-1][1]]
                        elif t == 0: ibi_before = s[t] - startreplytime
                        else: ibi_before = s[t] - s[t-1]
                        all_bursts_IBI.setdefault(neuron, []).append(ibi_before)
                        
                        if 'IBI_low' in all_features_set or 'IBI_high' in all_features_set:
                            if ibi_before <= IBI_threshold:
                                ibi_low_val, ibi_high_val = ibi_before, None 
                                nr_ibis_low += 1
                            else:
                                ibi_low_val, ibi_high_val = None, ibi_before
                                nr_ibis_high += 1 
                            burst_dict['IBI_low']   = ibi_low_val   # IBI_low 
                            burst_dict['IBI_high']  = ibi_high_val  # IBI_high
                        
                        #####################################################################
                        # To burst matrix append: start and end indexes, start time instant,
                        # burst duration, spike nr., burst rate, accomodation by 3rd moment
                        # of burst timestamps, inter burst interval, ibsr, 
                        # kurtosis by 4th of burst timestamps, auto-correlation function
                        # applied to binned timestamps data. 
                        # Giving to SN a zero value, not None, because None is not right 
                        # for the following arithmetic calculations
                        #####################################################################
                        burst_dict['IBI']   = ibi_before    # IBI
                        burst_dict['BD']    = 0             # BD
                        
                        ###########################################################################
                        # Populating bursts with data about burst_features and desc_burst_features
                        ###########################################################################
                        _bursts_append([burst_dict[i] for i in begin_vars + all_features])
                    
                    t += 1
                    continue
                else:
                    if s[intra_t] - s[t] > max_t_end:
                        t += 1
                        continue
                    else:
                        while intra_t < l - 1 and s[intra_t + 1] - s[intra_t] <= max_t_end: intra_t += 1
                        if end_last_burst > 0 and s[t] - s[end_last_burst] < min_dur_interbursts:
                            t, merge_bursts_flag = start_last_burst, True
                        if s[intra_t] - s[t] >= min_dur_burst and \
                           intra_t + 1 - t >= min_nr_spikes_per_burst:
                            ############################################################
                            # last_burst_ts contains the delays of timestamps of spikes
                            # after the first one from the timestamp of the latter
                            ############################################################
                            ts_for_moments = s[t:intra_t + 1]
                            last_burst_ts = [lts - ts_for_moments[0] for lts in ts_for_moments[1:]]
                            _burst_ts_append(last_burst_ts)
                            all_bursts_ts.setdefault(neuron, []).extend(last_burst_ts)
                            isi = isi_return(ts_for_moments)
                            all_bursts_ISI.setdefault(neuron, []).extend(isi)
                            
                            ##########################################################
                            # Setting ibi_before or the IBI value preceding the burst
                            ##########################################################
                            if end_last_burst > 0:
                                if not(merge_bursts_flag): ibi_before = s[t] - s[bursts[-1][1]]
                                else: ibi_before = all_bursts_IBI[neuron][-1]                                    
                            elif t == 0: ibi_before = s[t] - startreplytime
                            else: ibi_before = s[t] - s[t-1]
                            all_bursts_IBI.setdefault(neuron, []).append(ibi_before)
                            
                            ##################################################
                            # Setting variables ibi_low_val and ibi_high_val, 
                            # based on the value fixed for IBI_threshold
                            ##################################################
                            if 'IBI_low' in all_features_set or 'IBI_high' in all_features_set:
                                if ibi_before <= IBI_threshold:
                                    ibi_low_val, ibi_high_val = ibi_before, None 
                                    nr_ibis_low += 1
                                else:
                                    ibi_low_val, ibi_high_val = None, ibi_before
                                    nr_ibis_high += 1
                                burst_dict['IBI_low']   = ibi_low_val   # IBI_low 
                                burst_dict['IBI_high']  = ibi_high_val  # IBI_high
                            
                            #####################################################################
                            # To burst matrix append: start and end indexes, start time instant,
                            # burst duration, spike nr., burst rate, accomodation by 3rd moment
                            # of burst timestamps, inter burst interval, ibsr, 
                            # kurtosis by 4th of burst timestamps, auto-correlation function
                            # applied to binned timestamps data. 
                            #####################################################################
                            burst_dict['start_index']       = t                 # start index
                            burst_dict['end_index']         = intra_t           # end index
                            burst_dict['start_timestamp']   = s[t]              # start timestamp
                            burst_dict['SN1']               = intra_t + 1 - t   # SN1 
                            burst_dict['SN']                = intra_t + 1 - t   # SN 
                            burst_dict['BD']                = s[intra_t] - s[t] # BD
                            burst_dict['IBSR']              = (intra_t + 1 - t) / (s[intra_t] - s[t])   # IBSR
                            burst_dict['IBI']               = ibi_before        # IBI
                                
                            ###################################################################
                            # Storing in the matrix burst only the values of burst_dict, 
                            # so leaving his dictionary data structure. Values are copied from 
                            # each dict burst_dict in the same order as that in the list 
                            # variable burst_features read from burst_detection_settings.txt.
                            # Populating bursts with data about all_features.
                            ###################################################################
                            if merge_bursts_flag: bursts[-1] = [burst_dict[i] for i in begin_vars + all_features]
                            else: _bursts_append([burst_dict[i] for i in begin_vars + all_features])

                            start_last_burst, end_last_burst = t, intra_t
                            t = intra_t + 1
                            continue
                        else:
                            t += 1
                            continue

            ##################################################################
            # Storing all bursts data just calculated for the related neuron,
            # by copying that from the matrix named bursts.
            # Then distinguish if burst matrix is empty or not in order to
            # calculate data to cluster.
            ##################################################################
            self.bursts_all_data[neuron] = bursts
            
    def __bursts_stat__(self, tpa_dict = None, startreplytime = None, endreplytime = None):
        """Calculating statistical index of features values obtained
        from bursts data of each neuron"""

        #########################################################
        # Building a dict with all_features as keys, in order to 
        # check faster if a feature is used, and Fixing SN1 flag 
        #########################################################
        #stat_feat = stat_features if len(stat_features) == len(burst_features) \
        #            else [meanval] * len(burst_features)
        #stat_desc_feat = stat_desc_features if len(stat_desc_features) == len(desc_burst_features) \
        #                 else [meanval] * len(desc_burst_features)
        stat_feat = [meanval] * len(burst_features)
        stat_desc_feat = [meanval] * len(desc_burst_features)
        stat_dict_all_feat = dict(zip(all_features, stat_feat + stat_desc_feat))
        
        if tpa_dict:

            ############################################################
            # Calculating ACFs data.
            # Calculating the first three principal components of ACF.
            # Calculating CCF for all pairs of neurons if requested.
            ############################################################
            self.__ACF__(tpa_dict)
            if set(['ACF_PC1', 'ACF_PC2', 'ACF_PC3']) & all_features_set: acf_pca_score = self.__PCA_ACF__()
            if CCF_on: self.__CCF__(tpa_dict)

            ####################################################################
            # Calculating ISIs data.
            # Calculating the first three principal components of ISI histogram 
            ####################################################################
            all_ISI = dict((neuron, numpy.asarray(isi_return(tpa_dict[neuron]))) for neuron in self.neuron_spiking)
            self.__ISI__(all_ISI)

            if set(['ISI_PC1', 'ISI_PC2', 'ISI_PC3']) & all_features_set:
                isi_pca_score = self.__PCA_ISI__()

            ##################################################################################
            # Creating variables for the first three principal components of histogram of IBI 
            ##################################################################################
            if 'IBI' in all_features_set and \
                set(['IBI_PC1', 'IBI_PC2', 'IBI_PC3']) & all_features_set:
                exp_step_ibi_bins, self.exp_step_ibis = numpy.arange(0, IBI_delay + IBI_bin, IBI_bin), {}


        for neuron, bursts in self.bursts_all_data.iteritems():

            if bursts == []:
    
                self.bursts_mean_tocluster[neuron]  = [None] * len(burst_features)
                self.desc_bursts_stat[neuron]       = [None] * len(desc_burst_features)
                self.neurons_no_burst_neither_SN1.add(neuron)
            
            else:
              
                bursts_dict_all_feat = dict.fromkeys(all_features)
                bursts_transpose = map(None, *bursts)
                len_brs = len(bursts)
                
                if startreplytime is not None and endreplytime is not None:

                    ##########################################################################################
                    # Spike rate SR as nr of spikes divided by the temporal duration of the experimental step
                    # For the burst rate BR we calculate the mean burst rate at minute (* 60 seconds)
                    ##########################################################################################
                    if 'SR' in all_features_set:
                        bursts_dict_all_feat['SR'] = self.spikes_nr[neuron] / (endreplytime - startreplytime)
                    if 'BR' in all_features_set:
                        bursts_dict_all_feat['BR'] = len_brs * 60. / (endreplytime - startreplytime)
                
                if tpa_dict:

                    ###################################################################
                    # 1st and/or 2nd, 3rd PCs of ACF values, ISI and IBI distributions
                    ###################################################################
                    acf_scores = acf_pca_score.get(neuron, [[0], [0], [0]])
                    if 'ACF_PC1' in all_features_set: bursts_dict_all_feat['ACF_PC1'] = acf_scores[0]
                    if 'ACF_PC2' in all_features_set: bursts_dict_all_feat['ACF_PC2'] = acf_scores[1]
                    if 'ACF_PC3' in all_features_set: bursts_dict_all_feat['ACF_PC3'] = acf_scores[2]
                    
                    if 'ISI_PC1' in all_features_set: bursts_dict_all_feat['ISI_PC1'] = isi_pca_score[neuron][0]
                    if 'ISI_PC2' in all_features_set: bursts_dict_all_feat['ISI_PC2'] = isi_pca_score[neuron][1]
                    if 'ISI_PC3' in all_features_set: bursts_dict_all_feat['ISI_PC3'] = isi_pca_score[neuron][2]
                    
                    if 'IBI' in all_features_set and \
                        set(['IBI_PC1', 'IBI_PC2', 'IBI_PC3']) & all_features_set:
                        IBI_index = all_features.index('IBI') + self.nr_begin_vars
                        self.exp_step_ibis[neuron] = numpy.histogram(bursts_transpose[IBI_index], \
                                                                     bins = exp_step_ibi_bins)[0]
                        if IBI_shape_normalize:
                            self.exp_step_ibis[neuron] = numpy.asarray(shape_normalize(self.exp_step_ibis[neuron]))
            
                    ############################################################
                    # Calculating two indexes from ACFs data of each neuron:
                    # - the ratio between the max and the min of ACF
                    # - the bin nr corresponding to half life period of ACF
                    ############################################################
                    if 'ratio_ACF' in all_features_set:
                        _max_ACF = max(self.exp_step_acf[neuron])
                        for _last_ACF in self.exp_step_acf[neuron][::-1]:
                            if _last_ACF: break
                        else: _last_ACF = 0.
                        bursts_dict_all_feat['ratio_ACF'] = _max_ACF / _last_ACF if _last_ACF else 0.
                        
                    if 'halflife_ACF' in all_features_set:
                        _half_max_ACF = max(self.exp_step_acf[neuron]) / 2.
                        for _i, _acf in enumerate(self.exp_step_acf[neuron]): 
                            if _acf >= _half_max_ACF: bursts_dict_all_feat['halflife_ACF'] = _i
             
                ###################################################################################
                # Calculating BK value for intra-burst timestamps.
                # Calculating SFA value for intra-burst timestamps.
                # Calculating CV value for intra-burst timestamps.
                # Calculating CV value for all ISIs/timestamps in the experimental step.
                # Calculating the kurtosis of ISIs distribution in the experimental step.
                # Calculating the skewness of ISIs distribution in the experimental step.
                # Calculating Fano factor value for intra-burst timestamps.
                # It's mandatory to transform the iterator in a list, when applying Fano function, 
                # otherwise iterator's data will be lost after its first processing. 
                # Calculating Fano factor value for all binned timestamps in the experimental step
                ###################################################################################
                if 'BK' in all_features_set: bursts_dict_all_feat['BK'] = kurtosis(isi)
                if 'SFA' in all_features_set: bursts_dict_all_feat['SFA'] = skewness(isi)
                if 'CV' in all_features_set: bursts_dict_all_feat['CV'] = coeff_variation(isi)
                if 'Fano' in all_features_set:
                    #fano_bins = [numpy.histogram(br_ts, \
                    #                             bins = numpy.arange(0, max(br_ts) + self.Fano_factor_bin, self.Fano_factor_bin))[0] \
                    #                             for br_ts in burst_ts]
                    #bursts_dict_all_feat['Fano'] = Fano_factor(list(itertools.chain(*fano_bins)))
                    bursts_dict_all_feat['Fano'] = Fano_factor(numpy.hstack([numpy.histogram(br_ts, \
                                                        bins = numpy.arange(0,max(br_ts)+self.Fano_factor_bin, self.Fano_factor_bin))[0] \
                                                        for br_ts in burst_ts]))
                
                if tpa_dict:

                    if 'ISI_K' in all_features_set: bursts_dict_all_feat['ISI_K'] = kurtosis(all_ISI[neuron])
                    if 'ISI_S' in all_features_set: bursts_dict_all_feat['ISI_S'] = skewness(all_ISI[neuron])
                    if 'CV_all' in all_features_set: bursts_dict_all_feat['CV_all'] = coeff_variation(all_ISI[neuron])
                    if 'CV_all_SQ' in all_features_set: bursts_dict_all_feat['CV_all_SQ'] = coeff_variation(all_ISI[neuron]) ** 2
                    if 'CV_all_cut' in all_features_set:
                        bursts_dict_all_feat['CV_all_cut'] = coeff_variation(all_ISI[neuron],
                                                                            upper_threshold = ISI_upper_threshold,
                                                                            lower_threshold = ISI_lower_threshold)
                    if 'CV2_all' in all_features_set:
                        bursts_dict_all_feat['CV2_all'] = coeff_variation2(all_ISI[neuron],
                                                                            upper_threshold = ISI_upper_threshold,
                                                                            lower_threshold = ISI_lower_threshold)
                
                if tpa_dict and startreplytime is not None and endreplytime is not None:
                    
                    if 'Fano_all' in all_features_set:
                        bursts_dict_all_feat['Fano_all'] = \
                            Fano_factor(numpy.histogram(tpa_dict[neuron], \
                                            bins = numpy.arange(startreplytime, endreplytime, self.Fano_factor_bin))[0])
                        #                    self.Fano_factor_bin))[0].tolist())
                    if 'Fano_all_cut' in all_features_set:
                        bursts_dict_all_feat['Fano_all_cut'] = \
                            Fano_factor(numpy.histogram(tpa_dict[neuron], \
                                            bins = numpy.arange(startreplytime, endreplytime, self.Fano_factor_bin))[0], \
                                        lower_threshold = FF_lower_threshold, upper_threshold = FF_upper_threshold)
                        #                    self.Fano_factor_bin))[0].tolist(), \
                
                #####################################
                # BDK: BD kurtosis, BDS: BD skewness 
                # SNK: SN kurtosis, SNS: SN skewness 
                #####################################
                if any(map(all_features.__contains__, ['BDK', 'BDS'])):
                    BD_index = all_features.index('BD') + self.nr_begin_vars
                if 'BDK' in all_features_set: bursts_dict_all_feat['BDK'] = kurtosis(bursts_transpose[BD_index])
                if 'BDS' in all_features_set: bursts_dict_all_feat['BDS'] = skewness(bursts_transpose[BD_index])
                
                if any(map(all_features.__contains__, ['SNK', 'SNS'])):
                    SN_index = all_features.index('SN') + self.nr_begin_vars
                if 'SNK' in all_features_set:  bursts_dict_all_feat['SNK'] = kurtosis(bursts_transpose[SN_index])
                if 'SNS' in all_features_set:  bursts_dict_all_feat['SNS'] = skewness(bursts_transpose[SN_index])

                ############################################################
                # Calculating RIBI value
                # Calculating statistical index for IBI_low and/or IBI_high
                ############################################################
                if 'RIBI' in all_features_set:
                    bursts_dict_all_feat['RIBI'] = nr_ibis_low / (nr_ibis_low + nr_ibis_high) \
                                                   if nr_ibis_low + nr_ibis_high > 0 else None
                if 'IBI_low' in all_features_set:
                    IBI_low_index = all_features.index('IBI_low') + self.nr_begin_vars
                    if len_brs > 1:
                        if not(any(bursts_transpose[IBI_low_index])): bursts_dict_all_feat['IBI_low'] = 0
                        else: bursts_dict_all_feat['IBI_low'] = stat_dict_all_feat['IBI_low'](bursts_transpose[IBI_low_index])
                    elif len_brs == 1:
                        if not(bursts_transpose[IBI_low_index]): bursts_dict_all_feat['IBI_low'] = 0
                        else: bursts_dict_all_feat['IBI_low'] = bursts_transpose[IBI_low_index]
                if 'IBI_high' in all_features_set:
                    IBI_high_index = all_features.index('IBI_high') + self.nr_begin_vars
                    if len_brs > 1:
                        if not(any(bursts_transpose[IBI_high_index])): bursts_dict_all_feat['IBI_high'] = 0
                        else: bursts_dict_all_feat['IBI_high'] = stat_dict_all_feat['IBI_high'](bursts_transpose[IBI_high_index])
                    elif len_brs == 1:
                        if not(bursts_transpose[IBI_high_index]): bursts_dict_all_feat['IBI_high'] = 0
                        else: bursts_dict_all_feat['IBI_high'] = bursts_transpose[IBI_high_index]

                #######################################################################
                # Calculating statistical index for ('BD', 'SN', 'SN1', 'IBSR', 'IBI')
                # for the ones present in all_features list
                #######################################################################
                for feat_to_stat in features_to_histo & set(all_features):
                    feat_to_stat_index = all_features.index(feat_to_stat) + self.nr_begin_vars
                    bursts_dict_all_feat[feat_to_stat] = \
                            stat_dict_all_feat[feat_to_stat](bursts_transpose[feat_to_stat_index])

                ################################################################################
                # Populating bursts_transpose and then bursts_data with features for which any
                # statistical index has been not yet calculated (so not the previous features).
                # The following count is necessary for the features present more than one time
                # in all_features, in order to well fill column data in *_burst_data.csv file.
                ################################################################################
                for feat_to_stat in (i for i in all_features if i not in features_to_histo):
                    feat_to_stat_indexes, nr_feat_count, previous_feat_to_stat_index = [], 0, 0
                    while nr_feat_count < all_features.count(feat_to_stat):
                        feat_to_stat_index = all_features.index(feat_to_stat, previous_feat_to_stat_index) + \
                                             self.nr_begin_vars
                        if feat_to_stat_index not in feat_to_stat_indexes:
                            feat_to_stat_indexes.append(feat_to_stat_index)
                        if len_brs > 1: bursts_transpose[feat_to_stat_index] = \
                                            [bursts_dict_all_feat[feat_to_stat]] * len_brs
                        else: bursts_transpose[feat_to_stat_index] = bursts_dict_all_feat[feat_to_stat]
                        nr_feat_count += 1
                        previous_feat_to_stat_index = feat_to_stat_index + 1 - self.nr_begin_vars

                ################################
                # Updating self.bursts_all_data
                ################################
                if len_brs > 1: self.bursts_all_data[neuron] = map(list, map(None, *bursts_transpose))
                else: self.bursts_all_data[neuron] = [bursts_transpose]

                ########################################################################
                # Filling the following dicts with neurons as keys with the statistical
                # indexes of features both from burst_features and desc_burst_features
                ########################################################################
                self.bursts_mean_tocluster[neuron]  = [bursts_dict_all_feat[i] for i in burst_features]
                self.desc_bursts_stat[neuron]       = [bursts_dict_all_feat[i] for i in desc_burst_features]
        
        #############################################################################
        # Updating self.neuron_spiking with neurons in self.neurons_no_burst_neither_SN1,
        # that is, that who didn't bursted neither spiked in SN1 mode
        #############################################################################
        self.neuron_spiking = [neuron for neuron in self.neuron_spiking \
                               if not(self.neurons_no_burst_neither_SN1.__contains__(neuron))]
        self.neuron_not_spiking -= self.neurons_no_burst_neither_SN1
        
        ####################################################
        # The first three PCs of IBI histogram distribution
        ####################################################
        if tpa_dict and 'IBI' in all_features_set and \
           set(['IBI_PC1', 'IBI_PC2', 'IBI_PC3']) & all_features_set:
            pca_exp_step_ibis = pca.pca([self.exp_step_ibis[neuron] for neuron in self.neuron_spiking])
            ibi_pca_score = dict(zip(self.neuron_spiking, pca_exp_step_ibis.SCORE)) 
            
            for ibi_pca_nr, feat_to_stat in enumerate(('IBI_PC1', 'IBI_PC2', 'IBI_PC3')):
                if feat_to_stat in burst_features:
                    feat_to_stat_index = burst_features.index(feat_to_stat)
                    for neuron, mean_tocluster in self.bursts_mean_tocluster.items():
                        if neuron in self.exp_step_ibis:
                            mean_tocluster[feat_to_stat_index] = ibi_pca_score[neuron][ibi_pca_nr]
                            bursts_transpose = map(None, *self.bursts_all_data[neuron])
                            len_brs = len(self.bursts_all_data[neuron])
                            bursts_transpose[feat_to_stat_index + self.nr_begin_vars] = \
                                                                [mean_tocluster[feat_to_stat_index]] * len_brs
                            ################################
                            # Updating self.bursts_all_data
                            ################################
                            self.bursts_all_data[neuron] = map(list, map(None, *bursts_transpose)) \
                                                            if len_brs > 1 else [bursts_transpose]
                
                if feat_to_stat in desc_burst_features:
                    feat_to_stat_index = desc_burst_features.index(feat_to_stat)
                    for neuron, desc_mean_tocluster in self.desc_bursts_stat.items():
                        if neuron in self.exp_step_ibis:
                            desc_mean_tocluster[feat_to_stat_index] = ibi_pca_score[neuron][ibi_pca_nr]
                            bursts_transpose = map(None, *self.bursts_all_data[neuron])
                            len_brs = len(self.bursts_all_data[neuron])
                            bursts_transpose[feat_to_stat_index + self.nr_begin_vars] = \
                                                                [desc_mean_tocluster[feat_to_stat_index]] * len_brs
                            ################################
                            # Updating self.bursts_all_data
                            ################################
                            self.bursts_all_data[neuron] = map(list, map(None, *bursts_transpose)) \
                                                            if len_brs > 1 else [bursts_transpose]
            
            ########################################
            # Storing IBI histogram for each neuron 
            ########################################
            if self.store_results:
                self.store_feat_data('IBI',
                                    self.step_name, 
                                    exp_step_ibi_bins[1:],
                                    self.exp_step_ibis)
            
            ##################################
            # Storing PCA latent data for IBI
            ##################################
            if store_PCAdata:
                pca_results_path = self.results_path + 'PCA_data\/'
                if not(os.path.exists(pca_results_path)): os.makedirs(pca_results_path, 0777)
                pca_writer = csv.writer(open(pca_results_path + self.step_name + '__IBI_latent.csv', 'wb'), **csv_options)
                pca_writer.writerow(list(pca_exp_step_ibis.latent))
        
        ########################################################################################
        # Calculating histogram data of features_to_histo for each neuron.
        # self.start_end_bin is dict which will contain for each feat_to_histo a tuple of the
        # following three values: the value of the first and the last bins, and bin's duration.
        # The dict self.all_feat_to_histo_data will contain the histogram data of all neurons
        # and for each feature, whcih is here acting as the key of the same dict.
        ########################################################################################
        if self.store_results: self.__histo_data_features__()

        self.neurons_keys = self.bursts_all_data.iterkeys()
        end = time.time()
    

    def __histo_data_features__(self):

        self.start_end_bin, self.all_feat_to_histo_data = {}, {}

        for feat_to_histo in features_to_histo & set(all_features):
            feat_to_histo_index = all_features.index(feat_to_histo) + self.nr_begin_vars
                
            ############################################################################
            # If the feature is 'IBSR' or 'SN', which are None when there're no bursts,
            # it is needed to filter values None from the data array.  
            ############################################################################
            if feat_to_histo == 'IBSR' or feat_to_histo == 'SN' :
                feat_to_histo_data = dict((neuron,
                                           numpy.asarray([ibsr for ibsr in map(operator.itemgetter(feat_to_histo_index), \
                                                                               self.bursts_all_data[neuron]) if ibsr is not None])) \
                                           for neuron in self.neuron_spiking)
            else: feat_to_histo_data = dict((neuron,
                                             numpy.asarray(map(operator.itemgetter(feat_to_histo_index),
                                                               self.bursts_all_data[neuron])))
                                             for neuron in self.neuron_spiking)
                
            #############################################################################################
            # When feature is 'IBSR', the value of the first bin is usually not zero, as for the others.
            # Creating the bins array from the values of self.start_end_bin[feat_to_histo].
            #############################################################################################
            if feat_to_histo_data:
                feat_to_histo_data_conc = numpy.concatenate(feat_to_histo_data.values())
                self.start_end_bin[feat_to_histo] = (numpy.floor(feat_to_histo_data_conc.min()) 
                                                     if feat_to_histo == 'IBSR' else 0,
                                                     feat_to_histo_data_conc.max() + 
                                                     features_to_histo_bin_dict[feat_to_histo],
                                                     features_to_histo_bin_dict[feat_to_histo])
            else: self.start_end_bin[feat_to_histo] = (0., 0., features_to_histo_bin_dict[feat_to_histo])

            feat_to_histo_bins = numpy.arange(self.start_end_bin[feat_to_histo][0],
                                              self.start_end_bin[feat_to_histo][1],
                                              self.start_end_bin[feat_to_histo][2])

            for neuron in self.neuron_spiking:
                feat_to_histo_data[neuron] = numpy.histogram(feat_to_histo_data[neuron],
                                                             feat_to_histo_bins)[0]
            
            self.all_feat_to_histo_data[feat_to_histo] = feat_to_histo_data
                
            if self.store_results and store_allburstsdata:
                self.store_feat_data(feat_to_histo,
                                     self.step_name,
                                     feat_to_histo_bins[:-1],
                                     feat_to_histo_data)
                                               

    def __repr__(self): return self.__str__()
    def __str__(self): return "maxinterval class"

    def __ACF__(self, tpa_dict):
        """Calculating ACFs data"""

        ###########################################################
        # We won't process neurons with empty timestamps trains by 
        # imposing that neurons comes from 'self.neuron_spiking'
        ###########################################################
        _dict_data = dict((neuron, tpa_dict[neuron]) for neuron in self.neuron_spiking)
        self.exp_step_acf = neurocorr.neurons_acf(_dict_data,
                                                  bin = self.ACF_bin,
                                                  delta = self.ACF_delay,
                                                  normalize = ACF_shape_normalize,
                                                  no_AC0 = ACF_no_AC0,
                                                  smoothing_on = ACF_smoothing,
                                                  pdf = ACF_pdf,
                                                  pdf_samples_nr = pdf_samples_nr)
        
        ###################################
        # Storing ACF data for each neuron
        ###################################
        if self.store_results and store_allburstsdata:
            self.store_feat_data('ACF',
                                 self.step_name, 
                                 numpy.arange(self.ACF_bin, self.ACF_delay + self.ACF_bin, self.ACF_bin), 
                                 self.exp_step_acf)

    def __PCA_ACF__(self):

        ###########################################################
        # We won't process neurons with empty timestamps trains by 
        # imposing that neurons comes from 'self.neuron_spiking'
        ###########################################################
        _acfs = self.exp_step_acf
        maxl = max(map(len, _acfs.itervalues()))
        pca_exp_step_acf = pca.pca([list(_acfs[i])+[0]*(maxl-len(_acfs[i])) for i in self.neuron_spiking])
        
        ##################################
        # Storing PCA latent data for ACF
        ##################################
        if store_PCAdata:
            pca_results_path = self.results_path + 'PCA_data\/'
            if not(os.path.exists(pca_results_path)): os.makedirs(pca_results_path, 0777)
            pca_writer = csv.writer(open(pca_results_path + self.step_name + '__ACF_latent.csv', 'wb'), **csv_options)
            pca_writer.writerow(list(pca_exp_step_acf.latent))
        
        return dict(zip(self.neuron_spiking, pca_exp_step_acf.SCORE))

    def __CCF__(self, tpa_dict):
        """Calculating CCFs data"""
            
        CCFs = neurocorr.neurons_ccf(dict((neuron, tpa_dict[neuron]) for neuron in self.neuron_spiking),
                                     bin = CCF_bin, 
                                     delta = CCF_delay)
        CCF_csvname = self.feat_storing_root_path + 'CCF-' + self.step_name + '.csv'
        CCF_csvwriter = csv.writer(open(CCF_csvname, 'wb'), **csv_options)
        CCF_csvwriter.writerows([i[0][0], i[0][1], ':'] + i[1].tolist() for i in sorted(CCFs.iteritems()))

    def __ISI__(self, all_ISI):
        """Calculating ISIs data"""

        ########################
        # Calculating ISIs data
        ########################
        exp_step_isi_bins = numpy.arange(0, ISI_delay + ISI_bin, ISI_bin)
            
        ##############################################################
        # We won't process neurons with empty timestamps trains by 
        # imposing that neurons comes from 'self.neuron_spiking'
        ##############################################################
        self.exp_step_isi_histo = dict((neuron, numpy.histogram(all_ISI[neuron], \
                                                                bins = exp_step_isi_bins)[0]) \
                                        for neuron in self.neuron_spiking)
          
        #############################################
        # Normalizing self.exp_step_acf if requested
        #############################################
        if ISI_shape_normalize:  self.exp_step_isi_histo = dict((i, numpy.asarray(shape_normalize(j))) \
                                                                for i, j in self.exp_step_isi_histo.iteritems())
                
        #########################################
        # Storing all ISI values for each neuron  
        #########################################
        if self.store_results and store_allburstsdata:
            self.store_feat_data('ISI',
                                self.step_name,
                                exp_step_isi_bins[1:],
                                self.exp_step_isi_histo)
    
    def __PCA_ISI__(self):

        ###########################################################
        # We won't process neurons with empty timestamps trains by 
        # imposing that neurons comes from 'self.neuron_spiking'
        ###########################################################
        pca_exp_step_isis = pca.pca([self.exp_step_isi_histo[neuron] for neuron in self.neuron_spiking])
        
        ##################################
        # Storing PCA latent data for ISI
        ##################################
        if store_PCAdata:
            pca_results_path = self.results_path + 'PCA_data\/'
            if not(os.path.exists(pca_results_path)): os.makedirs(pca_results_path, 0777)
            pca_writer = csv.writer(open(pca_results_path + self.step_name + '__ISI_latent.csv', 'wb'), **csv_options)
            pca_writer.writerow(list(pca_exp_step_isis.latent))

        return dict(zip(self.neuron_spiking, pca_exp_step_isis.SCORE))



class netbursts(BurstDetection_root):
    """Algorithm for the detection of network bursts, as explained in the reference:
       - Michael I. Ham, Luis M. Bettencourt, Floyd D. McDaniel, Guenter W. Gross. 
         Spontaneous coordinated activity in cultured networks: Analysis of multiple
         ignition sites primary circuits, and burst phase delay distributions.
         Journal of Computational Neuroscience, 24 (3), 2008, p.346.
       Data extracted about network bursts thanks to the algorithm above are:
       - index of the burst start timestamp
       - name of the burst leader neuron
       - burst start timestamps
       - index of the burst end timestamp
       - burst duration
       - number of neurons in the burst
       - timestamps and related neuron name of the first timestamps between leader's followers
    """

    @staticmethod
    def _neurons_set(L):
        """Return the set of neurons spiking in the list of tuples L,
        whose every tuple is of the kind (neuron, timestamps)
        """
        
        return set(zip(*t_interval)[0])
    
    @staticmethod
    def _neurons_ranking(L, max_rank = None):
        """Return the list of the first spike and its delay from
        the burst leader fired by the first max_rank neurons
        """

        if not(max_rank): max_rank = len(_neurons_set(L))
        neurons_ranked_set = set()
        ranking = []
        rank = 0
        for i, j in enumerate(L): 
            if j[0] not in neurons_ranked_set:
                ranking.append((j[0], rank, j[1] - L[0][1], i))
                neurons_ranked_set.add(j[0])
                rank += 1
            if rank == max_rank: break
        return ranking
     
    def _tstuples(self, T): return zip(itertools.repeat(T[0]), T[1:])

    def ts_in_interval(self, ts):
        
        if factory_cmp(self.start_time, self.end_time)(ts): return ts
        return False    
    
    def _filter_neuron_ts(self, T):
        
        bsearch = self.binsearch_ts_interval(self.start_time, self.end_time) 
        return zip(itertools.repeat(T[0]), bsearch(T[1:]))
    
    def all_neuron_ts_tuples(self, all_T):
        
        filter_neuron_ts = self._filter_neuron_ts
        return sorted(itertools.chain(*map(filter_neuron_ts, all_T)), key = operator.itemgetter(1))

    def __init__(self, start_time, end_time, all_ts):

        self.start_time, self.end_time = start_time, end_time

        #####################################################################
        # Storing timestamps data from transplall matrix.
        # Give attention to give ts strictly > or < in order to not exchange 
        # Null ts values with startreplytime or endreplytime == 0.0 seconds.
        # Data arrays all_ts_tuples and all_tstuple_per_neuron
        # are composed of tuples with as first element the neuron name and 
        # as the second the related timestamps.
        # all_ts_tuples is the same as all_tstuple_per_neuron,
        # but sorted based on timestamps.
        #####################################################################
        all_ts_tuples = self.all_neuron_ts_tuples(all_ts)
        
        self.bursts_container = []
        allts_nr = len(all_ts_tuples)
        z_begin, PBoB_counter, previous_Nethigh_threshold = 0, 0, False
        t_end_lastNB = -1000000
        
        while z_begin < allts_nr - Nethigh_threshold + 1:
            t_begin = all_ts_tuples[z_begin][1]

            #######################################################################
            # Scanning a single timebin, in order to stop the search when the isi 
            # is > Netburst_binlength , by using the Python construct for...else
            # The else clause is basilar, elsewhere the while cycle couldn't stop
            #######################################################################
            for z in xrange(z_begin + 1, allts_nr):
                if all_ts_tuples[z][1] - t_begin > Netburst_binlength: z_end = z - 1; break
            else: z_end = allts_nr - 1

            ###############################################################
            # If the end spike of the previous burst has been spiked more
            # than min_interburst_length seconds ago, clear the variables:
            # PBoB_counter, previous_Nethigh_threshold, PBoB_dict.
            ###############################################################
            if z_begin and t_begin - all_ts_tuples[z_begin - 1][1] > min_interburst_length \
                    or t_begin - t_end_lastNB > min_interburst_length:
                PBoB_counter, previous_Nethigh_threshold, PBoB_dict = 0, False, {}
                #z_begin += 1
                #continue
            
            #################################################################
            # If there is a nr of spikes > high_threshold during the timebin
            #################################################################
            if z_end - z_begin + 1 >= Nethigh_threshold:

                #####################################################################
                # If in the previous timebin there was nr of spikes > high_threshold
                #####################################################################
                if previous_Nethigh_threshold:
                    self.bursts_container[-1]['end_index']     = z_end
                    self.bursts_container[-1]['end_time']      = all_ts_tuples[z_end][1]
                    self.bursts_container[-1]['spikes_events'] = \
                            all_ts_tuples[self.bursts_container[-1]['start_index']:z_end + 1]
                #####################################################################################
                # If in the previous timebin there was low_threshold < nr of spikes < high_threshold
                #####################################################################################
                elif PBoB_counter > 0:
                    PBoB_dict['end_index']     = z_end
                    PBoB_dict['end_time']      = all_ts_tuples[z_end][1]
                    PBoB_dict['spikes_events'] = all_ts_tuples[PBoB_dict['start_index']:z_end + 1]
                    self.bursts_container.append(PBoB_dict)
                ####################################################################
                # If in the previous timebin there was nr of spikes < low_threshold
                ####################################################################
                else:
                    self.bursts_container.append({'start_index':   z_begin, \
                                                  'end_index':     z_end, \
                                                  'start_time':    t_begin, \
                                                  'end_time':      all_ts_tuples[z_end][1], \
                                                  'spikes_events': all_ts_tuples[z_begin:z_end + 1]})
                PBoB_counter, previous_Nethigh_threshold = 0, True
                t_end_lastNB =  all_ts_tuples[z_end][1]

            #################################################################################
            # If there is a low_threshold < nr of spikes < high_threshold during the timebin
            #################################################################################
            elif Netlow_threshold <= z_end - z_begin + 1 < Nethigh_threshold:
                
                #####################################################################
                # If in the previous timebin there was nr of spikes > high_threshold
                #####################################################################
                if previous_Nethigh_threshold:
                    self.bursts_container[-1]['end_index']     = z_end
                    self.bursts_container[-1]['end_time']      = all_ts_tuples[z_end][1]
                    self.bursts_container[-1]['spikes_events'] = \
                            all_ts_tuples[self.bursts_container[-1]['start_index']:z_end + 1]
                    t_end_lastNB =  all_ts_tuples[z_end][1]
                ####################################################################
                # If in the previous timebin there was nr of spikes < low_threshold
                ####################################################################
                elif PBoB_counter == 0:
                    PBoB_dict = ({'start_index':   z_begin, \
                                  'start_time':    t_begin, \
                                  'end_index':     None, \
                                  'end_time':      None, \
                                  'spikes_events': None})
                PBoB_counter += 1            

            z_begin = z_end + 1
        #print len(self.bursts_container)
