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


from gettimestamps import pneuro
import burstdetection
import burst_clusterizer
import netburst_processing
from storage_functions import burst_data_storage

from importlibs import *


__author__ = "Andrea Maffezzoli"
__revision__ = "$Revision: 11.62 - 28 Jun 2011$"
__name__ = "timeneuro_mod"
__all__ = ["timeneuro_processing", "burstsbox"]



class burstsbox(dict):
    """This class stores data of all burst from each neuron 
    in a dictionary with experimental steps as keys"""

    def __init__(self): pass


class timeneuro_processing(object):
    """Class performing all operations to get and save results
    about single neuron burst detection, obtained by the class
    'maxinterval' in module 'burstdetection'"""
    
    nr_of_inst      = 0
    processing_flow = []
    data            = None
    exp_steps       = None
    store_results   = False
    randomization   = False
    Netburst_detect = globals().get('Netburst_detection', True)
    store_single_neuron_bursts = False
    verbose         = False
    _blk_col_nr     = 6
    selection       = set()
   
    def __init__(self, path = None, name = None):
        """Saving name and path of the experiment to process
        Creating a list where to save steps processing times
        """
        
        self.__class__.nr_of_inst += 1
        self.exp_path = filter(bool, (path, '.'))[0]
        self.exp_name = filter(bool, (name, '%s_run_%s' % (self.__class__.__name__, self.__class__.nr_of_inst)))[0]
        if self.data:      self.get_data = doNothing
        if self.exp_steps: self.return_exp_steps = doNothing
        self.processing_flow = [self.set_paths, self.get_data,
                                self.return_exp_steps, self.set_ref_steps,
                                self.pre_burst_detection, self.clustering, self.burst_detection]
        
        ############################################################################
        # NB Hold the order of the following operations!
        ############################################################################
        if self.store_results:   self.processing_flow.append(self.write_results)
        if self.randomization:   self.processing_flow.append(self.randomize)
        if self.Netburst_detect: self.processing_flow.append(self.network_burst)
        
        self.processing_flow.append(self.print_proc_flow_times)
        self.proc_flow_times = []; self.print_devices = []; self.devices_open = False 
        if self.store_single_neuron_bursts: self.processing_flow.append(self.write_burst_results)
        
        self.BreakOut = False

    def __call__(self):
        
        if self.BreakOut: return
        for func in self.processing_flow: func()

    def pre_burst_detection(self, **kw):
        
        ############################################################################
        # PERFORM BURST DETECTION a' la NEUROEXPLORER MAX-INTERVAL
        # INSTANTIATION OF BURST'S DATA OBJECTS 
        # Initializing burst detection data containers. Perforimg burst detection
        # and storing in memory bursts data for each reference experimental time
        # interval. Initializing burst detection data containers and then storing
        # burst detection data in bbox, bn_to_cluster and burst_data_2csv.
        ############################################################################
        maxint_root                = burstdetection.BurstDetection_root
        maxint_root.results_path   = self.results_path
        maxint_root.exp_name       = self.exp_name
        maxint_root.burst_features = burst_features
        maxint_root.exp_steps_dict = {}
        maxint_root.store_results  = self.store_results
        
        self.bbox = burstsbox()
        for _i in self.ref_steps:
            print 'Performing burst detection on reference', _i, 'experimental step..'
            for i, j, k in self.exp_steps:
                if i == _i: break
            self.bbox[i] = burstdetection.maxinterval(j, k, self.data, step_name = i, **kw)
        
    def burst_detection(self):
        
        t_start = time.clock()
        for i,j,k in self.exp_steps:
            if i in self.bbox: continue
            print 'Performing burst detection on', i, 'experimental step..'
            self.bbox[i] = burstdetection.maxinterval(j, k, self.data, step_name = i, **self.last_feat_kw)

        ############################################################################
        # SELECTION OF NEURONS SPIKING ON ALL REF.STEPS AND COLLECTING THEIR DATA
        ############################################################################
        self.bursts_data_neurons = list(set(deep_chain_iter(self.bbox[_step].bursts_all_data.iterkeys() \
                                       for _step in self.ref_steps)))
        self.proc_flow_times.append(('Single neuron burst detection', time.clock() - t_start))

    
    def clustering(self):

        def printhalflife_ACF(D):        
            try:
                clsData = zip(*sorted((len(c), numpy.array(c)) for c in D.itervalues()))[1]
                ix = burst_features.index("halflife_ACF")
            except: print "Some error occurred when printing Half Life ACF values" ; return
            ixs = [(step, ix + i*len(burst_features)) for i, step in enumerate(self.ref_steps)]
            for step, ix in ixs:
                for j, c in enumerate(clsData):
                    print "halflife_ACF at step %s in cluster %d: %.3f" % (step, j+1, numpy.mean(c[:, ix]))

        def printfeat(D):        
            try:
                clsData = zip(*sorted((len(c), numpy.array(c)) for c in D.itervalues()))[1]
            except: print "Some error occurred when printing Half Life ACF values" ; return
            for ix, feat in enumerate(burst_features):
                ixs = [(step, ix + i*len(burst_features)) for i, step in enumerate(self.ref_steps)]
                for step, ix in ixs:
                    for j, c in enumerate(clsData):
                        print "%s at step %s in cluster %d: %.3f" % (feat, step, j+1, numpy.mean(c[:, ix]))
        
        def prerun0(_bn2Clust):
            clustered_data = clustcls(_bn2Clust, results_path = self.results_path,
                                      exp_name = self.exp_name, exp_step = self.ref_steps, **csv_options)
            clustered_data(True, 2) ; printfeat(clustered_data.clusters_data)
            clustered_data(True, 3) ; printfeat(clustered_data.clusters_data)
        
        def prerun():
            _bn2Clust = burstdetection.maxinterval()
            for step in self.ref_steps: _bn2Clust.add(self.bbox[step])
            prerun0(_bn2Clust)
            while True:
                _repeatClust = askparam()
                if not(_repeatClust): return
                self.pre_burst_detection(**self.last_feat_kw)
                del _bn2Clust ; _bn2Clust = burstdetection.maxinterval()
                for step in self.ref_steps: _bn2Clust.add(self.bbox[step])
                print "ACF_delay=%.2f, Fano_factor_bin=%.2f, Mahalanobis_threshold=%.2f" \
                        % (self.last_feat_kw['ACF_delay'], self.last_feat_kw['Fano_factor_bin'], \
                           clustcls.kmeans_dict['Mahalanobis_threshold'])
                prerun0(_bn2Clust)
                toBreak = raw_input("To stop pre-run digit 'Y', otherwise only Enter: ")
                if toBreak == 'Y': return

        def askparam():
            _repeatClust = False
            try:
                new_ACFd = float(raw_input("To change ACF_delay digit its new value, if not only Enter: "))
                self.last_feat_kw['ACF_delay'] = new_ACFd; _repeatClust = True
            except: pass
            try:
                new_FFbin = float(raw_input("To change Fano_factor_bin digit its new value, if not only Enter: "))
                self.last_feat_kw['Fano_factor_bin'] = new_FFbin; _repeatClust = True
            except: pass
            try:
                new_Mt = float(raw_input("To change Mahalanobis_threshold digit its new value, if not only Enter: "))
                self.last_feat_kw['Mahalanobis_threshold'] = \
                        clustcls.kmeans_dict['Mahalanobis_threshold'] = new_Mt; _repeatClust = True
            except: pass
            return _repeatClust
        
        clustcls = burst_clusterizer.clusterize
        self.last_feat_kw = {'ACF_delay': ACF_delay, 'Fano_factor_bin': Fano_factor_bin,
                             'Mahalanobis_threshold': Mahalanobis_threshold} 
        if 'Release' in globals():
            do_prerun = raw_input("To perform prerun, digit 'Y' then Enter, otherwise only Enter: ")
            if do_prerun == 'Y':
                for s in ("ACF_delay", "Fano_factor_bin", "Mahalanobis_threshold"):
                    print 'Initial value for %s: %.2f' % (s, globals()[s])
                prerun()
                clustcls.kmeans_dict['Mahalanobis_threshold'] = globals()['Mahalanobis_threshold']
                new_ref_steps = raw_input("To change reference step(s) digit Y, if not only Enter: ")
                if new_ref_steps == 'Y':
                    self.set_ref_steps(verbose = False) ; self.pre_burst_detection(**self.last_feat_kw)
        
        bn2Clust = burstdetection.maxinterval()
        for step in self.ref_steps: bn2Clust.add(self.bbox[step])
        clustered_data = clustcls(bn2Clust, results_path = self.results_path,
                                  exp_name = self.exp_name, exp_step = self.ref_steps, **csv_options)
        clustered_data()
        self.clustered_data = clustered_data
                
        ############################################################################
        # self.clusters_names_sorted_x_nr contains:
        # - cluster ordinal index before sorting;
        # - neurons nr in a cluster, used to sort the clusters;
        # - neurons name/id in a cluster;
        # - cluster ordinal index after sorting, it is like a range(...). 
        # self.clusters_nr: number of clusters. 
        ############################################################################
        clusters_sorted_x_nr = sorted(zip(xrange(self.clustered_data.clusters_nr), \
                                          self.clustered_data.nr_samples_for_cluster, \
                                          self.clustered_data.clusters_names.itervalues()), \
                                      key = operator.itemgetter(1), reverse = True)
        while True:
            trashthreshold = raw_input("Set the max nr of neurons a cluster should have to be discarded: ")
            if trashthreshold == '': trashthreshold = 0; break
            try: trashthreshold = int(trashthreshold); break
            except ValueError: pass
        clusters_sorted_x_nr = [tup for tup in clusters_sorted_x_nr if tup[1] > trashthreshold]
        if trashthreshold:
            print "After applying threshold above, final clustering results are:"
            for i, tup in enumerate(clusters_sorted_x_nr): print "Cluster %d contains %d neurons" % (i+1, tup[1])
        self.clusters_nr = len(clusters_sorted_x_nr)
        self.clusters_names_sorted_x_nr = zip(*clusters_sorted_x_nr) + [range(1, len(clusters_sorted_x_nr)+1)]

    def set_paths(self):

        ############################################################################
        # PATH INSTANTIATION
        # Setting paths for all results and other data, and creating if necessary
        # the directories where to store results files
        ############################################################################
        self.root_results_path  = self.exp_path + '\/MEAdata_cluster_results\/'
        self.results_path       = self.root_results_path + 'single_neuron_burst_results\/'
        self.feat_storing_path  = self.results_path + 'bursts_data_extracted\/'
        self.Net_path           = '%s%s' % (self.root_results_path, 'Network_burst_results\/')
        for path in (self.root_results_path, self.results_path, self.feat_storing_path, self.Net_path):
            if not(os.path.exists(path)): os.mkdir(path, 0777)
        self.csv_path           = '%sclusters_data' % self.results_path
        self.settings_path      = '%s%s%s%s' % (self.exp_path, '\/MEAdata_settings\/', self.exp_name, '.txt')
    
    def get_data(self):
        
        print "Experiment name: %s" % self.exp_name
        while True:
            minspikesNr = raw_input("Set minimum spikes number per neuron for the entire record or Enter if none: ")
            try:
                self.minspikes_nr = 0 if minspikesNr == '' else int(minspikesNr)
                break
            except ValueError: pass
        try:
            t_start = time.clock()
            
            self.pneuro = pneuro(self.exp_name, '%s\/MEAdata_plx\/' % self.exp_path,
                                 min_spikes_nr = self.minspikes_nr)
            self.data = self.pneuro.transplall
            if self.selection:
                self.data = [row for row in self.data if row[0] in self.selection]
            
            if neurons_excluded:
                self.data = [row for row in self.data if row[0] not in neurons_excluded]

            self.proc_flow_times.append(('Timestamps data retrieving', time.clock() - t_start))
        except IOError, (errno, strerror):
            print 'File csv of experimental timestamps is absent or corrupted'
            print 'Close the program, correct such issue if possible, then re-launch it'
            raise IOError

    def return_exp_steps(self):
        """Reading the experimental timestamps from *.txt file placed in
        settings_path, and collecting them in the list 'exp_steps',
        of tuples ('step name', start time, end time).
        """

        try:
            self.exp_steps, self.BreakOut = [], True
            with open(self.settings_path, 'r') as f:
                for s in f:
                    if s.startswith('#') or s.count('\t') != 2: continue
                    try:
                        _t = s.translate(None, '\n\r').split('\t')
                        _t = tuple([_t[0]] + map(float, _t[1:3]))
                        if _t[1] >= self.pneuro.maxt:
                            print "Interval %s [%.1fs,%.1fs] is empty (data end at %.1fs)" % (_t + (self.pneuro.maxt, )) 
                            continue
                        if len(_t) == 3: self.exp_steps.append(_t); self.BreakOut = False
                    except ValueError: pass
        except IOError, (errno, strerror):
            print 'Impossible to find the file: %s' % self.settings_path
            print "I/O error(%s): %s" % (errno, strerror)
            raise
        except:
            print 'Please check the format of file: %s' % self.settings_path
            print "Unexpected error:", sys.exc_info()[0]
            raise

    def set_ref_steps(self, verbose = True):
        
        steps_names = set(s[0] for s in self.exp_steps)
        if verbose:
            print "The following experimental intervals will be processed:"
            for ix, t in enumerate(self.exp_steps, start = 1):
                i,j,k = t
                print "[%d] %s:\t%s-%s" % (ix, i.ljust(1), repr(j).rjust(7), repr(k).ljust(8))
        while True:
            ref = raw_input("Digit index of clustering reference interval(s) (comma separated if many): ")
            if not(ref):
                if 'Release' not in globals(): self.ref_steps = [self.exp_steps[0][0]]; break
                continue
            self.ref_steps = []
            for i in ref.split(','):
                s =  self.exp_steps[int(i)-1][0]
                if s not in steps_names: continue
                self.ref_steps.append(s)
            if not(self.ref_steps): continue
            break

    def io_print_flow_devices(self, s = None, tup = None, closeMe = False):
        
        if closeMe:
            if getattr(self, 'f_flow', False):
                self.f_flow.write('\n')
                if getattr(self, 'f_flow_olddata', False):
                    self.f_flow.writelines(self.f_flow_olddata)
            for dev in self.print_devices:
                if dev != sys.stdout:
                    try: dev.close()
                    except: pass
            self.devices_open = False
            return
        if not(self.devices_open):
            if self.verbose: self.print_devices.append(sys.stdout)
            if 'Release' in globals() or self.verbose:
                f_flowFile = '%s%s__analysis_times.txt' % (self.root_results_path, self.exp_name)
                if not(os.path.exists(f_flowFile)): self.f_flow = open(f_flowFile, 'w')
                else:
                    self.f_flow = open(f_flowFile, 'r+')
                    self.f_flow_olddata = self.f_flow.readlines()
                    self.f_flow.seek(0)
                self.f_flow.write("%s%s" % ('-' * 70, '\n'))
                self.print_devices.append(self.f_flow)
            self.devices_open = True
        for dev in self.print_devices:
            if tup: print >> dev, s % tup
            else: print >> dev, s

    def print_proc_flow_times(self):

        ############################################################################
        # PRINT ON GIVEN DEVICES TIMES OF PROCESSING FLOW 
        ############################################################################
        datetimeNow = datetime.datetime.now; tupNow = datetimeNow().timetuple()
        tupNowtostr = tuple(str(z).zfill(2) for z in tupNow[:3][::-1] + tupNow[3:6])
        self.io_print_flow_devices("Date: %s/%s/%s  h: %s:%s:%s", tupNowtostr)
        self.io_print_flow_devices("Amount of time occurred for each processing step:")
        tot_time = 0; _print_flow = self.io_print_flow_devices
        for tup in self.proc_flow_times: _print_flow(" %s: %.3fs", tup); tot_time += tup[1]
        _print_flow("Total amount of time occurred:  %.3fs", (tot_time, ))
        _print_flow(closeMe = True)

    def randomize(self, write_randomized_results = False):
        
        ############################################################################
        # CALL RANDOMIZATION PROCEDURE (OPTIONAL)
        ############################################################################
        t_start = time.clock()

        ############################################################################
        # Storing original info and changing paths for randomized results 
        ############################################################################
        info_before_random = (self.results_path, self.feat_storing_path, self.exp_name, self.csv_path)
        self.results_path = '%s%s' % (self.results_path[:self.results_path[:-1].rfind('/') + 1], \
                                      'random_timeneuro_processing_results\/')
        self.feat_storing_path = self.results_path + 'bursts_data_extracted\/'
        for path in (self.results_path, self.feat_storing_path):
            if not(os.path.exists(path)): os.mkdir(path, 0777)
        self.exp_name += '__random_'
        
        ############################################################################
        # Now randomizing burst data in self.bbox
        ############################################################################
        for bn in self.bbox.itervalues(): bn.random_burst()
        
        ############################################################################
        # Now storing all results of previous processing
        ############################################################################
        if self.store_results:
            print 'Writing results from randomized burst data..'
            self.write_results()
            self.write_burst_results()
        
        ############################################################################
        # Finally rolling back to the original paths
        ############################################################################
        self.results_path, self.feat_storing_path, self.exp_name, self.csv_path = info_before_random
        
        self.proc_flow_times.append(('Burst randomization', time.clock() - t_start))

    def network_burst(self, *args, **kwargs):
    
        ############################################################################
        # Storing tables of % of Exc/Inh neurons respect classes from networks activity 
        ############################################################################
        def storeCmpTable(neurClust, set_cl, path):

            sortingClasses = sorted((len(tup[0]),i,tup[0]) for i,tup in enumerate(neurClust))[::-1]
            tableCmp = [[" ", "Exc", "Inh"]]
            for l, i, neurs in sortingClasses:
                tableCmp.append(["Class %d" % (i+1), \
                                100.*len(set(neurs) & set_cl[0])/len(set_cl[0]), \
                                100.*len(set(neurs) & set_cl[1])/len(set_cl[1])] + \
                                [" ", " ", "NeurNr:", l, "Neurons:", neurs])
            csvwriterows_safe('table.csv', tableCmp, '%s%s' % (self.Net_path, path))
        
        t_start = time.clock()
        npn = netburst_processing.netburst
        npn.data_path, npn.store_results = self.Net_path, self.store_results
        npn.trash_poor_states = self.trash_poor_states = trash_poor_states
        npn.Net_features_bin = self.Net_features_bin = Net_features_bin
        npn.Net_features = self.Net_features = Net_features
        while True:
            ############################################################################
            # CALL NETWORK BURSTS PROCEDURE (OPTIONAL)
            # If choosen to perform Network burst detection on clusters previously found
            ############################################################################
            sortingInfos = [['step','NBnr','IBImean','IBIsem']]
            for i in xrange(1, self.Net_features_bin['states_nr']+1):
                sortingInfos[0].extend(["IBI%d"    % i, "sem",
                                        "NBnr%d%c" % (i, '%'), "Outnr%d%c" % (i, '%'),
                                        "NN%d"     % i, "sem",
                                        "ExcNN%d"  % i, "sem", "ExcNN%d%s" % (i, '%'), "sem",
                                        "InhNN%d"  % i, "sem", "InhNN%d%s"% (i, '%'), "sem",
                                        "SN%d"     % i, "sem",
                                        "ExcSN%d"  % i, "sem", "ExcSN%d%s" % (i, '%'), "sem",
                                        "InhSN%d"  % i, "sem", "InhSN%d%s"% (i, '%'), "sem"])
            sortingInfos[0].extend(["nrofSilent", "muNNExc%", "muNNInh%",
                                    "muSNExc%", "muSNExc%", "asynCumProbExc", "asynCumProbInh"]) 
            dict2clustDelays = collections.defaultdict(list)
            dict2clustCounts = collections.defaultdict(list)
            print "Collecting network burst data to find neurons classes.."
            for _step in self.exp_steps:
                if _step[0] not in self.ref_steps: continue
                netbrs = npn(_step_info = _step)
                netbrs.do_Ttest = False
                netbrs(self.data, clusters_members = self.clusters_names_sorted_x_nr)
                for neur, xs in netbrs.delaysNeur.iteritems(): dict2clustDelays[neur].extend(xs)
                for neur, xs in netbrs.probNeurE.iteritems(): dict2clustCounts[neur].extend(xs)
                for neur, xs in netbrs.probNeurI.iteritems(): dict2clustCounts[neur].extend(xs)
            
            ############################################################################
            # Cluster neurons histograms counts from ref steps
            ############################################################################
            print "Finding classes based on neurons histograms spike counts"
            netbrs.do_Ttest = False
            _histoFunc = netbrs.histo(max(map(max, dict2clustCounts.itervalues())), 1)
            _sort = netbrs.sortingFunc(clust_nr = self.Net_features_bin['classesCounts_nr'],
                                       outl = self.Net_features_bin['rm_outliers_classesCounts'],
                                       Mahal = self.Net_features_bin['Mahalanobis_thresh_classesCounts'])
            _Neurs = numpy.array(dict2clustCounts.keys())
            countsClasses = tuple((list(_Neurs[x]), list(_Neurs[y])) for x, y in \
                                  _sort(map(_histoFunc, dict2clustCounts.itervalues())).itervalues())
            storeCmpTable(countsClasses, netbrs.set_cl, 'countshistoNeurClasses\/')
            
            ############################################################################
            # Cluster neurons delays from ref steps
            ############################################################################
            if store_delaysNeurClasses:
                print "Finding classes based on neurons spike delays"
                netbrs.do_Ttest = False
                _histoFunc = netbrs.histo()
                _sort = netbrs.sortingFunc(clust_nr = self.Net_features_bin['classesDelays_nr'],
                                           outl = self.Net_features_bin['rm_outliers_classesDelays'],
                                           Mahal = self.Net_features_bin['Mahalanobis_thresh_classesDelays'])
                _Neurs = numpy.array(dict2clustDelays.keys())
                delaysClasses = tuple((list(_Neurs[x]), list(_Neurs[y])) for x, y in \
                                      _sort(map(_histoFunc, dict2clustDelays.itervalues())).itervalues())
                storeCmpTable(delaysClasses, netbrs.set_cl, 'delaysNeurClasses\/')
            
            ############################################################################
            # Cluster neurons histograms counts from ref steps
            ############################################################################
            netdata2ref, netdata2, ixSR, = [], [], desc_burst_features.index('SR')
            neursSilent = []
            NetStates, NetStatesNN, NetStatesNNc, NetStatesEX, clustersEI, clustersEI_xbins = [], [], [], [], [], []
            NetStates34, NetStatesNN34, NetStatesNNc34, NetStatesEX34 = [], [], [], []
            HisCounts, CumCounts, SpikeCounts, FF = [], [], [], []
            for i, _step in enumerate(self.exp_steps, start=1):
                print 'Performing network burst analysis on step', _step[0]
                netbrs = npn(_step_info = _step)
                netbrs(self.data, clusters_members = self.clusters_names_sorted_x_nr)
                netbrs.storeCountsNeurClasses(countsClasses)
                if store_delaysNeurClasses:
                    netbrs.storeDelaysNeurClasses(delaysClasses)
                neursSilent.append([_step[0],":",len(netbrs.neursSilent)," "]+list(netbrs.neursSilent))
                sortingInfosTemp = netbrs.infoHeader
                for k in netbrs.sortinfo: sortingInfosTemp.extend(k)
                sortingInfosTemp.extend([" "]*24*(netbrs.states_nr-len(netbrs.sortinfo))) 
                sortingInfosTemp.append(len(netbrs.neursSilent))
                sortingInfosTemp.extend(netbrs.infoTail)
                sortingInfosTemp.extend([netbrs.EI_only[3][-1], netbrs.EI_only[7][-1]])
                sortingInfos.append(sortingInfosTemp)
                clustersEI.extend(netbrs.EI_only)
                clustersEI_xbins.extend(netbrs.EI_only_xbins)
                set_cl0, set_cl1 = self.clusters_names_sorted_x_nr[2][:2]
                SRE = filter(lambda i: i != None, [self.bbox[_step[0]].desc_bursts_stat[neur][ixSR] for neur in set_cl0])
                SRI = filter(lambda i: i != None, [self.bbox[_step[0]].desc_bursts_stat[neur][ixSR] for neur in set_cl1])
                netdata2Temp = netbrs.netdata2 + [netbrs._mean(SRE), netbrs._sem(SRE), netbrs._mean(SRI), netbrs._sem(SRI)]
                netdata2.append(netdata2Temp)
                if _step[0] in self.ref_steps: netdata2ref.append(netdata2Temp)
                if (i % 15) == 1:
                    NetStates.extend(netbrs.NetStates)
                    NetStatesNN.extend(netbrs.NetStatesNN)
                    NetStatesNNc.extend(netbrs.NetStatesNNc)
                    NetStatesEX.extend(netbrs.NetStatesEX)
                    if Net_features_bin['states_nr'] > 2:
                        NetStates34.extend(netbrs.NetStates34)
                        NetStatesNN34.extend(netbrs.NetStatesNN34)
                        NetStatesNNc34.extend(netbrs.NetStatesNNc34)
                        NetStatesEX34.extend(netbrs.NetStatesEX34)
                    if store_NeuronsCounts:
                        HisCounts.extend(netbrs.HisCounts)
                        CumCounts.extend(netbrs.CumCounts)
                        SpikeCounts.extend(netbrs.SpikeCounts)
                    if store_FF: FF.extend(netbrs.FF)
                else:
                    NetStates.extend(netbrs.NetStates[1:])
                    NetStatesNN.extend(netbrs.NetStatesNN[1:])
                    NetStatesNNc.extend(netbrs.NetStatesNNc[1:])
                    NetStatesEX.extend(netbrs.NetStatesEX[1:])
                    if Net_features_bin['states_nr'] > 2:
                        NetStates34.extend(netbrs.NetStates34[1:])
                        NetStatesNN34.extend(netbrs.NetStatesNN34[1:])
                        NetStatesNNc34.extend(netbrs.NetStatesNNc34[1:])
                        NetStatesEX34.extend(netbrs.NetStatesEX34[1:])
                    if store_NeuronsCounts:
                        HisCounts.extend(netbrs.HisCounts[1:])
                        CumCounts.extend(netbrs.CumCounts[1:])
                        SpikeCounts.extend(netbrs.SpikeCounts[1:])
                    if store_FF: FF.extend(netbrs.FF[1:])
                if not(i % 15):
                    #csvwriterows_safe('HisCount15_Old_%d.csv' % ((i+15-1)//15), map(None, *clustersEI),
                    #                  '%scountshistoNeurClasses\/' % netbrs.data_path) 
                    #csvwriterows_safe('SpikesCount15_Old_%d.csv' % ((i+15-1)//15), map(None, *clustersEI_xbins),
                    #                  '%scountshistoNeurClasses\/' % netbrs.data_path) 
                    csvwriterows_safe('NetStates15_%d.csv' % ((i+15-1)//15), map(None, *NetStates),
                                      '%snb_clusters\/' % netbrs.data_path) 
                    csvwriterows_safe('NetStatesNN15_%d.csv' % ((i+15-1)//15), map(None, *NetStatesNN),
                                      '%snb_clusters\/' % netbrs.data_path) 
                    csvwriterows_safe('NetStatesEX15_%d.csv' % ((i+15-1)//15), map(None, *NetStatesEX),
                                      '%snb_clusters\/' % netbrs.data_path) 
                    if store_NeuronsCounts:
                        csvwriterows_safe('HisCount15_%d.csv' % ((i+15-1)//15), map(None, *HisCounts),
                                          '%scountshistoNeurClasses\/' % netbrs.data_path) 
                        csvwriterows_safe('CumCount15_%d.csv' % ((i+15-1)//15), map(None, *CumCounts),
                                          '%scountshistoNeurClasses\/' % netbrs.data_path) 
                        csvwriterows_safe('SpikesCount15_%d.csv' % ((i+15-1)//15), map(None, *SpikeCounts),
                                          '%scountshistoNeurClasses\/' % netbrs.data_path) 
                    if store_NetStatesNNcum:
                        csvwriterows_safe('NetStatesNN15cum_%d.csv' % ((i+15-1)//15), map(None, *NetStatesNNc),
                                          '%snb_clusters\/' % netbrs.data_path) 
                    if store_FF:
                        csvwriterows_safe('NetStatesFF15_%d.csv' % ((i+15-1)//15), map(None, *FF),
                                          '%snb_clusters\/' % netbrs.data_path) 
                    if Net_features_bin['states_nr'] > 2:
                        csvwriterows_safe('NetStates15_34_%d.csv' % ((i+15-1)//15), map(None, *NetStates34),
                                          '%snb_clusters\/' % netbrs.data_path) 
                        csvwriterows_safe('NetStatesNN15_34_%d.csv' % ((i+15-1)//15), map(None, *NetStatesNN34),
                                          '%snb_clusters\/' % netbrs.data_path) 
                        csvwriterows_safe('NetStatesEX15_34_%d.csv' % ((i+15-1)//15), map(None, *NetStatesEX34),
                                          '%snb_clusters\/' % netbrs.data_path) 
                        if store_NetStatesNNcum:
                            csvwriterows_safe('NetStatesNN15cum_34_%d.csv' % ((i+15-1)//15), map(None, *NetStatesNNc34),
                                              '%snb_clusters\/' % netbrs.data_path) 
                    NetStates, NetStatesNN, NetStatesNNc, NetStatesEX, clustersEI, clustersEI_xbins = [], [], [], [], [], []
                    NetStates34, NetStatesNN34, NetStatesNNc34, NetStatesEX34 = [], [], [], []
                    HisCounts, CumCounts, SpikeCounts, FF = [], [], [], []
            if i % 15:
                #csvwriterows_safe('HisCount15_Old_%d.csv' % ((i+15-1)//15), map(None, *clustersEI),
                #                  '%scountshistoNeurClasses\/' % netbrs.data_path) 
                #csvwriterows_safe('SpikesCount15_Old_%d.csv' % ((i+15-1)//15), map(None, *clustersEI_xbins),
                #                  '%scountshistoNeurClasses\/' % netbrs.data_path) 
                csvwriterows_safe('NetStates15_%d.csv' % ((i+15-1)//15), map(None, *NetStates),
                                  '%snb_clusters\/' % netbrs.data_path) 
                csvwriterows_safe('NetStatesNN15_%d.csv' % ((i+15-1)//15), map(None, *NetStatesNN),
                                  '%snb_clusters\/' % netbrs.data_path) 
                csvwriterows_safe('NetStatesEX15_%d.csv' % ((i+15-1)//15), map(None, *NetStatesEX),
                                  '%snb_clusters\/' % netbrs.data_path) 
                if store_NeuronsCounts:
                    csvwriterows_safe('HisCount15_%d.csv' % ((i+15-1)//15), map(None, *HisCounts),
                                      '%scountshistoNeurClasses\/' % netbrs.data_path) 
                    csvwriterows_safe('CumCount15_%d.csv' % ((i+15-1)//15), map(None, *CumCounts),
                                      '%scountshistoNeurClasses\/' % netbrs.data_path) 
                    csvwriterows_safe('SpikesCount15_%d.csv' % ((i+15-1)//15), map(None, *SpikeCounts),
                                      '%scountshistoNeurClasses\/' % netbrs.data_path) 
                if store_NetStatesNNcum:
                    csvwriterows_safe('NetStatesNN15cum_%d.csv' % ((i+15-1)//15), map(None, *NetStatesNNc),
                                      '%snb_clusters\/' % netbrs.data_path) 
                if store_FF:
                    csvwriterows_safe('NetStatesFF15_%d.csv' % ((i+15-1)//15), map(None, *FF),
                                      '%snb_clusters\/' % netbrs.data_path) 
                if Net_features_bin['states_nr'] > 2:
                    csvwriterows_safe('NetStates15_34_%d.csv' % ((i+15-1)//15), map(None, *NetStates34),
                                        '%snb_clusters\/' % netbrs.data_path) 
                    csvwriterows_safe('NetStatesNN15_34_%d.csv' % ((i+15-1)//15), map(None, *NetStatesNN34),
                                        '%snb_clusters\/' % netbrs.data_path) 
                    csvwriterows_safe('NetStatesEX15_34_%d.csv' % ((i+15-1)//15), map(None, *NetStatesEX34),
                                        '%snb_clusters\/' % netbrs.data_path) 
                    if store_NetStatesNNcum:
                        csvwriterows_safe('NetStatesNN15cum_34_%d.csv' % ((i+15-1)//15), map(None, *NetStatesNNc34),
                                            '%snb_clusters\/' % netbrs.data_path) 
            csvwriterows_safe('netstate_data.csv', sortingInfos, netbrs.data_path) 
            netdata2 = numpy.array(netdata2)
            netdata2ref = numpy.array(netdata2ref)
            netdata2refmean = netdata2ref.mean(axis=0)
            netdata2refmean[1+numpy.arange(0,len(netdata2refmean),2)] = netdata2refmean[numpy.arange(0,len(netdata2refmean),2)]
            netdata2 = numpy.hstack((netdata2, netdata2 / netdata2refmean))
            netdata2final = [['step','IBI','sem','SNE','sem','SNI','sem','NNE','sem','NNI','sem',
                              'EXE','sem','EXI','sem','REXEI','sem','SRE','sem','SRI','sem',
                              'IBIn','sem','SNEn','sem','SNIn','sem','NNEn','sem','NNIn','sem',
                              'EXEn','sem','EXIn','sem','REXEIn','sem','SREn','sem','SRIn','sem']]
            for i, _step in enumerate(self.exp_steps): netdata2final.append([_step[0]] + list(netdata2[i])) 
            csvwriterows_safe('netstate_data2.csv', netdata2final, netbrs.data_path) 
            csvwriterows_safe('neurSilent.csv', neursSilent, netbrs.data_path+'countshistoNeurClasses\/') 

            reply = raw_input("To exit digit N, to change settings Enter: ")
            if reply == 'N': break
            try:
                reply = raw_input("Input new trash_poor_states or Enter: ")
                if reply != '':
                    npn.trash_poor_states = self.trash_poor_states = int(reply)
                    npn.data_path = self.Net_path = '%s_trash%s\/' % (npn.data_path[:-2],
                                                                      self.trash_poor_states)
                reply = raw_input("Input new Net_features(comma separated) or Enter: ")
                if reply != '':
                    npn.Net_features = self.Net_features = reply.split(',')
                    npn.data_path = self.Net_path = '%s_feat%s\/' % (npn.data_path[:-2], reply)
                reply = raw_input("Input new bin or Enter: ")
                if reply != '':
                    npn.Net_features_bin['bin'] = self.Net_features_bin['bin'] = float(reply)
                    npn.data_path = self.Net_path = '%s_bin%s\/' % (npn.data_path[:-2],
                                                                    self.Net_features_bin['bin'])
                reply = raw_input("Input new binstates or Enter: ")
                if reply != '':
                    npn.Net_features_bin['binstates'] = self.Net_features_bin['binstates'] = float(reply)
                    npn.data_path = self.Net_path = '%s_binstates%s\/' % (npn.data_path[:-2],
                                                                          self.Net_features_bin['binstates'])
                reply = raw_input("Input new delay or Enter: ")
                if reply != '':
                    npn.Net_features_bin['delay'] = self.Net_features_bin['delay'] = float(reply)
                    npn.data_path = self.Net_path = '%s_delay%s\/' % (npn.data_path[:-2],
                                                                      self.Net_features_bin['delay'])
                reply = raw_input("Input new pvalue or Enter: ")
                if reply != '':
                    npn.Net_features_bin['pvalue'] = self.Net_features_bin['pvalue'] = float(reply)
                    npn.data_path = self.Net_path = '%s_pvalue%s\/' % (npn.data_path[:-2],
                                                                       self.Net_features_bin['pvalue'])
                reply = raw_input("Input new nr_pvalue or Enter: ")
                if reply != '':
                    npn.Net_features_bin['nr_pvalue'] = self.Net_features_bin['nr_pvalue'] = int(reply)
                    npn.data_path = self.Net_path = '%s_nr_pvalue%s\/' % (npn.data_path[:-2],
                                                                          self.Net_features_bin['nr_pvalue'])
            except ValueError: break
        
        self.proc_flow_times.append(('Network burst detection', time.clock() - t_start))       

    def write_results(self):

        ############################################################################
        # WRITE ALL RESULTS
        ############################################################################
        print 'Writing results data..'
        t_start = time.clock()
        
        self._blk_col_nr -= self.clusters_nr * (self.clusters_nr <= self._blk_col_nr)
        csv_path = '%sclusters_data' % self.results_path
        if self.clustered_data.Outliers_detect:
            csv_path += '__without_Outliers'
            outliers_names = self.clustered_data.all_outliers_names
        else: outliers_names = {}
        self.write_clustering_burst_results(csv_path, all_outliers_names = outliers_names)
        self.proc_flow_times.append(('Results storage', time.clock() - t_start))
    
    def write_clustering_burst_results(self, csv_path, all_outliers_names = []):
        """Writing all results of clustering processing in the *clusters_data*.csv file,
        storing results about features values in each cluster, for each experimental stage.
        At first we calculate the number of blank columns to insert in the results data file. 
        If  'Release' in globals(), we replace part of the first row of *clusters_data*.csv
        with the nr. of neurons and outliers in each cluster, experiment step name(s),
        burst features as needed in order to plot them in Origin with all
        available informations.
        """

        ############################################################################
        # All data to write will be appended to list _to_clusters_data in the next
        ############################################################################
        _to_clusters_data = []
        _to_clusters_norm_data = []

        first_row = ['experimental stage'] + \
                    deep_chain_new([[[feat+' cluster '+repr(i), feat+' sem'] \
                    for i in self.clusters_names_sorted_x_nr[3]] + ['']*2*self._blk_col_nr \
                    for feat in desc_burst_features])
                
        ############################################################################
        # Over-writing the first row with the informations used in Origin
        ############################################################################
        info_to_Origin = deep_chain_new([(i, len(all_outliers_names[j])) \
                                        if j in all_outliers_names else (i, None) \
                                        for i, j in zip(self.clusters_names_sorted_x_nr[1], \
                                                        self.clusters_names_sorted_x_nr[0])])
        info_to_Origin.extend(['']*(12-len(info_to_Origin)))

        first_row[1:len(info_to_Origin) + 1] = info_to_Origin
        _to_clusters_data.append(first_row)
            
        ############################################################################
        # Collecting all clustering results, then saving them in clusters_data.csv
        # Add blank columns if clusters nr. < self._blk_col_nr 
        ############################################################################
        results = dict((step, map(centroid_plus_sem, \
                                ([bbox_step_data.desc_bursts_stat[neuron] for neuron in cl] \
                                for cl in self.clusters_names_sorted_x_nr[2]))) \
                        for step, bbox_step_data in self.bbox.iteritems())
        
        to_norm = []
        to_norm_ref = []
        for k in self.exp_steps:
            step = k[0]
            data_2_csv = [step]
            for j in xrange(len(desc_burst_features)):
                data_2_csv.extend([fpformat.fix(z, 4) if z else '' \
                                   for d in results[step] for z in d[j]] + \
                                   ['']*2*self._blk_col_nr)
            _to_clusters_data.append(data_2_csv)
            # Start parte nuova 28/06/2011
            to_norm.append(deep_chain_new(results[step]))
            if step in self.ref_steps: to_norm_ref.append(deep_chain_new(results[step]))
        #print to_norm_ref
        
        to_norm = numpy.array(to_norm)
        to_norm_ref = numpy.array(to_norm_ref)
        to_norm_ref = numpy.mean(to_norm_ref, axis=0)
        res_2_norm = []
        for i, k in enumerate(self.exp_steps):
            step = k[0]
            if step in self.ref_steps: to_norm[i][::2] = 1
            else: to_norm[i][::2] -= to_norm_ref[::2]-1
            res_temp = [step]
            _j = 2*len(desc_burst_features)
            for j in xrange(len(desc_burst_features)):
                for z in xrange(min(6,self.clusters_nr)):
                    res_temp.extend(list(to_norm[i][2*j + z*_j:2*j + z*_j + 2]))
                res_temp.extend(['']*2*self._blk_col_nr)
            res_2_norm.append(res_temp)
        # End parte nuova 28/06/2011

        ############################################################################
        # Writing neurons Outliers names in each cluster, if Outliers
        # detection was applied, inside file *_clusters_data__no_Outliers.csv
        ############################################################################
        if self.clustered_data.Outliers_detect:
            _to_clusters_data.extend([[' ']] + [['Outliers in cluster ' + str(i) + ':'] + all_outliers_names[j] \
                                    for i, j in zip(self.clusters_names_sorted_x_nr[3], \
                                                    self.clusters_names_sorted_x_nr[0]) \
                                    if j in all_outliers_names])

        _to_clusters_data.extend([[], []])
        for n, c in zip(self.clusters_names_sorted_x_nr[3], self.clusters_names_sorted_x_nr[2]):
            _to_clusters_data.append([])
            data_2_csv = ['cluster ' + repr(n)] + sorted(c, cmp = sort_neuro_names)
            _to_clusters_data.append(data_2_csv)
            for k in self.exp_steps:
                data_2_search = self.bbox[k[0]].desc_bursts_stat
                data_2_csv = [k[0]]
                for i in sorted(c, cmp = sort_neuro_names):
                    if data_2_search[i] != [None] * len(desc_burst_features):
                        data_2_csv.append(','.join(map(str, data_2_search[i])))
                    else: data_2_csv.append(None)
                _to_clusters_data.append(data_2_csv)
        
        _to_clusters_data[0][5] = ','.join([str(self.last_feat_kw['ACF_delay']),
                                            str(self.last_feat_kw['Fano_factor_bin']),
                                            str(self.last_feat_kw['Mahalanobis_threshold'])])
        _to_clusters_data[0][6] = str(Net_features_bin['pvalue'])
        _to_clusters_data[0][7] = str(Net_features_bin['Mahalanobis_threshold'])
        _to_clusters_data[0][8] = str(min_interburst_length)
        _to_clusters_data[0][9] = ','.join(Net_features)
        _to_clusters_data[0][10] = str(Net_features_bin['binstates'])
        _to_clusters_data[0][11] = str(Netlow_threshold)
        _to_clusters_data[0][12] = str(Netburst_binlength)
        _to_clusters_data[0][13] = str(Nethigh_threshold) 
        _to_clusters_data[0][14] = str(self.exp_name) 
        _to_clusters_data[0][15] = ','.join(burst_features) 
        _to_clusters_data[0][16] = str(Net_features_bin['nr_pvalue'])
        _to_clusters_data[0][17:17+len(self.ref_steps)] = self.ref_steps
        _to_clusters_data[0][20] = str(trash_poor_states)
        _to_clusters_data[0][23] = str(Net_features_bin['bin'])
        _to_clusters_data[0][24] = str(Net_features_bin['binFF'])
        
        csvwriterows_safe(self.csv_path + '.csv', _to_clusters_data)    
        
        _to_clusters_norm_data = _to_clusters_data
        _to_clusters_norm_data[1:len(self.exp_steps)+1] = res_2_norm
        csvwriterows_safe(self.csv_path + '_norm.csv', _to_clusters_norm_data)    

        ############################################################################
        # Storing histograms of ACF, ISI, IBI and other features calculated in
        # burstdetection.py, grouping them by known clusters.
        ############################################################################
        for step in self.exp_steps:
            self.store_data_feat_to_pca(step[0], self.clustered_data.clusters_names.values())
            if store_allburstsdata:
                self.store_desc_feat_histogram(step[0], self.clustered_data.clusters_names.values())
    
    def write_burst_results(self):
        """This function in invoked to store all bursts data for each neuron during a
        given time interval on *_bursts_data.csv and optionally on a sqlite database.
        """
        
        print 'Writing single neuron bursts data..'
        burst_data_2csv = {}
        for _step in self.exp_steps:
            for k, k_data in self.bbox[_step[0]].bursts_all_data.iteritems():
                burst_data_2csv.setdefault(k, []).extend([[k] + b[2:] for b in k_data]) 
        
        bsf = burst_data_storage(self.bursts_data_neurons,
                                 burst_data_2csv,
                                 self.results_path)
        bsf.csv_format()
        
        ############################################################################
        # Storing bursts data for each neuron in each cluster
        ############################################################################
        burst_data_2csv = {}
        for _step in self.exp_steps:
            for k, k_data in self.bbox[_step[0]].bursts_all_data.iteritems():
                burst_data_2csv.setdefault(k, []).extend([[k] + b[2:] for b in k_data]) 
        
        try:
            for i, j in zip(self.clusters_names_sorted_x_nr[3], self.clusters_names_sorted_x_nr[2]):
                bursts_extracted_writer = csv.writer(open(bursts_extracted_path + \
                                        'bursts_data_cluster_' + repr(i) + '.csv', 'wb'), **csv_options)
                for k in sorted(j, cmp = sort_neuro_names):
                    bursts_extracted_writer.writerows([burst_data_2csv[k][0]] + \
                                            sorted(burst_data_2csv[k][1:], key = operator.itemgetter(1)))
        except IOError, (errno, strerror):
            print "I/O error(%s): %s" % (errno, strerror)
            raise

    def store_feat_data(self,
                        feat_mean_sem_x_cluster,
                        feat_name, 
                        feat_bin,
                        feat_delay,
                        feat_data, 
                        cluster_nr,
                        cluster_nr_str,
                        clusters_names,
                        step_name = None):
        """Filling the matrix feat_mean_sem_x_cluster with the values 
        of mean+s.e.m. of histogram data grouped by cluster"""
        
        ############################################################################
        # If feat_bin is an iterable. Note that it is expected to be a numpy.array
        # else if feat_bin is a number, we have to create the bins array
        ############################################################################
        if hasattr(feat_bin, '__iter__'): bins_column = [['bin'] + feat_bin.tolist()]
        else: bins_column = [['bin'] + numpy.arange(feat_bin, feat_delay + feat_bin, feat_bin).tolist()]
        
        feat_cluster_data   = numpy.asarray([feat_data[neuron] for neuron in clusters_names[cluster_nr] \
                                        if neuron in feat_data and feat_data[neuron].__len__()])
        len_feat_data       = feat_cluster_data.__len__().__float__()
        
        ############################################################################
        # Before to calculate the cumulative sum 'numpy.cumsum', apply
        # shape_normalize to feat_cluster_data_shape_norm_mean data, otherwise the
        # cumulative sum will not join an unitary value, because of all zeros array,
        # if present, lowering the sum.
        ############################################################################
        if len_feat_data:
            
            sqrt_len_feat_data  = numpy.sqrt(len_feat_data)
            tot_nr_event        = numpy.concatenate(feat_cluster_data).sum(dtype = numpy.float32)
            feat_cluster_data_shape_norm = numpy.asarray(shape_normalize(feat_cluster_data))
            feat_cluster_data_shape_norm_mean = numpy.mean(feat_cluster_data_shape_norm, axis = 0, dtype = numpy.float32)
            
            feat_mean_sem_x_cluster.extend(bins_column + \
                                        [[feat_name + cluster_nr_str] + \
                                        numpy.mean(feat_cluster_data, axis = 0).tolist()] + \
                                        [['sem' + feat_name + cluster_nr_str] + \
                                        numpy.divide(numpy.std(feat_cluster_data, axis = 0), \
                                                               sqrt_len_feat_data).tolist()] + \
                                        [[feat_name + 'n' + cluster_nr_str] + \
                                        feat_cluster_data_shape_norm_mean.tolist()] + \
                                        [['sem' + 'n' + feat_name + cluster_nr_str] + \
                                        numpy.divide(numpy.std(feat_cluster_data_shape_norm, axis = 0), \
                                                               sqrt_len_feat_data).tolist()] + \
                                        [[feat_name + 'c' + cluster_nr_str] + \
                                        numpy.cumsum(shape_normalize(feat_cluster_data_shape_norm_mean)).tolist()] + \
                                        [[''] * (int(len_feat_data) + 1)])
        else:

            feat_mean_sem_x_cluster.extend(bins_column + \
                                        [[feat_name + cluster_nr_str], feat_cluster_data.tolist(), \
                                        ['sem'] + feat_cluster_data.tolist()] * 3)

        ############################################################################
        # Adding the name name of the experimental step
        ############################################################################
        feat_mean_sem_x_cluster[6][1] = step_name

    def store_data_feat_to_pca(self,
                               exp_steps,
                               clusters_names):
        """Storing data of ACF, ISI or IBI histograms if present, grouped by cluster,
        for all steps in exp_steps together"""
            
        try:
            if not(hasattr(exp_steps, '__iter__')): exp_steps = [exp_steps]
            
            ############################################################################
            # Writing histogram data of ACF, if used, in a single file
            ############################################################################
            feat_to_pca_filename = "%sACF_%s.csv" % (self.feat_storing_path, '-'.join(exp_steps))
            feat_to_pca_x_cluster_writer = csv.writer(open(feat_to_pca_filename, 'wb'), **csv_options)
            
            feat_mean_sem_x_cluster = []
            _ACF_delay = self.last_feat_kw['ACF_delay']
            _ACF_bin = _ACF_delay / 100.
            for i, j in zip(self.clusters_names_sorted_x_nr[0], self.clusters_names_sorted_x_nr[3]):
                for step in exp_steps:
                    cluster_nr_str = repr(j)
                    self.store_feat_data(feat_mean_sem_x_cluster,
                                         'ACF', 
                                         _ACF_bin,
                                         _ACF_delay,
                                         self.bbox[step].exp_step_acf, 
                                         i,
                                         cluster_nr_str,
                                         clusters_names,
                                         step)

            feat_to_pca_x_cluster_writer.writerows(map(None, *feat_mean_sem_x_cluster))            
            
            if store_allburstsdata:

                ############################################################################
                # Writing only ISI histogram data
                ############################################################################
                feat_to_pca_filename = self.feat_storing_path + 'ISI_' + '-'.join(exp_steps) + '.csv'
                feat_to_pca_x_cluster_writer = csv.writer(open(feat_to_pca_filename, 'wb'), **csv_options)
                feat_mean_sem_x_cluster = []
                for i, j in zip(self.clusters_names_sorted_x_nr[0], self.clusters_names_sorted_x_nr[3]):
                    for step in exp_steps: self.store_feat_data(feat_mean_sem_x_cluster,
                                                                'ACF', #It's ISI,write ACF in file after request  
                                                                ISI_bin,
                                                                ISI_delay,
                                                                self.bbox[step].exp_step_isi_histo, 
                                                                i,
                                                                repr(j),
                                                                clusters_names)
                feat_to_pca_x_cluster_writer.writerows(map(None, *feat_mean_sem_x_cluster))          
                
                ############################################################################
                # Writing IBI, if applied, histogram data in a single file
                ############################################################################
                if set(['IBI_PC1', 'IBI_PC2', 'IBI_PC3']) <= set(all_features) and 'IBI' in all_features:
                    feat_to_pca_filename = self.feat_storing_path + 'IBI_' + '-'.join(exp_steps) + '.csv'
                    feat_to_pca_x_cluster_writer = csv.writer(open(feat_to_pca_filename, 'wb'), **csv_options)
                
                    feat_mean_sem_x_cluster = []
                    for i, j in zip(self.clusters_names_sorted_x_nr[0], self.clusters_names_sorted_x_nr[3]):
                        for step, start, end in exp_steps:
                            cluster_nr_str = repr(j)
                            self.store_feat_data(feat_mean_sem_x_cluster,
                                                 'IBI',
                                                 IBI_bin, 
                                                 IBI_delay, 
                                                 self.bbox[step].exp_step_ibis,
                                                 i,
                                                 cluster_nr_str,
                                                 clusters_names)
                    
                    feat_to_pca_x_cluster_writer.writerows(map(None, *feat_mean_sem_x_cluster))

        except IOError, (errno, strerror):
            print "I/O error(%s): %s" % (errno, strerror)
            raise
        except:
            print "Unexpected error:", sys.exc_info()[0]
            raise


    def store_desc_feat_histogram(self,
                                  exp_steps,
                                  clusters_names):
        """Storing features_to_histo histogram data grouped by cluster
        and for all steps in exp_steps together"""
            
        try: 
            if not(hasattr(exp_steps, '__iter__')): exp_steps = [exp_steps]
            feat_to_histo_filename = ''.join([self.feat_storing_path, 'HIS_', '-'.join(exp_steps), '.csv'])
            feat_to_histo_x_cluster_writer = csv.writer(open(feat_to_histo_filename, 'wb'), **csv_options)
            
            feat_mean_sem_x_cluster = []
            for i, j in zip(self.clusters_names_sorted_x_nr[0], self.clusters_names_sorted_x_nr[3]):
                for step in exp_steps:
                    cluster_nr_str = repr(j)
                    for feat_to_histo in features_to_histo & set(all_features):
                        self.store_feat_data(feat_mean_sem_x_cluster, 
                                            feat_to_histo,
                                            numpy.arange(self.bbox[step].start_end_bin[feat_to_histo][0] + \
                                                         self.bbox[step].start_end_bin[feat_to_histo][2],
                                                         self.bbox[step].start_end_bin[feat_to_histo][1],
                                                         self.bbox[step].start_end_bin[feat_to_histo][2]),
                                            None,
                                            self.bbox[step].all_feat_to_histo_data[feat_to_histo], 
                                            i,
                                            cluster_nr_str,
                                            clusters_names,
                                            step)
            feat_to_histo_x_cluster_writer.writerows(map(None, *feat_mean_sem_x_cluster))
        except IOError, (errno, strerror):
            print "I/O error(%s): %s" % (errno, strerror)
            raise
        except:
            print "Unexpected error:", sys.exc_info()[0]
            raise


class timeneuro_mod_run(object):
    
    def __get_exp_to_process_file(self): return 'experiment_to_process.txt'
    exp_to_process_file = property(__get_exp_to_process_file)
    
    @repoRelease
    @catchException
    @addMEAdata
    def __init__(self):

        #################################################################################
        #####                                                                       #####
        #####       Starting the (Max Interval) single neuron burst detection       #####
        #####       and Network Burst Detection (if enabled in the setting file)    #####
        #####                                                                       #####
        #################################################################################
        for path, name, selection in return_exp_to_process(self.exp_to_process_file):
            
            self.moreinfo = '%s%s' % ('Experiment name: ', name)
            self.stderr_file = path + '\/MEAdata_cluster_results\/logError.txt'
            
            timeneuro_processing.selection          = selection
            timeneuro_processing._blk_col_nr        = globals().get('blank_col_nr', False)
            timeneuro_processing.store_results      = ('Release' in globals()) | globals().get('store_results', False)
            timeneuro_processing.Netburst_detection = ('Release' in globals()) | globals().get('Netburst_detection', False)  
            timeneuro_processing.randomization      = globals().get('burst_randomization', False) 
            timeneuro_processing.verbose            = globals().get('verbose_screen', False) 
            timeneuro_processing.store_single_neuron_bursts = globals().get('store_single_neuron_bursts', False)  
            timeneuro_processing(path, name)()
            print "%s%s" % ("-" * 70, "\n")


################################################################################################
#####                                                                                      #####
#####                              '__main__' of timeneuro_mod                             #####      
#####                                                                                      #####
################################################################################################
if __name__ == "timeneuro_mod":
    
    timeneuro_mod_run()
