try: 
    import psyco
    psyco.full()
except ImportError: pass

import math
import operator
import itertools
import collections                    



__author__ = "Andrea Maffezzoli"
__version__ = "$Revision: 2.31 - 01 Sep 2010$"
__name__ = "statfuncs"
__all__ = ["factorial", "array_pruning", "merge_dict", 
           "meanval", "standdev", "sem", "variance", "skewness", "kurtosis",
           "stat_moment_mean", "second_moment", "third_moment", "fourth_moment",
           "coeff_variation", "coeff_variation2", "Fano_factor",
           "centroid", "centroid_plus_sem",
           "Mahalanobis_dist", "Mahalanobis_outliers",
           "poisson_pd", "isi_return", "isi_return_gen",
           "convolve_with_pdf", "normalize", "shape_normalize",
           "combinations", "dotproduct", "dotproduct_normalized",
           "acf", "acf_binned", "ac_isi","acf_OLD_VER",
           "ccf", "ccf_binned", "mi_binned",
           "tuplize", "tuplize_iter",
           "deep_chain", "deep_chain_new", "deep_chain_iter",
           "factory_cmp", "makeIterable", "makeNonStringIterable",
           "isIterable", "isNonStringIterable","concatsubseq",
           "myXrange", "purgeISI"]


#########################################################
# We have to append here to __all__ the modules imported
#########################################################
try:
    import numpy
    __all__.append('numpy')
except ImportError:
    raise ImportError, "Impossible to import numpy libraries"
try:
    import scipy
    __all__.append('scipy')
    import scipy.stats
except ImportError:
    raise ImportError, "Impossible to import scipy/scipy.stats libraries"
try:
    import pylab
    __all__.append('pylab')
except ImportError:
    pass #raise ImportError, "Impossible to import pylab/matplotlib module"


##############################################################################
# From http://stackoverflow.com/questions/1464028/what-is-the-pythonic-way-of-checking-if-an-object-is-a-list
# See also A.Martelli's recipes in Python Cookbook
##############################################################################
def isIterable(x): return isNonStringIterable(x, False)

def isNonStringIterable(x, string_as_scalar = True):
    """
    Read:
    - http://stackoverflow.com/questions/1549801/differences-between-isinstance-and-type-in-python
    - http://stackoverflow.com/questions/2158395/flatten-an-irregular-list-of-lists-in-python
    """
    
    try:
        iter(x)
        if isinstance(x, basestring):
            if string_as_scalar: return False
            if len(x) == 1: return False
            return True
    except: return False
    else: return True

def makeIterable(x): return makeNonStringIterable(x, False)

def makeNonStringIterable(x, string_as_scalar = True):
    
    if string_as_scalar: 
        if isinstance(x, basestring): return (x,)
    try: return iter(x)
    except: return (x,)



def factory_cmp(a, b):
            
    def _inner_cmp(x): return operator.le(x, b) and operator.ge(x, a)
    return _inner_cmp



def tuplize(L):
    """This function transforms recursively all iterables in input
       to a tuple, e.g. because list is an unhashable object and
       so it can't become a key of a dictionary
    """

    if not(hasattr(L, '__iter__')): return L
    return tuple(tuplize(t) for t in L) #tuple(map(tuplize, L))



def tuplize_iter(L):
    """This function transforms recursively all iterables in input
    to a tuple, e.g. because list is an unhashable object and so
    it can't become a key of a dictionary
    """

    if not(hasattr(L, '__iter__')): yield L
    for l in tuple(tuplize(t) for t in L): yield l #tuple(map(tuplize, L))



def deep_chain(L, chain_string = False):
    """This function concatenates recursively all iterables in input.
    Input itself is concatenated, so there is no output.
    The input chain_string gives the choice to consider iterable the
    strings too, through the check iter("a string"), which is True
    for strings, whereas hasattr("a string", '__iter__') is False.
    At regard see function isNonStringIterable above.
    The check isinstance(l, collections.Iterator) allows to chain
    iterators objects too, after have saved them as static objects
    (a tuple in the present case).
    Read:
    - http://stackoverflow.com/questions/2158395/flatten-an-irregular-list-of-lists-in-python
    """
    
    itercheck = lambda x: isNonStringIterable(x, not(chain_string))
    for i, l in enumerate(L):
        if itercheck(l):
            if isinstance(l, collections.Iterator): l = tuple(l)
            deep_chain(l, chain_string); L[i:i+1] = l


def deep_chain_new(L, new = False, chain_string = False):
    """This function concatenates recursively all iterables and not
    in input, returning the new concatenated list"""

    if new: deep_chain(L, chain_string); return L

    if not(hasattr(L, '__iter__')): return L
    R = []
    for l in L:
        if hasattr(l, '__iter__'): R.extend(deep_chain_new(l))
        else: R.append(l)
    return R


def deep_chain_iter(L, new = False, chain_string = False):
    """This function concatenates recursively all iterables and not
    given in input and it works as a generator, giving as result
    an iterator. Inspired by ActiveState Code's recipe at:
    http://code.activestate.com/recipes/576719/
    """

    if new:
        deep_chain(L, chain_string)
        for l in L: yield l

    if not(hasattr(L, '__iter__')): yield L
    for l in L:
        if hasattr(l, '__iter__'):
            for y in deep_chain_iter(l): yield y
        else: yield l



def pairwise(iterable):
    """From example code in itertools module of Python 2.6
    (http://docs.python.org/library/itertools.html#module-itertools):
    s -> (s0,s1), (s1,s2), (s2, s3), ...
    It could be slightly modified by using zip instead of izip
    in order to return a list of iterables samples, and I added
    the case of an empty iterable in input.
    """

    if not(iterable): return []
    a, b = itertools.tee(iterable)
    b.next()
    return itertools.izip(a, b)



def merge_dict(dict1, dict2):
    """A function able to merge the two dictionaries given in input"""
    
    if not(isinstance(dict1, dict)) or not(isinstance(dict2, dict)): return None
    for x, y in dict2.items():
        if x in dict1:
            if not isinstance(dict1[x], list): dict1[x] = [dict1[x]]
            if not isinstance(dict2[x], list): dict2[x] = [dict2[x]]
            dict1[x] += dict2[x]
        else: dict1[x] = dict2[x]
    return dict1



def factorial(n, factorials_dict = {}):
    """Factorial of the number in input.
    This last version of the function uses memoization
    in order to calculate factorials as fast as possible
    """

    if factorials_dict: i_start = max(factorials_dict.iterkeys()) + 1
    else:
        factorials_dict.update({0: 1})
        i_start = 1
    for i in xrange(i_start, n + 1): factorials_dict[i] = i * factorials_dict.get(i - 1)
    return factorials_dict[n]



def array_pruning(L):
    """Removing None values from an array.
    NB: Use tests for object identity when appropriate: if x is not None rather than if x != None. 
          It is much more efficient to test objects for identity than equality, because identity only
          checks their address in memory. Two objects are identical if they are the same object in the same
          physical location and not their actual data. From http://jaynes.colorado.edu/PythonIdioms.html)
    """

    L_pruned = [i for i in L if i is not None] # filter(lambda x: x is not None, L)
    if not(L_pruned): return None
    return L_pruned



def meanval(L):
    """Mean values of values in a list or list of lists"""

    if not(hasattr(L, '__iter__')): return None
    L_pruned = array_pruning(L)
    if not(L_pruned): return None
    if type(L_pruned) == type(numpy.array([])): return numpy.mean(L_pruned, dtype = numpy.float32)
    return sum(L_pruned) / float(len(L_pruned))



def stat_moment_mean(L, order):
    """The statistical moment about the mean"""

    if not(isinstance(order, int)): raise Exception, "The order of the statistical moment has to be an integer"
    ####################################################################################################
    # Use tests for object identity when appropriate: if x is not None rather than if x != None. 
    # It is much more efficient to test objects for identity than equality, because identity only
    # checks their address in memory. Two objects are identical if they are the same object in the same
    # physical location and not their actual data. From http://jaynes.colorado.edu/PythonIdioms.html)
    ####################################################################################################
    if not(hasattr(L, '__iter__')): return None
    L_pruned = array_pruning(L)
    if not(L_pruned): return None
    m = meanval(L_pruned)
    return sum([pow((i - m), order) for i in L_pruned]) / float(len(L_pruned))



def second_moment(L):
    """Second order statitical moment of values of a random variable"""
    
    if type(L) == type(numpy.array([])): return numpy.var(L, dtype = numpy.float32)
    return stat_moment_mean(L, 2)



def third_moment(L):
    """Third order statitical moment of values of a random variable"""

    return stat_moment_mean(L, 3)



def fourth_moment(L):
    """Fourth order statitical moment of values of a random variable"""

    return stat_moment_mean(L, 4)



def standdev(L):
    """Stanadard deviation of values of a random variable"""
    
    if type(L) == type(numpy.array([])): return numpy.std(L, dtype = numpy.float32)
    sec_mom = second_moment(L)
    if sec_mom is None: return None
    return math.sqrt(sec_mom)



def sem(L):
    """Stanadard error of the mean of values of a random variable"""
    
    stdv = standdev(L)
    if stdv is None: return None
    return stdv / math.sqrt(len(L))



def variance(L):
    """Variance of the mean of values of a random variable"""
    
    return second_moment(L)



def coeff_variation(isis, upper_threshold = None, lower_threshold = None):
    """Coefficient of variation (CV).
    It is equal to 1 for a homogeneous Poisson process.
    To be applied on ISI data."""
    
    if upper_threshold: isis = [i for i in isis if i <= upper_threshold]
    if lower_threshold: isis = [i for i in isis if i >= lower_threshold]
    m = meanval(isis)
    if not(m): return None
    stdv = standdev(isis)
    if stdv is None: return None
    return stdv / m



def coeff_variation2(isis, upper_threshold = None, lower_threshold = None):
    """Statistical index applied to estimate spike regularity.
    At regard see the following reference:
        G. R. Holt, W. R. Softky, C. Koch, and R. J. Douglas
        Comparison of discharge variability in vitro and in vivo in cat visual cortex neurons
        J Neurophysiol 75: 1806-1814, 1996.
    To be applied on ISI data."""
    
    if upper_threshold: isis = [i for i in isis if i <= upper_threshold]
    if lower_threshold: isis = [i for i in isis if i >= lower_threshold]
    pairwise_isis = pairwise(isis)
    return meanval((2 * math.fabs(i - j) / (i + j) for i, j in pairwise_isis))



def Fano_factor_OLD(C, upper_threshold = None, lower_threshold = None):
    """Fano factor.
    It is equal to 1 for a homogeneous Poisson process, and, for any
    renewal process, the Fano factor evaluated for long time intervals 
    approaches the value of CV**2.
    To be applied on spike counts per bin data."""

    if upper_threshold: C = [i for i in C if i <= upper_threshold]
    if lower_threshold: C = [i for i in C if i >= lower_threshold]
    m = meanval(C)
    if not(m): return None
    var = variance(C)
    if var is None: return None
    return var / m


def Fano_factor(C, upper_threshold = None, lower_threshold = None):
    """Fano factor.
    It is equal to 1 for a homogeneous Poisson process, and, for any
    renewal process, the Fano factor evaluated for long time intervals 
    approaches the value of CV**2.
    To be applied on spike counts per bin data."""

    if upper_threshold is None and lower_threshold is None:
        m = numpy.mean(C)
        if m == 0: return None
        v = numpy.var(C)
    else:
        if lower_threshold is None: condition = C <= upper_threshold
        elif upper_threshold is None: condition = C >= lower_threshold
        else: condition = numpy.logical_and(C >= lower_threshold, C>= upper_threshold)
        m = numpy.mean(numpy.extract(condition, C))
        if m == 0: return None
        v = numpy.var(numpy.extract(condition, C))
    return v / m


def skewness(L):
    """Skewness of values of a random variable
    To apply preferably on ISI data."""
    
    if isinstance(L, int) or isinstance(L, float): return 0.0
    std = standdev(L)
    if not(std): return None
    third_mom = third_moment(L)
    if third_mom is None: return None
    return third_mom / pow(std, 3)



def kurtosis(L):
    """Kurtosis of values of a random variable"""
    
    if isinstance(L, int) or isinstance(L, float): return 0.0
    second_mom = second_moment(L)
    if not(second_mom): return None
    fourth_mom = fourth_moment(L)
    if fourth_mom is None: return None
    return fourth_mom / pow(second_mom, 2) - 3



def centroid(M):
    """This function return the centroid of data in M, having shape: 
    [[var1, var2, ...varm] for sample 1, ...[var1, ...varm] for sample n"""
    
    l = len(M)
    if l > 1:
        if isinstance(M[0], int) | isinstance(M[0], float):
            ######################################
            # When working with only one variable
            ######################################
            return meanval(M)
        M_transpose = map(None, *M)
        return [meanval(i) for i in M_transpose]
    ########################################################
    # In the case the matrix M is composed by only one list
    ########################################################
    elif l == 1: return list(M[0])
    ###########################################
    # In the case the matrix M has zero lenght
    ###########################################
    else: return [None]



def centroid_plus_sem(M):
    """This function return the centroid plus standard deviation 
    of error (SEM) of data in M, having shape: 
    [[var1, var2, ...varm] for sample 1, ...[var1, ...varm] for sample n]
    """
    
    if not(hasattr(M, '__iter__')): return [None, None]
    l = len(M) 
    if l > 1:
        ##############################################
        # If M is an iterable of not iterable objects
        ##############################################
        if not(hasattr(M[0], '__iter__')): return [meanval(M), sem(M)]
        ##########################################
        # If M is an iterable of iterable objects
        ##########################################
        else:
            M_transpose = map(None, *M)
            return [[meanval(i), sem(i)] for i in M_transpose]
    ########################################################
    # In the case the matrix M is composed by only one list
    ########################################################
    elif l == 1: return [[i, 0] for i in M[0]]
    ###########################################
    # In the case the matrix M has zero lenght
    ###########################################
    else: return [None, None]



def Mahalanobis_dist(M, rowvar_input = 0):
    """It returns a dict of the Mahahalanobis distances between all samples of M,
    which will be the keys of the dict, and the centroid of M"""
    
    try:
        cluster_inv_cov_matrix = numpy.linalg.inv(numpy.cov(M, rowvar = rowvar_input))
        cluster_centroid = numpy.array(centroids(M))
        
        if not(isinstance(M, dict)): M = dict(enumerate(M))

        return dict([(i, numpy.dot( \
                                numpy.dot(numpy.array(M[i]) - cluster_centroid, \
                                          cluster_inv_cov_matrix), \
                                numpy.transpose(numpy.array(M[i]) - cluster_centroid))) \
                        for i in M])
    except LinAlgError:
        raise LinAlgError, 'Matrix in input is singular: impossible to invert its covariance matrix.'



def Mahalanobis_outliers(M, std_threshold = 2.0):
    
    if not(isinstance(M, dict)): M = dict(enumerate(M))
    Mahalanobis_d = Mahalanobis_dist(M)
    mean_Mahalanobis_d = meanval(Mahalanobis_d.values())
    std_Mahalanobis_d = standdev(Mahalanobis_d.values())
    return [(i, j) for i, j in M.items() if math.fabs(j - mean_Mahalanobis_d) > std_Mahalanobis_d * std_threshold]



def poisson_pd(Freq, time_interval, n):
    """Poisson Homogeneous probability function"""
    
    d = Freq * time_interval
    return pow(d, n) * math.exp(-d) / factorial(n)



def isi_return(T):
    """This function return a list of isis from the timestamps list in input"""

    return map(operator.sub, T[1:], T[:-1])


def isi_return_gen(T):
    """This function return a list of isis from the timestamps list in input"""
    
    T1, T0 = itertools.tee(T)
    T1.next()
    for t in T1: yield t - T0.next()


def acf_OLD_VER(T,
                bin,
                delta, 
                start = 0,
                two_side = True, 
                no_bias = True,
                no_AC_zero = True,
                time_normalized = True,
                bin_normalized = False):
    """
    I HOLD THIS ONLY FOR RETRO-COMPATIBILITI ISSUES BEACUSE USED IN TIMENEUROMOD
    Assume data T already sorted.
    Calculating autocorrelation for timestamps data T, as explained in:
    - P Dayan, LF Abbott. Theoretical Neuroscience (MIT Press). 2001,
    - the manual of Neuroexplorer software.
    The timebin equal to bin and the max time delay equal to delta,
    which give the rightmost edge until to calculate the ACF.
    The max time delay is equal to delta and it gives the righmost edge 
    until to calculate the ACF. The following check on it wants to assure 
    that results come from quasi-stationary data, so that a too high delay
    doesn't force to work on cutted data.
    For two distinct timestamps we have to add 2, in order to count
    their opposite temporal disposition, as for permutations with repetition.
    """
    
    if not T: return []
    len_T               = len(T)
    max_delay           = int(math.floor(delta / bin))
    _ceil_start_bin     = int(math.ceil(start / bin))
    _ceil_start         = _ceil_start_bin * bin
    start_delay         = int(_ceil_start_bin)
    delay_multiplier    = 1 if two_side else 0
    rightmost_time      = T[-1]
    leftmost_time       = T[0]
    T_duration          = float(rightmost_time - leftmost_time)
    nr_of_used_spikes   = 0
    _floor              = math.floor
    
    if delta == 0:
        if time_normalized: return [len_T / T_duration]
        return [len_T]
    if delta * (delay_multiplier + 1) >= rightmost_time - _ceil_start:
        print "The delay requested is too high, there are no enough data available."
        return []
    
    bins = list(itertools.chain(*(xrange(-max_delay * delay_multiplier, \
                                    -_ceil_start_bin * delay_multiplier), \
                             xrange(_ceil_start_bin, max_delay))))
    acf_dict = dict((b, 0) for b in bins)

    for t in T:
        if rightmost_time - t <= delta: break
        elif two_side and t - leftmost_time <= delta: continue 
        nr_of_used_spikes += 1
        for j in T:
            t_dist = j - t
            if t_dist >= delta: break
            elif not(two_side) and t_dist < 0: continue
            elif -delta < t_dist <= -_ceil_start or \
                 _ceil_start <= t_dist < delta: acf_dict[_floor(t_dist / bin)] += 1
    
    if no_AC_zero: acf_dict[0] -= nr_of_used_spikes

    if not(no_bias): acf = [acf_dict[k] for k in bins]
    else:
        bias = pow(nr_of_used_spikes, 2) * bin / T_duration
        acf = [acf_dict[k] - bias for k in bins]
    
    if bin_normalized: acf = [ac / bin for ac in acf]
    if time_normalized: return [ac / T_duration for ac in acf]
    else: return acf


def acf_OLD(T,
        bin,
        delta, 
        start = 0,
        two_side = True, 
        no_bias = True,
        no_AC0 = True,
        time_normalized = True,
        bin_normalized = False,
        fill_pads = True):

    """We assume that T is already sorted.
    Calculating autocorrelation for timestamps data T, as explained in:
    - P Dayan, LF Abbott. Theoretical Neuroscience (MIT Press). 2001,
    - the manual of Neuroexplorer software.
    The timebin equal to bin and the max time delay equal to delta,
    which give the rightmost edge until to calculate the ACF.
    The max time delay is equal to delta and it gives the righmost edge 
    until to calculate the ACF. The following check on it wants to assure 
    that results come from quasi-stationary data, so that a too high delay
    doesn't force to work on cutted data.
    For two distinct timestamps we have to add 2, in order to count
    their opposite temporal disposition, as for permutations with repetition.
    """
    
    len_T               = len(T)
    max_delay           = int(math.floor(delta / bin))
    _ceil_start_bin     = int(math.ceil(start / bin))
    _ceil_start         = _ceil_start_bin * bin
    start_delay         = int(_ceil_start_bin)
    delay_multiplier    = 1 if two_side else 0
    rightmost_time      = T[-1]
    leftmost_time       = T[0]
    T_duration          = float(rightmost_time - leftmost_time)
    nr_of_used_spikes   = 0
    _floor              = math.floor
    
    if delta == 0:
        if time_normalized: return [len_T / T_duration]
        return [len_T]
    if delta * (delay_multiplier + 1) >= rightmost_time - _ceil_start:
        print "The delay requested is too high, there are no enough data available."
        return []
    
    bins = list(itertools.chain(*(xrange(-max_delay * delay_multiplier, \
                                    -_ceil_start_bin * delay_multiplier), \
                             xrange(_ceil_start_bin, max_delay))))
    acf_dict = dict((b, 0) for b in bins)

    for t in T:
        if rightmost_time - t <= delta: break
        elif two_side and t - leftmost_time <= delta: continue 
        nr_of_used_spikes += 1
        for j in T:
            t_dist = j - t
            if t_dist >= delta: break
            elif not(two_side) and t_dist < 0: continue
            #elif not(t_dist) and no_AC0: continue
            elif -delta < t_dist <= -_ceil_start or \
                 _ceil_start <= t_dist < delta: acf_dict[_floor(t_dist / bin)] += 1
    
    if no_AC0 and not(start): acf_dict[0] -= nr_of_used_spikes

    if not(no_bias): acf = [acf_dict[k] for k in bins]
    else:
        bias = pow(nr_of_used_spikes, 2) * bin / T_duration
        acf = [acf_dict[k] - bias for k in bins]
    
    if bin_normalized: acf = [ac / bin for ac in acf]
    if time_normalized: return [ac / T_duration for ac in acf]
    else: return acf



def acf(T,
        bin,
        delta, 
        start = 0,
        two_side = True,
        no_bias = True,
        no_AC0 = True,
        time_normalized = True,
        bin_normalized = False,
        fill_pads = False):

    """Auto-Correlation Function (ACF) for timestamps in T (assumed already sorted).
    Some of the inputs are:
        bin - timebin
        delta - the time below which calc. ACF.
        start - the time above which calc. ACF.
        two_side - the specualar of ACF at negative delta values
        no_AC0 - if True, ACF(dt=0), shortly AC0, is set to zero 
    Reference:
    - P Dayan, LF Abbott. Theoretical Neuroscience. MIT Press, 2001.
    - usage manual of Neuroexplorer software.
    NB The ACF(dt=0) is treated as a case apart
    """

    delta_bin_floored   = int(delta / bin)
    delta_floored       = bin * delta_bin_floored
    start_bin_ceiled    = int(math.ceil(start / bin))
    start_ceiled        = start_bin_ceiled * bin
    leftmost_time, rightmost_time = T[0], T[-1]
    T_duration          = float(rightmost_time - leftmost_time)
   
    if delta < 0 or start < 0:
        print "Delay and/or start time can't be negative!"
        return []
    elif delta < start:
        print "Delay is < time from which ACF should be calculated"
        return []
    elif delta == 0:
        len_T = len(T)
        if time_normalized: return [len_T / T_duration]
        return [len_T]
    elif delta >= rightmost_time - start_ceiled:
        print "Delay is too high, no enough data available."
        return []

    spikes_nr = 0
    acf = [0] * int(delta_bin_floored + 1 - start_bin_ceiled) # +1 -> ACF(dt=0)
    Tit = iter(T)
    for t in T:
        if rightmost_time - t <= delta_floored: break
        spikes_nr += 1
        Tit, Tit2 = itertools.tee(Tit)
        Tit2.next() # shifting from t to t+1, ACF(dt=0) left out by now
        for j in Tit2:
            dt = j - t
            if start_ceiled <= dt < delta_floored:
                acf[int(dt / bin) + 1 - start_bin_ceiled] += 1; continue
            if dt >= delta_floored: break
        Tit.next()
     
    bias = pow(spikes_nr, 2) * bin / T_duration if no_bias else 0.
    acf_factor_Den = bin if bin_normalized else 1.
    if time_normalized: acf_factor_Den *= T_duration
    acf[1:] = [(ac - bias) / acf_factor_Den for ac in acf[1:]]
    
    if fill_pads: acf = [0] * start_bin_ceiled + acf
    elif not(no_AC0) and not(start): acf[0] = (spikes_nr - bias) / acf_factor_Den
    if two_side: acf = acf[::-1][:-1] + acf
    
    return acf


def ccf(T1, T2,
        bin,
        delta, 
        start = 0,
        no_bias = True,
        time_normalized = True,
        bin_normalized = False,
        fill_pads = False):

    """Cross-Correlation Function (CCF) between timestamps in T1 and T2 (assumed already sorted).
    Some of the inputs are:
        bin - timebin
        delta - the time below which calc. ACF.
        start - the time above which calc. ACF.
    Reference: see references of ACF function above.
    NB The CCF(dt=0) is treated as a case apart
    """
    
    delta_bin_floored   = int(delta / bin)
    delta_floored       = bin * delta_bin_floored
    start_bin_ceiled    = int(math.ceil(start / bin))
    start_ceiled        = start_bin_ceiled * bin
    rightmost_time, leftmost_time = T1[-1], T1[0]
    T_duration          = float(rightmost_time - leftmost_time)
    start_to_delta_bin  = delta_bin_floored - start_bin_ceiled
   
    if delta < 0 or start < 0:
        print "Delay and/or start time can't be negative!"
        return []
    elif delta < start:
        print "Delay is < time from which ACF should be calculated"
        return []
    elif delta >= rightmost_time - start_ceiled:
        print "Delay is too high, no enough data available."
        return []
    
    Tit1, Tit1_pre = itertools.tee(T1)
    for t in Tit1_pre:
        if t - leftmost_time > delta_floored: break
        Tit1.next()
    ccf = [0] * (2 * start_to_delta_bin + 1)
    spikes_nr = 0
    _fl, _cl = math.floor, math.ceil
    Tit = iter(T2)
    for t in Tit1:
        if rightmost_time - t <= delta_floored: break
        spikes_nr += 1
        Tit, Tit2 = itertools.tee(Tit)
        for j in Tit2:
            dt = j - t
            if not(dt): ccf[start_to_delta_bin] += 1; continue
            if -delta_floored <= dt < -start_ceiled:
                ccf[start_to_delta_bin + start_bin_ceiled + int(_fl(dt / bin))] += 1; continue
            if start_ceiled <= dt < delta_floored:
                ccf[int(dt / bin) + 1 + start_to_delta_bin - start_bin_ceiled] += 1; continue
            if dt < -delta_floored: Tit.next(); continue 
            if dt >= delta_floored: break
     
    bias = pow(spikes_nr, 2) * bin / T_duration if no_bias else 0.
    ccf_factor_Den = bin if bin_normalized else 1.
    if time_normalized: ccf_factor_Den *= T_duration
    ccf = [(cc - bias) / ccf_factor_Den for cc in ccf]
    
    if start_bin_ceiled: ccf[start_to_delta_bin] = 0.
    if fill_pads: ccf = ccf[:start_to_delta_bin] + ([0] * start_bin_ceiled) + \
                        [ccf[start_to_delta_bin]] + ([0] * start_bin_ceiled) + \
                        ccf[start_to_delta_bin + 1:]
    
    return ccf


def ccf_OLD(T, P,
            bin, delta, 
            start = 0,
            no_bias = True,
            time_normalized = True,
            bin_normalized = False,
            two_side = True,
            no_CC0 = True,
            fill_pads = False):

    """We assume that T is already sorted.
    Calculating autocorrelation for timestamps data T, as explained in:
    - P Dayan, LF Abbott. Theoretical Neuroscience (MIT Press). 2001,
    - the manual of Neuroexplorer software.
    The timebin equal to bin and the max time delay equal to delta,
    which give the rightmost edge until to calculate the ACF.
    The max time delay is equal to delta and it gives the righmost edge 
    until to calculate the ACF. The following check on it wants to assure 
    that results come from quasi-stationary data, so that a too high delay
    doesn't force to work on cutted data.
    For two distinct timestamps w ehave to add 2, in order to count
    their opposite temporal disposition, as for permutations with repetition.
    """
    
    len_T               = len(T)
    max_delay           = int(math.floor(delta / bin))
    _ceil_start_bin     = int(math.ceil(start / bin))
    _ceil_start         = _ceil_start_bin * bin
    start_delay         = int(_ceil_start_bin)
    delay_multiplier    = 1
    rightmost_time      = T[-1]
    leftmost_time       = T[0]
    T_duration          = float(rightmost_time - leftmost_time)
    if delta * (delay_multiplier + 1) >= rightmost_time - _ceil_start:
        print "The delay requested is too high, there are no enough data available."
        return []
    
    bins = itertools.chain(*(xrange(-max_delay, -_ceil_start_bin), \
                             xrange(_ceil_start_bin, max_delay)))
    ccf_dict = dict((b, 0) for b in bins)
    _floor = math.floor
    nr_of_used_spikes = 0
    for t in T:
        if rightmost_time - t <= delta: break
        if t - leftmost_time <= delta: continue 
        nr_of_used_spikes += 1
        for p in P:
            t_dist = p - t
            if t_dist >= delta: break
            if -delta < t_dist <= -_ceil_start or \
                 _ceil_start <= t_dist < delta: ccf_dict[_floor(t_dist / bin)] += 1
        
    if not(no_bias): ccf = [ccf_dict[k] for k in sorted(ccf_dict.iterkeys())]
    else:
        bias = pow(nr_of_used_spikes, 2) * bin / T_duration
        ccf = [ccf_dict[k] - bias for k in sorted(ccf_dict.iterkeys())]
    
    if bin_normalized: ccf = [cc / bin for cc in ccf]
    if time_normalized: return [cc / T_duration for cc in ccf]
    else: return ccf


def acf_binned(T_binned, bin, delta, start = 0):
    """Calculating autocorrelation for binned timestamps data T,
    where bin is the timebin lenght by which we divided point process data
    obtaining their histogram T_binned. Correlation is here expressed as the
    number of all possible combinations between spikes in bins from the two
    time data series.
    As example the result given at the first time bin is equal to:
        sum([i ** 2 for i in T_binned[:-int(math.ceil(delta / bin))]])
    The max time delay is equal to delta and it gives the rightmost edge 
    until to calculate the ACF. The following check on it wants to assure 
    that results come from quasi-stationary data, so that a too high delay
    doesn't force to work on cutted data.
    """

    corr_function = numpy.correlate
    len_T_binned, max_delay, start_index = len(T_binned), int(math.floor(delta / bin)), int(math.ceil(start / bin))
    if max_delay >= len_T_binned - start_index:
        print "The delay requested is too high, not enough data available."
        return []
    return [corr_function(T_binned[start_index:len_T_binned - max_delay],
                          T_binned[start_index + d :len_T_binned + d - max_delay],
                          mode = 'valid')[0] for d in xrange(max_delay)]


def ccf_binned(X, Y, delta, bin = 1.):
    """Calculating Cross-Correlation between binned timestamps vectors X and Y,
    with bin gives the timebin lenght and delta gives the time delay."""

    ccf = lambda x, y: numpy.correlate(x, y, mode = 'valid', old_behavior = True)[0]
    l, delay = min(len(X), len(Y)), int(math.ceil(delta / bin))
    if delay >= l: raise Exception, "The delay requested is too high, not quite much data available" 
    P, Q = numpy.asarray(X[delay:l - delay]), numpy.asarray(Y)
    ccf_data = numpy.asarray([ccf(P, Q[d:l + d - delay]) for d in xrange(0, 2*delay)])
    #                              for d in xrange(-delay, delay)])
    return tuple((ccf_data - ccf_data.mean()) / (l - 2*delay))


def ac_isi(T, delta, mode = 'valid'):
    """Calculating autocorrelation for isi's from data L,
    with a max delay of delta nr. of steps"""
    
    corr_function = numpy.correlate
    l = len(T)
    if delta >= l: raise Exception, "The delay requested is too high, not quite much data available" 
    return [corr_function(T[:l - delta], T[d:l + d - delta], mode)[0] for d in range(delta + 1)]


def mi_binned(T, bin, delta):
    """Calculating mutual information for binned timestamps data L,
    with timebin lenght equal to bin and a max time delay equal to delta"""

    try:
        import statinfo
        __all__.append('statinfo')
    except ImportError:
        raise ImportError, "Impossible to import statinfo module"

    mi = statinfo.mutualinfo
    l, max_delay = len(T),  int(math.ceil(delta / bin))
    if max_delay >= l: raise Exception, "The delay requested is too high, not quite much data available" 
    if max_delay == 0: return statinfo.entropy(T)
    return [mi(T[:l - max_delay], T[d:l + d - max_delay]) for d in range(max_delay + 1)]


def concatsubseq(L1, L2, cmplen = 2):

    L2cmpseq = tuple(L2[:cmplen])
    L1it = itertools.tee(iter(L1[::-1]), cmplen)
    for i in xrange(1, cmplen):
        for j in xrange(i): L1it[j].next()      
    i = 0     
    for t in itertools.izip(*L1it):
        t = tuple(t)
        if t == L2cmpseq: return L1[:-i-cmplen]+L2
        i += 1
    return L1+L2


def combinations(iterable, r):
    """Return r length subsequences of elements from the input iterable.
    Combinations are emitted in lexicographic sort order. So, if the input iterable is sorted, 
    the combination tuples will be produced in sorted order.
    Elements are treated as unique based on their position, not on their value.
    So if the input elements are unique, there will be no repeat values in each combination.
    New in version 2.6. 
    Examples:
    combinations('ABCD', 2) --> AB AC AD BC BD CD
    combinations(range(4), 3) --> 012 013 023 123

    Code from http://docs.python.org/library/itertools.html#module-itertools
    available in the built-in module itertools of Python 2.6, but not in 2.5.
    It has been slightly changed to cover the case where len(iterable) < r.
    """

    pool = tuple(iterable)
    n = len(pool)
    indices = range(r)
    if n < r: return
    yield tuple(pool[i] for i in indices)
    while True:
        for i in reversed(xrange(r)):
            if indices[i] != i + n - r: break
        else: return
        indices[i] += 1
        for j in xrange(i+1, r): indices[j] = indices[j-1] + 1
        yield tuple(pool[i] for i in indices)


def dotproduct(vec1, vec2):
    """The following code comes from the url:
    http://docs.python.org/library/itertools.html#module-itertools
    _square = lambda i: pow(i, 2)
    return sum(itertools.imap(operator.mul, vec1, vec2))
    But this is not as performant as numpy.inner function
    """

    return numpy.inner(vec1, vec2)



def dotproduct_normalized(vec1, vec2):
    """Code inspired to the one at the url:
    http://docs.python.org/library/itertools.html#module-itertools
    _square = lambda i: pow(i, 2)
    return sum(itertools.imap(operator.mul, vec1, vec2)) / \
    (math.sqrt(sum(itertools.imap(_square, vec1))) * \
    math.sqrt(sum(itertools.imap(_square, vec2))))
    But this is not as performant as numpy.inner and numpy.linalg.norm functions
    """

    return numpy.inner(vec1, vec2) / (numpy.linalg.norm(vec1) * numpy.linalg.norm(vec2))



def convolve_with_pdf(T, 
                      _pdf_arange = [],
                      mode_conv = 'full',
                      pdf = 'norm',
                      pdf_start = -1,
                      pdf_end = 1,
                      pdf_bin = 1,
                      cut_to_inputs_length = True,
                      **kwargs):
    """Calculate the convolution between T1 and a list of nr_samples_pdf from a continuous 
    or discrete probability distribution, given in input as pd. Such pd is picked 
    between those defined in scipy.stats module (scipy 0.6 at the moment, as from
    http://docs.scipy.org/doc/scipy/reference/stats.html).
    In order to cut boundary effects given by _pdf_arange, use mode = 'full' for the
    numpy.convolve function and take the slice of results long up to len(T).
    """
    
    if not len(T): return numpy.array([])
    if _pdf_arange == []:

        _pdf =  eval('scipy.stats.' + pdf)(**kwargs).pdf
        _pdf_arange = _pdf(numpy.arange(pdf_start, pdf_end + pdf_bin, pdf_bin))
        
        ######################################################################
        # Normalizing _pdf_arange in order to make it a discrete prob. distr.  
        ######################################################################
        _inverse_sum_pdf_arange = pow(sum(_pdf_arange), -1) 
        _pdf_arange = map(_inverse_sum_pdf_arange.__mul__, _pdf_arange)

    if not(cut_to_inputs_length) or mode_conv == 'valid': return numpy.convolve(T, _pdf_arange, mode = mode_conv)
    else:  return numpy.convolve(T, _pdf_arange, mode = mode_conv)[:len(T)]



def normalize(M, vars_in_row = True):
    """Normalization consists to perform each given data x in: 
    (x - mean([..., x, ...])) / standard_deviation([..., x, ...])"""
    
    if not(hasattr(M, '__iter__')): raise Exception, 'Object in input is not iterable!'
    
    if not(hasattr(M[0], '__iter__')):
        m, sd = meanval(M), standdev(M)
        if m == 0: return M
        if sd == 0: return [i - m for i in M]
        if type(M) == type(numpy.array([])): return numpy.divide([i - m for i in M], sd)
        else: return [(i - m) / sd for i in M]
    
    if vars_in_row: M_temp, variables_nr = map(None, *M), len(M[0])
    else: M_temp, variables_nr = M, len(M)
    all_mean, all_std = map(meanval, M_temp), map(standdev, M_temp)
    if None in all_mean:
        print 'There is a None mean value!'
        return [[None] * len(M_temp[0])] * len(M_temp)
    if not(all(all_std)):
        print 'There is a zero or None standard deviation value!'
        return [[None] * len(M_temp[0])] * len(M_temp)
    M_normalized_temp = [[(i - m) / sd for i in L] for (L, m, sd) in zip(M_temp, all_mean, all_std)]
    if vars_in_row:
        if variables_nr == 1: return [map(None, *M_normalized_temp)]
        return map(None, *M_normalized_temp)
    else: return M_normalized_temp



def shape_normalize(M, data_sum = None):
    """This function normalizes an array of data by dividing its values for 
    their sum, or by the value data_sum in input if it is not None"""
    
    if not(hasattr(M, '__iter__')): raise Exception, 'There is no an iterable object in input!'
    if not(len(M)): raise IndexError, 'Array is invalid: it is zero-length!'
    if hasattr(M[0], '__iter__'): return map(shape_normalize, M) 
    #raise Exception, 'This function is actually implemented only for arrays!'
    if type(M) == type(numpy.array([])):
        if data_sum is None: data_sum = M.sum(dtype = numpy.float32)
        if data_sum == 0: return M
        return numpy.divide(M, data_sum)
    else:
        if data_sum is None: data_sum = sum(M)
        if data_sum == 0: return M
        return map((1.0 / data_sum).__mul__, M)



def myXrange(i, stop = None, step = 1):
    
    if stop == None:
            while True: yield i; i += step
    while i < stop: yield i; i += step


def purgeISI(ts, minISI = .003):

    ts_it = iter(ts) ; last_t = ts_it.next() ; _ts = [last_t]
    for t in ts_it:
        if t-last_t > minISI: _ts.append(t); last_t = t
    return _ts
    

"""
###########################
# Test
###########################
if __name__ == 'statfuncs':

    import timeit
"""
"""
    # Fano Factor tests
    X = numpy.random.rand(10**7)*100
    t = timeit.Timer("Fano_factor(X)", "from __main__ import Fano_factor, X")
    print t.timeit(number = 10), 's required'
    t = timeit.Timer("Fano_factor(X, 10, 90)", "from __main__ import Fano_factor, X")
    print t.timeit(number = 10), 's required'
    t = timeit.Timer("Fano_factor(X, 10)", "from __main__ import Fano_factor, X")
    print t.timeit(number = 10), 's required'
    t = timeit.Timer("Fano_factor(X, upper_threshold = 90)", "from __main__ import Fano_factor, X")
    print t.timeit(number = 10), 's required'
    print Fano_factor(X)    
"""

"""
    def testacf(T, new = False, bin = 0.001, delta = 0.01, start = 0.005, i = 1):
        if i == 0: kw = {}
        if i == 1: kw = {'two_side': True, 'no_AC0': True}
        if i == 2: kw = {'two_side': True, 'no_AC0': False}
        if i == 3: kw = {'two_side': True, 'no_AC0': True}
        if i == 4: kw = {'no_AC0': True, 'no_bias': False}
        if i == 5: kw = {'two_side': True, 'no_AC0': True, 'no_bias': True}
        if i == 6: kw = {'two_side': True, 'no_AC0': True, 'no_bias': False, 'normalized': False}
        if i == 7: kw = {'two_side': False, 'no_AC0': True, 'no_bias': False, 'normalized': False}
        if i == 8: kw = {'two_side': True, 'no_AC0': True, 'no_bias': True, 'start': start}
        if i == 9: kw = {'two_side': True, 'no_AC0': True, 'no_bias': True, 'start': start}
        if i == 10: kw = {'two_side': True, 'no_AC0': True, 'no_bias': True, 'start': start, 'fill_pads': True}
        if new: return acf(T, bin, delta, **kw)
        return acf_OLD(T, bin, delta, **kw)
    
    def testccf(T1, T2, new = False, bin = 0.001, delta = 0.01, start = 0.005, i = 1):
        if i == 1: kw = {}
        if i == 2: kw = {'start': start}
        if i == 3: kw = {'start': start, 'fill_pads': True}
        if new: return ccf(T1, T2, bin, delta, **kw)
        return ccf_OLD(T1, T2, bin, delta, **kw)
    
    try:    
        import pylab
        import numpy
        import timeit
    
        T = numpy.random.rand(1000) * 100; T.sort()
        
        ###################
        # Test the ACF
        ###################
        pylab.figure(1)
        t = timeit.Timer("testacf(T, i=1)", "from __main__ import testacf, T")
        print t.timeit(number = 1), 's required'
        pylab.plot(testacf(T, i = 9))

        pylab.figure(2)
        t = timeit.Timer("testacf(T, new = True, i=1)", "from __main__ import testacf, T")
        print t.timeit(number = 1), 's required'
        pylab.plot(testacf(T, new = True, i = 9))

        pylab.show()
"""
"""
        ###################
        # Test the CCF
        ###################
        T2 = numpy.random.rand(1000) * 100; T2.sort()
        #T = numpy.arange(1000); T2 = (T + .0075).tolist() + (T + .0025).tolist(); T2.sort()

        pylab.figure(1)
        t = timeit.Timer("testccf(T, T2, i = 1)", "from __main__ import testccf, T, T2")
        print t.timeit(number = 1), 's required'
        w = testccf(T, T2, i = 2)
        pylab.plot(w)

        pylab.figure(2)
        t = timeit.Timer("testccf(T, T2, new = True, i = 1)", "from __main__ import testccf, T, T2")
        print t.timeit(number = 1), 's required'
        w = testccf(T, T2, new = True, i = 2)
        pylab.plot(w)

        pylab.show()
    
    except ImportError: raise

"""
"""
    T = numpy.arange(0, 10, 0.01)
    Y_bin = [0] * len(T)
    for k in xrange(100): Y_bin[random.randint(0, len(T) - 1)] = 1 
    Y_bin[0] = 1
    Y_conv = convolve_with_pdf(Y_bin, pdf = 'expon', pdf_start = 0, \
                                pdf_end = 1, pdf_bin = 0.01)
    pylab.plot(T, Y_bin, 'x')
    pylab.plot(T, Y_conv, 'c')
    #print 'Sum of convoluted data =', sum(Y_conv), ' and sum of binned data =', sum(Y_bin)
    
    #A = [0.1, 0.3, 0.4, 0.66, 0.9, 1.0, 1.3, 1.9, 2.33, 3.5, 4.0, 5.1]
    #print acf_binned(A, 0.05, 2).AC
    #print ac_isi(A, 8).AC
    colors = ['b',  # blue
              'g',  # green
              'r',  # red
              'c',  # cyan
              'm',  # magenta
              'y',  # yellow
              'k',  # black
              'w']  # white
    
    #T = [[random.randint(0, 25) for i in xrange(1000)] for j in range(8)]
    T = [[random.randint(0, 25) for i in xrange(1000)] for j in range(50)]
    M = [mi_binned(Z, 5.0 / 1000, 0.5) for Z in T]
    pylab.figure(figsize = (5.6, 3.9))
    #for MI in M: pylab.plot(MI, marker = ',', color = colors[M.index(MI)])
    for MI in M: pylab.plot([p + 0.1 * M.index(MI) for p in MI if MI.index(p) > 0], \
                            color = '0.3', linewidth = 0.3, linestyle = '-')
    
    pylab.show()
"""
