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

from importlibs import *
import burstdetection
import statf.dataprocflow as procflow
from storage_functions import pd_cluster_data
from statf.signalprocessing import smooth


__author__ = "Andrea Maffezzoli - modified by Antonella Farinaccio"
__revision__ = "$Revision: 1.37 - 24 Jun 2011$"
__name__ = "netburst_processing"
__all__ = ["netburst"]
numpy.seterr(divide='raise')
numpy.seterr(invalid='raise')


class netburst(object):

    store_results, data_path = False, ""
    trash_poor_states = None
    Net_features_bin = {}
    Net_features = []
    
    def __init__(self, **kwargs):

        self.step_name, self.start_time, self.end_time = kwargs.get('_step_info', [None] * 3)
        self.bin, self.delay = self.Net_features_bin.get('bin'), self.Net_features_bin.get('delay')
        self.states_nr = self.Net_features_bin['states_nr']
        self.do_Ttest = self.Net_features_bin['do_Ttest']

    def detection(self, *args): return burstdetection.netbursts(*args)
    
    def __call__(self, all_ts, **kwargs):
        """Input:
           - all_ts: timestamps matrix for any recorded channel
           - clusters_members: a list of iterables giving neurons membership to clusters 
             (in clusters_members[2]), as instances of attribute clusters_names_sorted_x_nr
             from method __clustering__ of class single_neuron_burst from module timeneuro_mod.py.
        """
        
        self.__process__(all_ts, **kwargs)

    @staticmethod
    def _mean(x):
        
        if not(len(x)): return []
        return numpy.mean(x, axis=0).tolist()

    @staticmethod
    def _sem(x):
        
        if not(len(x)): return []
        return numpy.sqrt(numpy.var(x, axis=0) / len(x)).tolist()
    
    @staticmethod
    def _cum(x):
        mu = netburst._mean(x)
        if not(mu): return []
        csum = numpy.cumsum(mu, dtype='float')
        return (csum / csum[-1]).tolist()

    def histo(self, delay, bin=1):
        
        _bins = numpy.arange(0, delay + bin, bin)
        def func(x): return numpy.histogram(x, bins = _bins)[0]
        return func
    
    def sortingFunc(self, clust_nr, outl, Mahal, _norm = True):

        _sorting = procflow.sorting(min_clusters_nr = clust_nr, max_clusters_nr = clust_nr,
                                    rmvOutliers = outl, Mahalanobis_threshold = Mahal,
                                    normalize = _norm)
        def doSorting(Data): return _sorting.__clustering__(Data, range(len(Data)))
        return doSorting

    def sortDelays(self, delays_cl0, delays_cl1, delays_cl0Bin, delays_cl1Bin, delaysNN0, delaysNN1):
        """Finding states (or modes) of network bursting activity"""
        
        bin_div = self.bin if divide_per_bin else 1.        
        delays = numpy.sum([delays_cl0, delays_cl1], axis=0)
        delaysBin = numpy.sum([delays_cl0Bin, delays_cl1Bin], axis=0)
        delaysNN = numpy.sum([delaysNN0, delaysNN1], axis=0)
        if len(delays) == 0 or not(numpy.any(delays)):
            print "Attention: no Network Burst detected on step %s!" % self.step_name
            return
        
        ######################################################
        # sortingStates: for each state, a tuple ([NBs indexes],[outliers NBs ind.])  
        ######################################################
        def _sort(n): return self.sortingFunc(clust_nr = n,
                                              outl = self.Net_features_bin['remove_outliers'],
                                              Mahal = self.Net_features_bin['Mahalanobis_threshold'])
        data2csv = [[" "] + list(numpy.arange(0, self.delay + self.Net_features_bin.get('binstates'), 
                                              self.Net_features_bin.get('binstates')))]
        PcaDelays = pca.pca(delays).SCORE[:,:3]
        data2cluster = numpy.empty((self.NB_nr, len(self.Net_features)))
        for i, k in enumerate(self.Net_features):
            if k == 'Del1': data2cluster[:,i] = PcaDelays[:,0]
            if k == 'Del2': data2cluster[:,i] = PcaDelays[:,1]
            if k == 'Del3': data2cluster[:,i] = PcaDelays[:,2]
            if k == 'BD':   data2cluster[:,i] = self.durations
            if k == 'IBI':  data2cluster[:,i] = self.ibis
            if k == 'IBI2': data2cluster[:,i] = numpy.array(list(self.ibis[1:])+[0])
            if k == 'SN':   data2cluster[:,i] = self.spikes_nr[:,0]
            if k == 'NN':   data2cluster[:,i] = self.neurons_nr[:,0]
        
        sortingStates = _sort(self.states_nr)(data2cluster)
        sortingStates = sorted((len(c[0]), c) for c in sortingStates.itervalues())
        _trashNr = self.NB_nr * self.trash_poor_states / 100.
        sortingStates = zip(*sortingStates)[1][::-1]
        sortingStates = [tup for tup in sortingStates if len(tup[0]) > _trashNr]
        self.states_nr = min(self.states_nr, len(sortingStates))
        sortingStates34 = sortingStates[2:4] #if len(sortingStates)>2 else None
        if self.do_Ttest and self.states_nr>1:
            ##################################        
            # T-test on delays of Exc/Inh in the two clusters 
            ##################################
            ix0, ix1 = sortingStates[0][0], sortingStates[1][0]
            delExc0 = netburst._mean(delays_cl0[ix0])[:self.Net_features_bin['nr_pvalue']]
            delExc1 = netburst._mean(delays_cl0[ix1])[:self.Net_features_bin['nr_pvalue']]
            delInh0 = netburst._mean(delays_cl1[ix0])[:self.Net_features_bin['nr_pvalue']]
            delInh1 = netburst._mean(delays_cl1[ix1])[:self.Net_features_bin['nr_pvalue']]
            if scipy.stats.ttest_rel(delExc0, delExc1)[1] > self.Net_features_bin['pvalue'] and \
               scipy.stats.ttest_rel(delInh0, delInh1)[1] > self.Net_features_bin['pvalue']:
                #print "Going down from %d to %d network bursts states (controllo 12)" % (self.states_nr, (self.states_nr-1))
                self.states_nr -= 1
                sortingStates = [(sortingStates[0][0]+sortingStates[1][0],
                                  sortingStates[0][1]+sortingStates[1][1])]
                #print 'len(sortingStates) dopo il going down del controllo 12:', len(sortingStates)
            if len(sortingStates34) == 1: sortingStates += sortingStates34 
            elif len(sortingStates34) == 2:
                ix0, ix1 = sortingStates34[0][0], sortingStates34[1][0]
                delExc0 = netburst._mean(delays_cl0[ix0])[:self.Net_features_bin['nr_pvalue']]
                delExc1 = netburst._mean(delays_cl0[ix1])[:self.Net_features_bin['nr_pvalue']]
                delInh0 = netburst._mean(delays_cl1[ix0])[:self.Net_features_bin['nr_pvalue']]
                delInh1 = netburst._mean(delays_cl1[ix1])[:self.Net_features_bin['nr_pvalue']]
                if scipy.stats.ttest_rel(delExc0, delExc1)[1] > self.Net_features_bin['pvalue'] and \
                   scipy.stats.ttest_rel(delInh0, delInh1)[1] > self.Net_features_bin['pvalue']:
                    #print "Going down from %d to %d network bursts states (controllo 34)" % (self.states_nr, (self.states_nr-1))
                    self.states_nr -= 1
                    #print 'len(sortingStates) dopo il going down  del controllo 34 e prima di aggiornamento:', len(sortingStates)
                    #print 'sortingStates dopo il going down  del controllo 34 e prima di aggiornamento:', sortingStates
                    sortingStates += [(sortingStates34[0][0]+sortingStates34[1][0],
                                       sortingStates34[0][1]+sortingStates34[1][1])]
                    #print 'sortingStates dopo il going down  del controllo 34 e dopo di aggiornamento:', sortingStates
                    #print 'len(sortingStates) dopo il going down del controllo 34:', len(sortingStates)
        self.sortingStates = sortingStates
        self.states_nr = len(sortingStates)
        self.IxStates = IxStates = zip(*sorted((ix,s) for s, ixs in enumerate(sortingStates) for ix in ixs[0]))[1]
        ##################################
        # Prob states passages matrix 
        ##################################
        if store_probStates:
            probStates = numpy.zeros((self.states_nr, self.states_nr))
            itIxStates0, itIxStates1 = itertools.tee(IxStates) ; itIxStates1.next()
            for i, j in zip(itIxStates0, itIxStates1): probStates[i, j] += 1.
            probStates /= len(IxStates)-1
            csvwriterows_safe('probStates_%s.csv' % self.step_name, probStates.tolist(),
                              '%s%s' % (self.data_path, 'nb_clusters\/probStates\/'))

        _all_bins = numpy.arange(0, self.delay + self.bin, self.bin)
        data2csvNN = [[" "] + _all_bins.tolist()]
        data2csvSR = [[" "] + _all_bins.tolist()]
        sortinfo = [" ", self.step_name, str(self.NB_nr)]
        neurons_nr_mean, _ibis = [], self.ibis
        NNmeans = []
        self.infoHeader = sortinfo[1:] + \
                          [str(netburst._mean(_ibis)), str(netburst._sem(_ibis))]
        check_exc, check_inh, self.sortinfo = [], [], []
        
        def cumTuned(x):
            xs = numpy.array(netburst._cum(x))
            if not(len(xs)): return []
            xs -= xs[0]
            return list(xs)
        binsCounts = numpy.arange(0, self.Net_features_bin['maxCounts'] + self.Net_features_bin['binCounts'],
                                  self.Net_features_bin['binCounts'])
        binsCounts = binsCounts.tolist(); binsCounts[0] = 0.5
        hisCountsF = self.histo(self.Net_features_bin['maxCounts'], self.Net_features_bin['binCounts'])
        self.HisCounts, self.CumCounts = [[" "]+binsCounts], [[" "]+binsCounts]
        self.FF = [[" ", 0.001] + \
                    list(numpy.arange(self.Net_features_bin['binFF'], self.delay, self.Net_features_bin['binFF']))]
        def _FanoFactor(xs): return max(0, Fano_factor(xs))
        
        for i, tup in enumerate(sortingStates):
            NBix, OutliersNBix = tup
            if i < 3: sortinfo.extend(map(str, map(len, tup)))
            sortinfoTemp = ["%.2f" % netburst._mean(_ibis[NBix]), 
                            "%.2f" % netburst._sem(_ibis[NBix]),
                            "%.2f" % (100*float(len(NBix)) / self.NB_nr),
                            "%.2f" % (100*float(len(OutliersNBix)) / len(NBix))]
            _d, _d_cl0, _d_cl1 = delays[NBix]/bin_div, delays_cl0[NBix]/bin_div, delays_cl1[NBix]/bin_div
            _dNN, _dNN0, _dNN1 = delaysNN[NBix], delaysNN0[NBix], delaysNN1[NBix]
            _dSR  = [[p / q if q else 0. for p, q in zip(x, y)] for x, y in zip(delaysBin[NBix], delaysNN[NBix])]
            _dSR0 = [[p / q if q else 0. for p, q in zip(x, y)] for x, y in zip(delays_cl0Bin[NBix], delaysNN0[NBix])]
            _dSR1 = [[p / q if q else 0. for p, q in zip(x, y)] for x, y in zip(delays_cl1Bin[NBix], delaysNN1[NBix])]
            data2csv.extend([[" "] + netburst._mean(_d), [" "] + netburst._sem(_d),
                             [" "] + netburst._mean(_d_cl0), [" "] + netburst._sem(_d_cl0),
                             [" "] + netburst._mean(_d_cl1), [" "] + netburst._sem(_d_cl1)])
            data2csvNN.extend([[" "] + netburst._mean(_dNN), [" "] + netburst._sem(_dNN),
                               [" "] + netburst._mean(_dNN0), [" "] + netburst._sem(_dNN0),
                               [" "] + netburst._mean(_dNN1), [" "] + netburst._sem(_dNN1)])
            data2csvSR.extend([[" "] + netburst._mean(_dSR),  [" "] + netburst._sem(_dSR),
                               [" "] + netburst._mean(_dSR0), [" "] + netburst._sem(_dSR0),
                               [" "] + netburst._mean(_dSR1), [" "] + netburst._sem(_dSR1)])
            NN_mean, NN_sem = netburst._mean(self.neurons_nr[NBix]), netburst._sem(self.neurons_nr[NBix])
            SN_mean, SN_sem = netburst._mean(self.spikes_nr[NBix]), netburst._sem(self.spikes_nr[NBix])
            NNmeans.append(NN_mean)
            for j, tup in enumerate(zip(NN_mean, NN_sem)):
                NNval = ["%.1f" % _i for _i in tup]
                neurons_nr_mean.extend(NNval)
                sortinfoTemp.extend(NNval)
                if j == 0: continue
                NNvalperc = ["%.1f" % ((100.*tup[0])/len(self.set_cl[j-1])), \
                                       "%.1f" % ((100.*tup[1])/len(self.set_cl[j-1]))]
                sortinfoTemp.extend(NNvalperc)
            for j, tup in enumerate(zip(SN_mean, SN_sem)):
                SNval = ["%.1f" % _i for _i in tup]
                sortinfoTemp.extend(SNval)
                if j == 0: continue
                SNvalperc = ["%.1f" % ((100.*tup[0])/SN_mean[0]), \
                             "%.1f" %  ((100.*tup[1])/SN_mean[0])]
                sortinfoTemp.extend(SNvalperc)
            self.sortinfo.append(sortinfoTemp)
            if i > 1: continue 
            CountsE = [hisCountsF(xs[NBix])/self.NB_nr for xs in self.probNeurEs.values()] 
            CountsI = [hisCountsF(xs[NBix])/self.NB_nr for xs in self.probNeurIs.values()]
            Counts, CountsE, CountsI = map(numpy.array, (CountsE + CountsI, CountsE, CountsI))
            self.HisCounts.extend([["mu%d" %(i+1)] + netburst._mean(Counts), ["sem"] + netburst._sem(Counts),
                                   ["muExc%d" %(i+1)] + netburst._mean(CountsE), ["sem"] + netburst._sem(CountsE),
                                   ["muInh%d" %(i+1)] + netburst._mean(CountsI), ["sem"] + netburst._sem(CountsI)])
            self.CumCounts.extend([["mu%d" %(i+1)] + cumTuned(Counts), [" "],
                                   ["muExc%d" %(i+1)] + cumTuned(CountsE), [" "],
                                   ["muInh%d" %(i+1)] + cumTuned(CountsI), [" "]])
            if store_FF:
                FFE = [map(_FanoFactor, self.FFcounts[neur][NBix].T) for neur in self.set_cl[0]]
                FFI = [map(_FanoFactor, self.FFcounts[neur][NBix].T) for neur in self.set_cl[1]]
                FF  = FFE+FFI
                self.FF.extend([["mu%d" %(i+1)] + netburst._mean(FF), ["sem"] + netburst._sem(FF),
                                ["muExc%d" %(i+1)] + netburst._mean(FFE), ["sem"] + netburst._sem(FFE),
                                ["muInh%d" %(i+1)] + netburst._mean(FFI), ["sem"] + netburst._sem(FFI)])
        if len(sortingStates) < 3:
            data2csv.extend([" "]*2)
            data2csvNN.extend([" "]*2)
            data2csvSR.extend([" "]*2)
            sortinfo.extend([" "]*2)
            neurons_nr_mean.extend([" "]*2)
        sortinfo.extend(neurons_nr_mean)
        self.infoTail = []
        for k in (8,12,18,22):
            _s = 0
            try: _s += float(self.sortinfo[0][2]) * float(self.sortinfo[0][k])
            except IndexError: pass
            try: _s += float(self.sortinfo[1][2]) * float(self.sortinfo[1][k])
            except IndexError: pass
            self.infoTail.append(_s*.01)
        
        if store_FF:
            self.FF = map(list, map(None, *self.FF))
            self.FF[0][:] = sortinfo[:13]; self.FF = map(None, *self.FF)
        self.HisCounts = map(list, map(None, *self.HisCounts))
        self.HisCounts[0][:] = sortinfo[:13]; self.HisCounts = map(list, map(None, *self.HisCounts))
        self.CumCounts = map(list, map(None, *self.CumCounts))
        self.CumCounts[0][:] = sortinfo[:13]; self.CumCounts = map(list, map(None, *self.CumCounts))
        self.SpikeCounts = [xs[:] for xs in self.HisCounts]
        for i in xrange(1, 12, 2):
            _binsIt = iter(binsCounts)
            for j in xrange(1, len(self.SpikeCounts[i])):
                if self.SpikeCounts[i][j]: self.SpikeCounts[i][j] *= _binsIt.next()
                else: _binsIt.next()    
        
        if len(data2csv) != 1: data2csv[0][1] = .001; data2csv = map(None, *data2csv); data2csv[0] = sortinfo
        else: data2csv = [[x,' '] for x in data2csv[0]]
        if len(data2csvNN) != 1: data2csvNN[0][1] = .001; data2csvNN = map(None, *data2csvNN); data2csvNN[0] = sortinfo
        else: data2csvNN = [[x,' '] for x in data2csvNN[0]]
        if len(data2csvSR) != 1: data2csvSR[0][1] = .001; data2csvSR = map(None, *data2csvSR); data2csvSR[0] = sortinfo
        else: data2csvSR = [[x,' '] for x in data2csvSR[0]]
        
        self.NetStates   = map(None, *data2csv)[:13]    
        self.NetStates34 = map(list, map(None, *([self.NetStates[0]] + map(None, *data2csv)[13:13+12])))
        self.NetStates34[0][:13] = sortinfo[:13]; self.NetStates34[0][3:9] = sortinfo[7:9]+[" "]*4; self.NetStates34 = map(None, *self.NetStates34)
        
        self.NetStatesNN    = map(None, *data2csvNN)[:13]    
        self.NetStatesNN34  = map(list, map(None, *([self.NetStatesNN[0]] + map(None, *data2csvNN)[13:13+12])))    
        self.NetStatesNN34[0][:13] = sortinfo[:13]; self.NetStatesNN34[0][3:9]=sortinfo[7:9]+[" "]*4; self.NetStatesNN34=map(None, *self.NetStatesNN34)
        
        self.NetStatesEX    = map(None, *data2csvSR)[:13]   
        self.NetStatesEX34  = map(list, map(None, *([self.NetStatesEX[0]] + map(None, *data2csvSR)[13:13+12])))    
        self.NetStatesEX34[0][:13] = sortinfo[:13]; self.NetStatesEX34[0][3:9]=sortinfo[7:9]+[" "]*4; self.NetStatesEX34=map(None, *self.NetStatesEX34)
        
        self.NetStatesNNc   = map(list, map(None, *data2csvNN)[:13])
        self.NetStatesNNc34 = map(list, map(None, *([self.NetStatesNNc[0]] + map(None, *data2csvNN)[13:13+12])))    
        self.NetStatesNNc34[0][:13] = sortinfo[:13]; self.NetStatesNNc34[0][3:9]=sortinfo[7:9]+[" "]*4; self.NetStatesNNc34 = map(list, map(None,*self.NetStatesNNc34))
        for j in xrange(1, len(self.NetStatesNN), 2):
            _2cumsum = numpy.cumsum([i for i in self.NetStatesNNc[j][1:] if i != None and i != ''])
            if len(_2cumsum):
                _2cumsum /= _2cumsum[-1] + 1.*(_2cumsum[-1]==0) 
                if j == 3:  _2cumsum *= NNmeans[0][1] / len(self.set_cl[0])
                if j == 5:  _2cumsum *= NNmeans[0][2] / len(self.set_cl[1])
                if j == 9:  _2cumsum *= NNmeans[1][1] / len(self.set_cl[0])
                if j == 11: _2cumsum *= NNmeans[1][2] / len(self.set_cl[1])
            self.NetStatesNNc[j][1:1+len(_2cumsum)] = list(_2cumsum)
            self.NetStatesNNc[j+1][1:] = [''] * (len(self.NetStatesNNc[j+1])-1)
        for j in xrange(1, len(self.NetStatesNN34), 2):
            _2cumsum = numpy.cumsum([i for i in self.NetStatesNNc34[j][1:] if i != None and i != ''])
            if len(_2cumsum):
                _2cumsum /= _2cumsum[-1] + 1.*(_2cumsum[-1]==0) 
                try:
                    if j == 3:  _2cumsum *= NNmeans[2][1] / len(self.set_cl[0])
                    if j == 5:  _2cumsum *= NNmeans[2][2] / len(self.set_cl[1])
                    if j == 9:  _2cumsum *= NNmeans[3][1] / len(self.set_cl[0])
                    if j == 11: _2cumsum *= NNmeans[3][2] / len(self.set_cl[1])
                except IndexError: pass
            self.NetStatesNNc34[j][1:1+len(_2cumsum)] = list(_2cumsum)
            self.NetStatesNNc34[j+1][1:] = [''] * (len(self.NetStatesNNc34[j+1])-1)
        
        ################################
        # Storing NBs data and NBstarts
        ################################
        if store_netstates:
            csvwriterows_safe('netstates_%s.csv' % self.step_name, data2csv, '%s%s' % (self.data_path,'nb_clusters\/'))
            csvwriterows_safe('netstatesNN_%s.csv' % self.step_name, data2csvNN, '%s%s' % (self.data_path,'nb_clusters\/'))
        if store_NBstarts:
            _NBstarts = [self.NBstarts[j] + (i+1,) for i, tup in enumerate(sortingStates) for j in tup[0]]
            _NBstarts.sort()
            csvwriterows_safe('NBstarts_%s.csv' % self.step_name,
                              [['start', 'neuron', 'end', 'state']] + _NBstarts,
                              '%s%s' % (self.data_path, 'nb_clusters\/NBstarts\/'))

        ######################################################
        # Storing NBs features data
        ######################################################
        if self.store_results:
            feat2store = []
            for feat, featDelay, data in (('BD','BD_max',self.durations), ('IBI','IBI_max',self.ibis),
                                          ('IBI2','IBI2_max',numpy.asarray(list(self.ibis[1:])+[0])),
                                          ('NN','NN_max',self.neurons_nr), ('SN', 'SN_max',self.spikes_nr)):
                bin, delay = self.Net_features_bin.get(feat, None), self.Net_features_bin.get(featDelay, None)
                if bin == None or delay == None: continue
                bins = numpy.arange(bin, delay+bin, bin)
                _histo = self.histo(delay, bin)
                _histoData = numpy.array(_histo(data))
                _histoDataNorm = _histoData / (_histoData.sum()+1^(_histoData.sum()!=0))
                _histoCumData = numpy.cumsum(_histoData)
                _histoCumData /= (_histoCumData.sum()+1^(_histoCumData.sum()!=0))
                feat2store.append(['bins'] + list(bins))
                feat2store.append(['%s' % feat] + list(_histoData))
                feat2store.append(['%sn' % feat] + list(_histoDataNorm))
                feat2store.append(['%sc' % feat] + list(_histoCumData))

                for i, tup in enumerate(sortingStates, start=1):
                    NBix, OutliersNBix = tup
                    if i > 4: break
                    _histoData = numpy.array(_histo(data[NBix]))
                    if _histoData.sum() != 0: _histoDataNorm = _histoData / _histoData.sum()
                    _histoCumData = numpy.cumsum(_histoData)
                    if _histoCumData.sum() != 0: _histoCumData /= _histoCumData.sum()
                    
                    feat2store.append(['%s%d'  % (feat, i)] + list(_histoData))
                    feat2store.append(['%sn%d' % (feat, i)] + list(_histoDataNorm))
                    feat2store.append(['%sc%d' % (feat, i)] + list(_histoCumData))
                feat2store.extend([[' ']] * 3 * (5 - len(sortingStates)))

            ######################################################
            # Storing the following Network burst features -> not used anymore, so it is commented
            ######################################################           
            #csvwriterows_safe(self.step_name + '_HIST_netstate_features.csv', map(None, *feat2store),
            #                  '%s%s' % (self.data_path, 'netstate_features\/'))

    def __process__(self, all_ts, **kwargs):

        def divCheck(X, Y): return ((Y!=0) * X) / (Y + 1*(Y==0))             
        
        ######################################################
        # At first NBs detection, then collecting NBs data.
        # delaysNeur: dict of delays (=timestamp-NBstart_t) of each neuron
        ######################################################
        netbr = self.detection(self.start_time, self.end_time, all_ts)
        
        self.NB_nr = float(len(netbr.bursts_container))
        clusters_members = kwargs.get('clusters_members', [(), (), ([], [], []), []])
        self.set_cl = set_cl = set_cl0, set_cl1 = map(set, clusters_members[2][:2])
        set_cl01 = set_cl0 | set_cl1
        neurons_nr = numpy.zeros((self.NB_nr, 3))
        spikes_nr  = neurons_nr.copy()
        durations  = numpy.zeros(self.NB_nr) 
        firstTs = collections.defaultdict(list)
        delaysNeur = firstTs.copy()
        probNeurE  = firstTs.copy(); probNeurI = firstTs.copy()
        self.probNeurEs = dict((neur, numpy.zeros(self.NB_nr)) for neur in set_cl0)  
        self.probNeurIs = dict((neur, numpy.zeros(self.NB_nr)) for neur in set_cl1)  
        NBstarts   = []
        ibis = durations.copy()
        delays_cl0 = numpy.zeros((self.NB_nr, math.ceil(self.delay/self.Net_features_bin.get('binstates'))))
        delays_cl1 = delays_cl0.copy()
        delays_cl0Bin = numpy.zeros((self.NB_nr, math.ceil(self.delay/self.bin)))
        delays_cl1Bin = delays_cl0Bin.copy()
        delaysNN0 = numpy.zeros((self.NB_nr, math.ceil(self.delay/self.bin)))
        delaysNN1 = delaysNN0.copy()
        end_last_burst = 0.
        self.neursSilent = set()
        _mean, _sem = netburst._mean, netburst._sem
        _histo = self.histo(self.delay, self.Net_features_bin.get('binstates'))
        _histoBin = self.histo(self.delay, self.bin)
        binFF = self.Net_features_bin.get('binFF')
        self.FFcounts = dict((neur, numpy.zeros((self.NB_nr, self.delay/binFF))) for neur in set_cl01)

        for i, br in enumerate(netbr.bursts_container):
            neurons_in_burst, ts_in_bursts = zip(*br['spikes_events'])
            _it = iter(br['spikes_events'])
            ts_cl0, ts_cl1 = [], []
            sn_all, sn_exc, sn_inh = 0, 0 ,0
            delaysNeurtemp = collections.defaultdict(list)
            leader, t0 = br['spikes_events'][0]
            ibis[i] = t0-end_last_burst
            if store_FF and leader in set_cl01: self.FFcounts[leader][i][0] += 1
            for neur, t in _it:
                if neur in set_cl01:
                    dt = t-t0
                    delaysNeurtemp[neur].append(dt)
                    if store_FF and dt<self.delay: self.FFcounts[neur][i][dt//binFF] += 1
            durations[i] = t-t0
            NBstarts.append((t0, leader, t))
            end_last_burst = t
            for neur, ts in delaysNeurtemp.iteritems():
                delaysNNtemp = numpy.zeros(math.ceil(self.delay/self.bin), dtype='bool')
                delaysNNtemp[[int(_t / self.bin) for _t in ts if _t < self.delay]] |= 1
                l = len(ts)
                if neur in set_cl0:
                    delaysNeur[neur].extend(ts); ts_cl0.extend(ts)
                    delaysNN0[i] += delaysNNtemp
                    sn_all += l; sn_exc += l
                    probNeurE[neur].append(l); self.probNeurEs[neur][i] += l
                    if ts[0] < delay_netpath: firstTs[neur].append(ts[0])
                else:
                    delaysNeur[neur].extend(ts); ts_cl1.extend(ts) 
                    delaysNN1[i] += delaysNNtemp
                    sn_all += l; sn_inh += l
                    probNeurI[neur].append(l); self.probNeurIs[neur][i] += l
                    if ts[0] < delay_netpath: firstTs[neur].append(ts[0])
            for neur in (set_cl0):
                if neur not in delaysNeurtemp: probNeurE[neur].append(0)
            for neur in (set_cl1):
                if neur not in delaysNeurtemp: probNeurI[neur].append(0)
            delays_cl0[i] = _histo(ts_cl0)[:delays_cl0.shape[1]]
            delays_cl1[i] = _histo(ts_cl1)[:delays_cl1.shape[1]]
            delays_cl0Bin[i] = _histoBin(ts_cl0)[:delays_cl0Bin.shape[1]]
            delays_cl1Bin[i] = _histoBin(ts_cl1)[:delays_cl1Bin.shape[1]]
            neurons_nr[i][:] = map(len, (delaysNeurtemp, \
                            set(delaysNeurtemp) & set_cl0, set(delaysNeurtemp) & set_cl1))
            spikes_nr[i][:] = [sn_all, sn_exc, sn_inh] 
        
        for neur, xs in probNeurE.iteritems():
            if sum(xs) == 0: self.neursSilent.add(neur)
        for neur, xs in probNeurI.iteritems():
            if sum(xs) == 0: self.neursSilent.add(neur)
        self.delaysNeur, self.probNeurE, self.probNeurI = delaysNeur, probNeurE, probNeurI
        ibis[0] = 0 ; self.ibis = ibis
        self.durations  = durations  ; self.NBstarts  = NBstarts
        self.neurons_nr = neurons_nr ; self.spikes_nr = spikes_nr
        
        ###############################################
        # Filling netdata2
        ###############################################
        self.netdata2 = ([_mean(ibis), _sem(ibis),
                          _mean(spikes_nr[:, 1]), _sem(spikes_nr[:, 1]),
                          _mean(spikes_nr[:, 2]), _sem(spikes_nr[:, 2]),
                          _mean(neurons_nr[:, 1]), _sem(neurons_nr[:, 1]),
                          _mean(neurons_nr[:, 2]), _sem(neurons_nr[:, 2]),
                          _mean(divCheck(spikes_nr[:, 1], neurons_nr[:, 1])),
                          _sem(divCheck(spikes_nr[:, 1], neurons_nr[:, 1])),
                          _mean(divCheck(spikes_nr[:, 2], neurons_nr[:, 2])),
                          _sem(divCheck(spikes_nr[:, 2], neurons_nr[:, 2])),
                          _mean(divCheck(spikes_nr[:, 1] * neurons_nr[:, 2], neurons_nr[:, 1] * spikes_nr[:, 2])),
                          _sem(divCheck(spikes_nr[:, 1] * neurons_nr[:, 2], neurons_nr[:, 1] * spikes_nr[:, 2]))])
        
        _netpath = '%s%s' % (self.data_path, 'netpaths\/')
        if store_netpaths:
            isMEA256 = tuple(set_cl0)[0][0] in MEA256columns
            stMatFirst = numpy.zeros((16 << isMEA256, 16 << isMEA256))
            stMatFirstE = stMatFirst.copy() ; stMatFirstI = stMatFirst.copy()
            stMatAll = stMatFirst.copy() ; stMatAllE = stMatFirst.copy() ; stMatAllI = stMatFirst.copy()
            stMatArr, stMatArrA = [], []
            coordsX,  coordsY, coordsXE, coordsYE, coordsXI, coordsYI = [], [], [], [], [], []
            for neur, xs in firstTs.iteritems():
                x, y = mea_coord(neur, augmented=True)[:2]
                muFirst = numpy.mean(xs)
                muAlldt = [dt for dt in delaysNeur[neur] if dt < delay_netpath]
                if muAlldt == []: continue
                muAll = numpy.mean(muAlldt)
                if neur in set_cl0:
                    coordsX.append(x); coordsY.append(y);
                    coordsXE.append(x); coordsYE.append(y)
                    stMatArr.append(muFirst) ; stMatArrA.append(muAll)
                if neur in set_cl1:
                    coordsX.append(x); coordsY.append(y);
                    coordsXI.append(x); coordsYI.append(y)
                    stMatArr.append(muFirst) ; stMatArrA.append(muAll)
            stMatArr, stMatArrA = numpy.asarray(stMatArr), numpy.asarray(stMatArrA)
            stMatArr  /= numpy.amax(stMatArr)+.001  ; stMatArr  -= 1; stMatArr *= -1
            stMatArrA /= numpy.amax(stMatArrA)+.001 ; stMatArrA -= 1; stMatArrA *= -1
            stMatFirst[coordsX, coordsY] = stMatArr ; stMatAll[coordsX, coordsY] = stMatArrA
            stMatFirstE[coordsXE, coordsYE] = stMatFirst[coordsXE, coordsYE]
            stMatFirstI[coordsXI, coordsYI] = stMatFirst[coordsXI, coordsYI]
            stMatAllE[coordsXE, coordsYE] = stMatAll[coordsXE, coordsYE] 
            stMatAllI[coordsXI, coordsYI] = stMatAll[coordsXI, coordsYI]
            if not(normalize_maps):
                stMatFirstE /= numpy.amax(stMatFirstE)
                stMatFirstI /= numpy.amax(stMatFirstI)
                stMatAllE   /= numpy.amax(stMatAllE)
                stMatAllI   /= numpy.amax(stMatAllI)
            
            csvwriterows_safe('f_%s.csv'    % self.step_name, stMatFirst.tolist(),  _netpath)
            csvwriterows_safe('a_%s.csv'    % self.step_name, stMatAll.tolist(),    _netpath)
            csvwriterows_safe('fExc_%s.csv' % self.step_name, stMatFirstE.tolist(), _netpath)
            csvwriterows_safe('aExc_%s.csv' % self.step_name, stMatAllE.tolist(),   _netpath)
            csvwriterows_safe('fInh_%s.csv' % self.step_name, stMatFirstI.tolist(), _netpath)
            csvwriterows_safe('aInh_%s.csv' % self.step_name, stMatAllI.tolist(),   _netpath)
            
        ######################################################
        # Storing features and states info
        # Storing delays data grouped by clusters
        # Storing the following Network burst features 
        ######################################################
        if store_HIST_Delays:
            _Net_data_to_store = [([delaysNeur.get(neur,[]) for neur in cl],'bin','delay') for cl in set_cl]
            Delays_HIST_path = '%s%s' % (self.data_path, 'delays\/')
            if not(os.path.exists(Delays_HIST_path)): os.mkdir(Delays_HIST_path, 0777)   
            self.store_Net_feat_data(_Net_data_to_store,
                                     Delays_HIST_path + self.step_name + '_HIST_Delays.csv', 
                                     to_clusters = True, no_sem = False, **kwargs)
        
        ######################################################
        # Finding Network neurons classes and bursts states
        ######################################################
        self.sortDelays(delays_cl0, delays_cl1, delays_cl0Bin, delays_cl1Bin, delaysNN0, delaysNN1)

    def storeDelaysNeurClasses(self, neurClust = tuple()):
        
        sortingClasses = sorted((len(tup[0]),i,tup[0]) \
                                for i,tup in enumerate(neurClust))[::-1]
        _histoFunc = self.histo(self.delay, self.bin)
        _bins = list(numpy.arange(self.bin, self.delay + self.bin, self.bin))
        data2write, data2writeEI = [], []
        for l, i, neurs in sortingClasses:
            _hData = numpy.array([_histoFunc(self.delaysNeur.get(neur, []))/self.NB_nr for neur in neurs])
            _hDataE = numpy.array([_histoFunc(self.delaysNeur.get(neur, []))/self.NB_nr \
                                    for neur in neurs if neur in self.set_cl[0]])
            _hDataI = numpy.array([_histoFunc(self.delaysNeur.get(neur, []))/self.NB_nr \
                                    for neur in neurs if neur in self.set_cl[1]])
            data2write.extend([["bins"] + _bins, ["mucl%d" % (i+1)] + netburst._mean(_hData),
                               ["sem"] + netburst._sem(_hData), ["cumprob"] + netburst._cum(_hData)])
            data2writeEI.extend([["bins"] + _bins,
                                 ["muExcCl%d" % (i+1)] + netburst._mean(_hDataE),
                                 ["sem"] + netburst._sem(_hDataE), ["cumprob"] + netburst._cum(_hDataE),
                                 ["muInhCl%d" % (i+1)] + netburst._mean(_hDataI),
                                 ["sem"] + netburst._sem(_hDataI), ["cumprob"] + netburst._cum(_hDataI) ])
        data2write[0][0] = data2writeEI[0][0] = self.step_name
        csvwriterows_safe('Delays_%s.csv' % self.step_name, map(None, *data2write),
                          '%s%s' % (self.data_path, 'delaysNeurClasses\/'))
        csvwriterows_safe('EIDelays_%s.csv' % self.step_name, map(None, *data2writeEI),
                          '%s%s' % (self.data_path, 'delaysNeurClasses\/'))

    def storeCountsNeurClasses(self, neurClust = tuple()):
        
        def cumTuned(x):
            xs = numpy.array(netburst._cum(x))
            if not(len(xs)): return []
            xs -= xs[0]
            return list(xs)
        
        sortingClasses = sorted((len(tup[0]),i,tup[0]) for i,tup in enumerate(neurClust))[::-1]
        probNeur = self.probNeurE ; probNeur.update(self.probNeurI)
        maxSpikesNr = self.Net_features_bin['maxCounts'] #max(map(max, probNeur.itervalues()))
        #_bins = numpy.arange(0, max(maxSpikesNr, self.Net_features_bin['maxCounts']) + \
        _bins = numpy.arange(0, maxSpikesNr + self.Net_features_bin['binCounts'],
                             self.Net_features_bin['binCounts'])
        _histoFunc = self.histo(max(maxSpikesNr, self.Net_features_bin['maxCounts']), self.Net_features_bin['binCounts'])
        _bins = _bins.tolist(); _bins[0] = 0.5
        data2write, data2writeEI, Exc_only, Inh_only = [], [], [], []
        for l, i, neurs in sortingClasses:
            _hData, _hDataE, _hDataI = [], [], []
            for neur in neurs:
                _hDataNeur = _histoFunc(probNeur.get(neur, []))/self.NB_nr
                if neur in self.set_cl[0]: _hData.append(_hDataNeur); _hDataE.append(_hDataNeur); Exc_only.append(_hDataNeur)
                if neur in self.set_cl[1]: _hData.append(_hDataNeur); _hDataI.append(_hDataNeur); Inh_only.append(_hDataNeur)
            _hData, _hDataE, _hDataI = map(numpy.array, (_hData, _hDataE, _hDataI))
            data2write.extend([["bins"]+_bins, ["muCl%d" % (i+1)] + netburst._mean(_hData),
                               ["sem"] + netburst._sem(_hData), ["cumprob"] + cumTuned(_hData)])
            data2writeEI.extend([["bins"]+_bins,
                                 ["muExcCl%d" % (i+1)] + netburst._mean(_hDataE),
                                 ["sem"] + netburst._sem(_hDataE), ["cumprob"] + cumTuned(_hDataE),
                                 ["muInhCl%d" % (i+1)] + netburst._mean(_hDataI),
                                 ["sem"] + netburst._sem(_hDataI), ["cumprob"] + cumTuned(_hDataI)])
        Exc_only, Inh_only = map(numpy.array, (Exc_only, Inh_only))
        self.EI_only = [["bins"]+_bins, ["muExc"] + netburst._mean(Exc_only),
                        ["sem"] + netburst._sem(Exc_only), ["cumprob"] + cumTuned(Exc_only),
                        ["bins"]+_bins, ["muInh"] + netburst._mean(Inh_only),
                        ["sem"] + netburst._sem(Inh_only), ["cumprob"] + cumTuned(Inh_only)]
        data2write[0][0] = data2writeEI[0][0] = self.EI_only[0][0] = self.step_name
        self.EI_only_xbins = [xs[:] for xs in self.EI_only]
        _binsIt = iter(_bins)
        for i in xrange(1, len(self.EI_only_xbins[1])): self.EI_only_xbins[1][i] *= _binsIt.next()
        self.EI_only_xbins[1][41:] = smooth(self.EI_only_xbins[1][41:], window_len=5, window='flat')
        self.EI_only_xbins[1][1] = 0
        _binsIt = iter(_bins)
        for i in xrange(1, len(self.EI_only_xbins[5])): self.EI_only_xbins[5][i] *= _binsIt.next()
        self.EI_only_xbins[5][41:] = smooth(self.EI_only_xbins[5][41:], window_len=5, window='flat')
        self.EI_only_xbins[5][1] = 0
        if store_countshistoNeurClasses:
            csvwriterows_safe('counts_%s.csv' % self.step_name, map(None, *data2write),
                              '%s%s' % (self.data_path, 'countshistoNeurClasses\/'))
            csvwriterows_safe('EIcounts_%s.csv' % self.step_name, map(None, *data2writeEI),
                              '%s%s' % (self.data_path, 'countshistoNeurClasses\/'))
            csvwriterows_safe('clustersEI_%s.csv' % self.step_name,
                              map(None, *self.EI_only),
                              '%s%s' % (self.data_path, 'countshistoNeurClasses\/'))
            csvwriterows_safe('clustersEI_xbins_%s.csv' % self.step_name,
                              map(None, *self.EI_only_xbins),
                              '%s%s' % (self.data_path, 'countshistoNeurClasses\/'))
    
    def store_Net_feat_data(self,
                            _Net_feat_data_tuples,
                            filename,
                            no_sem = True,
                            to_clusters = False,
                            no_histo = False,
                            **kwargs):
        """In input:
            _Net_feat_data_tuples is a tuple of tuples, each one giving:
                - the data of a feature;
                - the name identifying the feature, and later used to get
                  from self.Net_features_bin the feature bin width of an histogram;
                - the name of the max value of the feature, used to set the
                  max value of the range of the feature data in the histogram.
            - the filename or the path where to store the histogram file;
            - parameter no_sem: if True the SEM column won't be written;
            - to_clusters defining if we are grouping histogram data
              by clusters;
            - no_histo defining if we want an histogram of data in input;
            - step_name, the name of the experimental step.
        """
        
        ######################################################
        # Filling data_to_write of data from dictionaries
        # passed in input in order to write them soon after.
        # Counter _it is useful only if to_clusters is True.
        # _Net_feat_data_tuples are usually nothing else
        # than self._Net_data_to_store. 
        ######################################################
        step_name, start_time, end_time = self.step_name, self.start_time, self.end_time
        data_to_write, _it = [], itertools.count(1)
        for _d_values, k, b in _Net_feat_data_tuples:

            ######################################################
            # Checking if d is a dict or a list
            # Checking if _d_values is empty and if it's 1-D array
            # Now building histogram of _Net_feat_data_tuples in input..
            ######################################################
            _bin   = self.Net_features_bin.get(k, .005)
            _delay = max(deep_chain_iter(_d_values)) if not(b) else self.Net_features_bin.get(b, 0.3)
            _all_bins = numpy.arange(0, _delay + _bin, _bin)
            _histo = self.histo(_delay, _bin)
            if no_histo: data_histo = numpy.asarray(_d_values)
            else:
                if not(_d_values): data_histo = []; print ' No Network bursts detected on step %s!' % step_name
                elif hasattr(_d_values[0], '__iter__'): data_histo = map(_histo, _d_values)
                else: data_histo = [_histo(_d_values)]

            ######################################################
            # Normalization parameters 
            # Specifically set for results in *_HIST_delays.csv files.
            ######################################################
            bin_div, time_norm, ens_norm, all_norm = False, False, True, True
            if k == 'bin':
                if divide_per_bin: bin_div = _bin
                ens_norm, all_norm, time_norm = False, False, ('bursts_nr', self.NB_nr)

            ######################################################
            # Grouping data by clusters, if requested
            # Passing and filling data_to_write by reference
            ######################################################
            if to_clusters:
                feat_name = ''.join([k + ' cluster ' + str(_it.next())])
                pd_cluster_data(data_to_write,
                                data_histo, 
                                _all_bins,
                                data_name = feat_name,
                                no_sem = no_sem,
                                bin_divide = bin_div,
                                ensemble_normalization = ens_norm,
                                time_normalization = time_norm,
                                all_normalization = all_norm)
            else: data_to_write = [_all_bins.tolist()]+[[n]+h for n,h in enumerate(data_histo)]
        
        ######################################################
        # Inserting the name of the experimental step
        ######################################################
        if step_name: data_to_write[6][1] = step_name

        ######################################################
        # Storing data_to_write on file
        ######################################################
        csvwriterows_safe(filename, map(None, *data_to_write))


