In [ ]:
# this function code has to be debugged further with matrix dimension matching

import numpy as np

def calculate_occurance_time(frequency, numFramesPerStim, spikeMat, binaryWords, numCoactive):
    

    # check occurence of high-order patterns
    idx = np.nonzeros(frequency >=2) # the frequency of spikings
    
    timeInMovieSTD = np.zeros((idx.shape[0],1))

    binEdges = range(-20,21)
    timeJitterDistribution = np.zeros((1, len(binEdges)))

    trialDist = []

    interval = []
    timeAll = []

    timeInMovieDist = np.zeros(numFramesPerStim, 1)
    occurance_time = np.zeros(len(frequency), 1)
    
    for i in range(len(idx)): # search through the neurons that fired at least twice
        # find out the time points when there are at least numCoactive neurons firing together
        times = np.nonzeros(spikeMat[binaryWords[i,:],:].sum()>=numCoactive)
        
        # kl: find out which the exact time frame the cell assembly respond to
        timeInMovie = times % numFramesPerStim
        # kl: the response at 0 is assumed to be responding to the last movie frame
        timeInMovie[timeInMovie==0] = numFramesPerStim
        # distribution of the time points when coactive firing
        timeInMovieDist[timeInMovie] = timeInMovieDist[timeInMovie]+1
        # what trials includes the coactive firings
        trial = (times-timeInMovie)/numFramesPerStim + 1
        # average the time points of coactive firings for all the frequency
        occurance_time[idx[i]] = np.mean(timeInMovie)
        # histogram of the coactive firings around average time points within binEdges
        counts = np.histogram(timeInMovie-np.mean(timeInMovie), range=binEdges)
        # normalize counts
        if counts.sum()>1 : 
            counts = counts/counts.sum()
        
        # the histogram of distribution of firing time points of coactive assembly around average time points
        timeJitterDistribution = timeJitterDistribution + counts
        timeInMovieSTD[i] = np.std(timeInMovie)

   
    # the histogram of distribution of firing time points of coactive assembly around average time points
    timeJitterDistribution = timeJitterDistribution/len(idx)