#!/usr/bin/python #----------------------------------------------------------------------- # File : coconad.py # Contents: COntinuous time ClOsed Neuron Assembly Detection # Author : Christian Borgelt # History : 2012.02.09 file created (with hyperedge finding only) # 2012.02.13 functions for finding matchings added # 2012.02.14 variable bounding function for branch and bound # 2012.02.18 structure of backtracking recursion improved # 2012.02.19 exploited precomputed hyperedge intersections # 2012.10.08 merged with hedges.py, train functions added # 2012.10.09 main function coconad() and its recursion added # 2012.10.10 main program with option evaluation added # 2013.01.24 program option evaluation improved # 2013.02.05 perfect extension pruning removed (incorrect) # 2013.02.06 function get_hedges() adapted to separate trains # 2013.02.07 all support variants made accessible via options # 2013.02.08 interface of function coconad() changed to trains # 2013.02.10 (safeguarded) perfect extension pruning added # 2013.02.11 item set reporting made a separate function # 2013.02.13 minimum support check added to function coconad() # 2013.02.14 replaced "while True" with "while 1" (faster) # 2013.02.17 trains lengths retrieved into variables (faster) # 2013.02.19 removed point and train filtering from pexcube() # 2013.02.21 bug in deleting added elements of elim[] fixed # 2013.02.22 interface made compatible with C library version # 2013.05.22 all frequent and maximal item set output added # 2013.11.22 input format with one train per line added # 2013.12.06 collection and output of pattern spectrum added # 2014.01.29 file reading made more general (separators etc.) # 2014.07.19 mining item sequences added (respect item order) # 2014.10.13 filtering with a support border added # 2014.10.28 influence map overlap based support added # 2014.12.12 better version of overlap based support added # 2014.12.15 pattern spectrum reporting corrected/extended # 2016.03.22 perfect extension flag changed to mode string #----------------------------------------------------------------------- from sys import argv, stderr, float_info, version_info from time import time #----------------------------------------------------------------------- # Constants #----------------------------------------------------------------------- oo = float('inf') # positive infinity as a constant huge = float_info.max # maximum float value as a constant heapmin = 7 # minimum number of trains for a heap #----------------------------------------------------------------------- # Train and Interval List Functions #----------------------------------------------------------------------- def train2ivals (train, width=0.003): '''train2ivals (train, width=0.003) Compute an interval list from a spike train and a window width. train (sorted) list of points (in time) width width of time window or influence region returns a (sorted) list of (disjoint) intervals [start,end]''' if not train: return [] # check for an empty train width *= 0.5 # compute half the window width out = [] # initialize the output list s = train[0] # get the first point a = s-width # compute start and end b = s+width # of the first interval for s in train[1:]: # traverse the points in the train c = s-width # compute start of new interval if c > b: # if the intervals are disjoint, out.append((a,b)) # store the preceding interval a = c # and start a new interval b = s+width # update end of current interval out.append((a,b)) # store the last interval return out # return the computed intervals #----------------------------------------------------------------------- def ivfilt (train, ivals): '''ivfilt (train, ivals) Filter a train with a time interval list. train (sorted) list of points/times, i.e. a train ivals (sorted) list of (disjoint) intervals, as pairs [start,end] returns a filtered train, with points/times in the intervals''' if not train or not ivals: # if train or interval list is empty, return [] # the filtered train is empty out = [] # initialize the output train and s = train[0]; i = 0; n = len(train) # get first point/time a,b = ivals[0]; k = 0; m = len(ivals) # and first interval while 1: # train filtering loop if s > b: # if time is after current interval k += 1 # skip the processed interval if k >= m: return out a,b = ivals[k] # get boundaries of the next interval else: # if time may be in current interval if s >= a: # if time is in current interval, out.append(s) # add it to the output list i += 1 # skip the processed point/time if i >= n: return out s = train[i] # get the next point/time #----------------------------------------------------------------------- def pfisect (ivals, train, width=0.003): '''pfisect (ivals, train, width=0.003) Intersect a time interval list with a train for point filtering. ivals (sorted) list of (disjoint) intervals, as pairs [start,end] train (sorted) list of points/times, i.e. a train width width of time window or influence region (width >= 0) returns a (sorted) list of (disjoint) intervals, resulting from an intersection of the given interval list with the intervals induced by the times of the train, i.e. [s-width,s+width]''' if not train or not ivals: # if train or interval list is empty, return [] # the intersection is empty out = []; o = [-oo,-oo] # initialize the output interval list s = train[0]; i = 0; n = len(train) # get first point/time a,b = ivals[0]; k = 0; m = len(ivals) # and first interval while 1: # interval list filtering loop if s > b: # if time is after current interval k += 1 # skip the processed interval if k >= m: return out a,b = ivals[k] # get boundaries of the next interval else: # if time may be in current interval if s >= a: # if time is in current interval, x = a if a > s-width else s-width # intersect y = b if b < s+width else s+width # intervals if x <= o[1]: o[1] = y else: o = [x,y]; out.append(o) i += 1 # skip the processed point/time if i >= n: return out s = train[i] # get the next point/time #----------------------------------------------------------------------- def isect (ivs1, ivs2): '''isect (ivs1, ivs2) Intersect a time interval list with a train for overlap based support. ivs1 (sorted) list of (disjoint) intervals, as pairs [start,end] ivs2 (sorted) list of (disjoint) intervals, as pairs [start,end] returns a (sorted) list of (disjoint) intervals, resulting from an intersection of the given interval lists''' if not ivs1 or not ivs2: # if one interval list is empty, return [] # the intersection is empty out = []; o = [-oo,-oo] # initialize the output interval list a,b = ivs1[0]; k = 0; m = len(ivs1) # get first interval x,y = ivs2[0]; i = 0; n = len(ivs2) # from each of the lists while 1: # intervals intersection loop if a < y and b > x: # if intersection is not empty s = a if a > x else x # compute the intersection t = b if b < y else y # of the two intervals if s <= o[1]: o[1] = t # add intersection to the output else: o = [s,t]; out.append(o) if b < y: # if first interval ends earlier k += 1 # skip the processed interval if k >= m: return out a,b = ivs1[k] # get boundaries of the next interval else: # if second interval ends earlier i += 1 # skip the processed interval if i >= n: return out x,y = ivs2[i] # get boundaries of the next interval #----------------------------------------------------------------------- # Binary Support Functions #----------------------------------------------------------------------- def sift (heap, lft, rgt): '''sift (heap, lft, rgt) Let an element sift down in a train heap. heap list of pairs (list of points/times, position) (a heap element may have additional entries, though) lft left border of the heap (index of sift element) rgt right border of the heap (index of last element)''' h = heap[lft] # note the sift element t,k = h[:2]; s = t[k] # and get its current time i = lft +lft +1 # compute index of first successor while i <= rgt: # sift loop t,k = heap[i][:2] # get the first successor x = t[k] # and from it the current point if i < rgt: # if there is a second successor t,k = heap[i+1][:2] # and it is smaller than the first y = t[k] # go to the second successor if y < x: x = y; i += 1 if x >= s: # if the successor is no less break # than the sift element, abort heap[lft] = heap[i] # let the successor ascend in heap lft = i # note index of current element and i += i+1 # compute index of first successor heap[lft] = h # store the sift element #----------------------------------------------------------------------- def support (trains, width=0.003): '''support (trains, width=0.003) Compute support of a list of trains. trains list of pairs (item id, list of points/times) width width of time window or influence region (width >= 0) returns the support of the item set underlying the trains (number of coincident events)''' k = len(trains) # get the number of trains if k <= 1: # if there is only one train, return len(trains[0][1])# simply return its length supp = 0 # initialize the support counter if k <= 2: # if there are only two trains s,t = [t for n,t in trains] n,m = len(s),len(t) # get the trains and their lengths i = k = 0 # traverse the two trains and while i < n and k < m: # count synchronous points/times if abs(s[i]-t[k]) <= width: i += 1; k += 1; supp += 1 elif s[i] < t[k]: i += 1 else: k += 1 return supp # return the computed support heap = [[t,0] for n,t in trains if len(t) > 0] n = len(heap) # build a heap/array of all trains if n < k: # and check for at least one train return 0 # in each of the trains if n >= heapmin: # if there are many trains for i in reversed(range(n >> 1)): sift(heap, i, n-1) # build a heap w.r.t. the first times tmax = max(t[0] for t,i in heap) while 1: # support computation loop t,i = heap[0]; tmin = t[i]; k = 0 if n < heapmin: # if there are only few trains for j in range(1,n):# find the train t,i = heap[j] # with the smallest time if t[i] < tmin: tmin = t[i]; k = j if tmax-tmin > width: # if the times are not synchronous, heap[k][1] += 1 # skip the processed time and t,i = heap[k] # check whether it was the last if i >= len(t): return supp if t[i] > tmax: tmax = t[i] # update the maximum time if n >= heapmin: sift(heap, 0, n-1) continue # let the new time sift down in heap supp += 1 # count the synchronous group for h in heap: # traverse all times in the array h[1] += 1 # skip the processed time and t,i = h # check whether it was the last if i >= len(t): return supp if t[i] > tmax: tmax = t[i] # update the maximum time if n >= heapmin: # if needed, rebuild the heap for i in reversed(range(n >> 1)): sift(heap, i, n-1) #----------------------------------------------------------------------- def seqsupp (trains, width=0.003): '''seqsupp (trains, width=0.003) Compute support of an ordered list of trains. trains list of pairs (item id, list of points/times) width width of time window/maximal distance (width >= 0) returns the coincident events of the item set underlying the trains''' k = len(trains) # get the number of trains if k <= 1: # if there is only one train, return len(trains[0][1])# simply return its length supp = 0 # initialize the support counter if k <= 2: # if there are only two trains s,t = [t for n,t in trains] n,m = len(s),len(t) # get the trains and their lengths i = k = 0 # traverse the two trains and while i < n and k < m: # while not at end of either train d = t[k]-s[i] # count synchronous points/times if d < 0: k += 1 elif d > width: i += 1 else: i += 1; k += 1; supp += 1 return supp # return the computed support curr = [[t,0,r] for r,t in trains if len(t) > 0] n = len(curr) # build a position array of all trains if n < k: return supp # and check for no empty trains a,z = curr[0],curr[-1] # get first and last train/position while a[1] < len(a[0]): # while not at end of first train for i in range(1,n): # traverse the following trains b,c = curr[i-1:i+1] # (or rather consecutive pairs) while c[1] < len(c[0]) and c[0][c[1]] < b[0][b[1]]: c[1] += 1 # ensure item succession in time if c[1] >= len(c[0]): return supp # if at the end of a train, abort if z[0][z[1]] -a[0][a[1]] <= width: supp += 1 # count the synchronous group for i in range(1,n): curr[i][1] += 1 # skip the counted points a[1] += 1 # get the next starting point return supp # return the computed support #----------------------------------------------------------------------- def ovlsupp (trains, width=0.003): '''ovlsupp (trains, width=0.003) Compute support of a list of trains based on influence map overlap. trains list of pairs (item id, list of points/times) width width of time window or influence region (width >= 0) returns the support of the item set underlying the trains''' if not trains: return 0 # check for at least one train w = 0.5*width # get half the window width and supp = 0 # initialize the support counter k = len(trains) # get the number of trains if k <= 1: # if there is only one train train = trains[0][1] # get this single train s = train[0] # get the first point of the train a = s-w; b = s+w # compute start/end of 1st interval for s in train[1:]: # traverse the points in the train c = s-w # compute start of new interval if c > b: # if the intervals are disjoint, supp += b-a # sum length of preceding interval a = c # and start a new interval b = s+w # update end of current interval return (supp +(b-a)) /float(width) heap = [[t,0] for n,t in trains if len(t) > 0] n = len(heap) # build a heap/array of all trains if n < k: # and check for at least one train return 0 # in each of the trains if n >= heapmin: # if there are many trains for i in reversed(range(n >> 1)): sift(heap, i, n-1) # build a heap w.r.t. the first times tmax = max(t[0] for t,i in heap) a = b = -huge # get max. point and init. 1st interval while 1: # support computation loop t,i = heap[0]; tmin = t[i]; k = 0 if n < heapmin: # if there are only few trains for j in range(1,n):# find the train t,i = heap[j] # with the smallest time if t[i] < tmin: tmin = t[i]; k = j if tmax-tmin < width: # if the times are synchronous c = tmax-w # compute start of new interval if c > b: # if the intervals are disjoint supp += b-a # sum length of preceding interval a = c # and start a new interval b = tmin+w # update end of current interval heap[k][1] += 1 # skip the processed point and t,i = heap[k] # check whether it was the last if i >= len(t): break # if yes, abort the loop if t[i] > tmax: tmax = t[i] # update the maximum time if n >= heapmin: sift(heap, 0, n-1) return (supp +(b-a)) /float(width) #----------------------------------------------------------------------- # Graded Support Functions #----------------------------------------------------------------------- def get_hedges (trains, width=0.003): '''get_hedges (trains, width=0.003) Collect hyperedges (overlap groups) for a given item set and a given influence region width from parallel train data. trains list of pairs (item id, list of points/times) width width of time window or influence region (width >= 0) returns a list of hyperedges, each a set of pairs (item id, time)''' k = len(trains) # get the number of trains if k <= 1: # if there is only one train return [(frozenset([s]),1.0) for s in trains[0]] hedges = [] # initialize the list of hyperedges if k <= 2: # if there are only two trains a,b = [n for n,t in trains] s,t = [t for n,t in trains] n,m = len(s),len(t) # get the trains and their lengths i = k = 0 # traverse the two trains while i < n and k < m: # while there is another time in both if s[i] < t[k]: # if smaller time in first train p = (a,s[i]) # get time common to the hyperedges for j in range(k,m): if t[j]-s[i] > width: break hedges.append(frozenset([p,(b,t[j])])) i += 1 # build hyperedges, skip processed time else: # if smaller time in second train p = (b,t[k]) # get time common to the hyperedges for j in range(i,n): if s[j]-t[k] > width: break hedges.append(frozenset([(a,s[j]),p])) k += 1 # build hyperedges, skip processed time return hedges # return the collected hyperedges heap = [[t,0,n] for n,t in trains if len(t) > 0] n = len(heap) # build an array of all trains if n < k: # and check for at least one time return [] # in each of the trains if n >= heapmin: # if there are many trains for i in reversed(range(n >> 1)): sift(heap, i, n-1) # build a heap w.r.t. the first times tmax = max(t[0] for t,i,k in heap) while 1: # hyperedge collection loop t,i = heap[0][:2]; tmin = t[i]; k = 0 if n < heapmin: # if there are only few trains for j in range(1,n): # find the train t,i = heap[j][:2] # with the smallest time if t[i] < tmin: tmin = t[i]; k = j if tmax-tmin <= width: # if there are synchronous times a = heap[k][2] # get item with smallest time h = [[(a,tmin)]] # initialize the hyperedge list for t,i,b in heap: # traverse the (other) trains if b == a: continue w = [] # collect times in the window for j in range(i,len(t)): if t[j]-tmin > width: break w.append((b,t[j])) h = [r+[s] for r in h for s in w] hedges += [frozenset(r) for r in h] heap[k][1] += 1 # skip the processed time and t,i = heap[k][:2] # check whether it was the last if i >= len(t): return hedges if t[i] > tmax: tmax = t[i] # update the maximum time if n >= heapmin: # if there are many trains (heap), sift(heap, 0, n-1) # let the new time sift down in heap #----------------------------------------------------------------------- def idwgt (x): '''idwgt (x) Identify function (simply returns its argument). x influence region overlap measure; 0 <= x <= 1 returns its argument, that is, x''' return x #----------------------------------------------------------------------- def epswgt (x, eps): '''epswgt (x, eps) Epsilon-insensitive weight function. x influence region overlap measure; 0 <= x <= 1 eps tolerance for which 1.0 is returned returns 1 for x >= 1-eps and x/(1-eps) otherwise''' return 1.0 if x >= 1-eps else x/(1-eps) if eps < 1 else 0.0 #----------------------------------------------------------------------- def binwgt (x): '''binwgt (x) Binary weight function (returns 1.0 for non-negative arguments and 0.0 otherwise). x influence region overlap measure; 0 <= x <= 1 returns 1 for x >= 0 and 0 otherwise''' return 1.0 if x >= 0 else 0.0 #----------------------------------------------------------------------- def wgt_hedges (hedges, width=0.003, wgtfn=idwgt): '''wgt_hedges (hedges, width=0.003, wgtfn=idwgt) Weight a set of hyperedges (overlap groups) for a given time window or influence region width and weighting function. hedges list of hyperedges, as lists of pairs (item id, time) width width of time window or influence region (width >= 0) wgtfn weighting function (e.g. idwgt or binwgt) returns a list of weighted hyperedges, as lists of pairs, each with a list of pairs (item id, time) and a weight''' out = [] # initialize the output list for h in hedges: # traverse the collected hyperedges t = [t for i,t in h] # compute the influence region overlap t = 1.0 if width <= 0 else (width-(max(t)-min(t))) /float(width) out.append((h, wgtfn(t))) return out # return the weighted hyperedges #----------------------------------------------------------------------- def elem (s): '''elem (s) Get an element from a set (slightly faster than iter(s).next()). s set from which to get an element returns an element of the set or "None" if the set is empty''' for e in s: return e # return an arbitrary element return None # return 'None' for empty set #----------------------------------------------------------------------- def ccrep (reps, n): '''ccrep (reps, n) Find representative of a node for connected components. reps dictionary of node representatives n node for which to find the representative returns the node representing the connected component that contains the given node n''' r = reps[n] # get representative of the node while r != n: # follow links to representatives n = r; r = reps[n] # until the node represents itself return r # return the found representative #----------------------------------------------------------------------- def get_ccomps (hedges): '''get_ccomps (hedges) Split a list of hyperedges into connected components. hedges list of weighted hyperedges (pairs of a node set and a weight) returns a list of connected components, each of which is a list of hyperedges (pairs of a node set and an overlap weight)''' reps = dict([(x,x) for x in set([n for h in hedges for n in h])]) for h in hedges: # find representatives of components r = [ccrep(reps, n) for n in h] for n in r: reps[n] = r[0] for n in reps: # flatten references to representatives reps[n] = ccrep(reps, n) ccomps = dict([(n,[]) for n in reps if reps[n] == n]) for h in hedges: # collect hyperegdes by con. component ccomps[ccrep(reps, elem(h))].append(h) return [ccomps[n] for n in ccomps] #----------------------------------------------------------------------- def match_greedy (hedges): '''match_greedy (hedges) Find a maximum weight k-partite matching in a greedy fashion. hedges list of weighted k-edges (pairs of a node set and a weight) returns a pair consisting of the list of k-edges in the found matching and the total weight of this matching''' match = [] # while there are selectable edges left, while len(hedges) > 0: # select one with the largest weight h = max(hedges, key = lambda k: k[1]) match.append(h) # exclude the overlapping edges hedges = [e for e in hedges if not e[0] & h[0]] return (match, sum([w for h,w in match])) # return selected edges #----------------------------------------------------------------------- def rec_btrack (best, match, hedges, nolap): '''rec_btrack (best, match, hedges, nolap) Recursively find a maximum weight k-partite matching (with simple backtracking, that is, full enumeration of all matchings). best best solution (best matching) found up to now match already selected k-edges (current matching) hedges remaining, selectable k-edges (extension candidates) nolap dictionary mapping k-edges to sets of overlapping k-edges returns the best matching found together with the total weight''' for i in range(len(hedges)):# traverse the remaining k-edges x = hedges[i] # include k-edge and exclude overlap r = [h for h in hedges[i+1:] if h in nolap[x]] w = match[1] +x[1] # compute the new matching weight if len(r) > 0: # if k-edges left, go into recursion best = rec_btrack(best, (match[0] +[x], w), r, nolap) elif w > best[1]: # if no k-edges left, update solution best = (match[0] +[x], w) return best # return the best matching found #----------------------------------------------------------------------- def match_btrack (hedges): '''match_btrack (hedges) Find a maximum weight k-partite matching with backtracking (full enumeration of all matchings, no branch and bound). hedges list of weighted k-edges (selectable hyperedges) returns a pair consisting of the list of k-edges in the found matching and the total weight of this matching''' n = len(hedges) # get the number of hyperedges and if n <= 0: return ([], 0) # handle trivial cases directly if n <= 1: return (hedges, hedges[0][1]) hedges = sorted(hedges, key = lambda k: k[1], reverse = True) nolap = dict([(h, set([e for e in hedges if not h[0] & e[0]])) for h in hedges]) # recursively find matching return rec_btrack(match_greedy(hedges), ([],0), hedges, nolap) #----------------------------------------------------------------------- def bound_num (hedges): '''bound_num (hedges) Find an upper bound for the total weight of a matching, based on the maximum number of possible edges (with maximum weight). hedges list of weighted k-edges (selectable hyperedges) returns an upper bound for the weight of a matching''' nodes = set([n for h,w in hedges for n in h]) n = int(len(nodes)/len(hedges[0][0])) # max. number of k-edges return sum([w for h,w in hedges[:n]]) # sum n highest weights #----------------------------------------------------------------------- def bound_wgt (hedges): '''bound_wgt (hedges) Find an upper bound for the total weight of a matching, based on the maximum edge weight contribution per node. hedges list of weighted k-edges (selectable hyperedges) returns an upper bound for the weight of a matching''' nodes = set([n for h,w in hedges for n in h]) maxws = [max(w for h,w in hedges if n in h) for n in nodes] return sum(maxws) /len(hedges[0][0]) #----------------------------------------------------------------------- def rec_bandb (best, match, hedges, nolap, bound): '''rec_bandb (best, match, hedges, nolap, bound) Recursively find a maximum weight k-partite matching (with branch and bound, using the given bound estimator). best best solution (best matching) found up to now match already selected k-edges (current matching) hedges remaining, selectable k-edges nolap dictionary mapping k-edges to sets of overlapping k-edges bound bounding function (yields upper bound for matching weight) returns the best matching found together with the total weight''' n = len(hedges) # get the number of hyperedges if n >= 4: # if at least four k-edges are left if match[1] +bound(hedges) <= best[1]: return best # bound result with bounding function for i in range(n): # traverse the selectable k-edges x = hedges[i] # include k-edge and exclude overlap # r = [h for h in hedges[i+1:] if not h[0] & x[0]] r = [h for h in hedges[i+1:] if h in nolap[x]] w = match[1]+x[1] # compute the new matching weight if len(r) > 0: # if k-edges left, go into recursion best = rec_bandb(best, (match[0]+[x], w), r, nolap, bound) elif w > best[1]: # if no k-edges left, update solution best = (match[0]+[x], w) return best # return the best matching found #----------------------------------------------------------------------- def match_bandb (hedges, bound): '''match_bandb (hedges, bound) Find a maximum weight k-partite matching with branch and bound. hedges list of weighted k-edges (pairs of a node set and a weight) bound either 'num' or 'wgt' for bounding functions that guarantee finding the optimal solution, or a number (also accepted as a string) by which the greedy solution is multiplied to obtain a bound for the weight of a matching (heuristic) returns a pair consisting of the list of k-edges in the found matching and the total weight of this matching''' n = len(hedges) # get the number of hyperedges and if n <= 0: return ([], 0) # handle trivial cases directly if n <= 1: return (hedges, hedges[0][1]) if bound == 'wgt': bfn = bound_wgt elif bound == 'num': bfn = bound_num else: # get the bounding function bound = float(bound) # (greedy bound is a heuristic) bfn = lambda h: bound *match_greedy(h)[1] hedges = sorted(hedges, key = lambda k: k[1], reverse = True) nolap = dict([(h, set([e for e in hedges if not h[0] & e[0]])) for h in hedges]) # recursively find matching return rec_bandb(match_greedy(hedges), ([], 0), hedges, nolap, bfn) #----------------------------------------------------------------------- def matching (hedges, algo='greedy', bound='num'): '''matching (hedges, algo='greedy', bound='num') Find a maximum weight k-partite matching. hedges list of weighted k-edges (pairs of a node set and a weight) algo algorithm to use; one of 'greedy', 'btrack', or 'bandb' bound bounding function indicator for the branch and bound search; either 'num' or 'wgt' for bounding functions that guarantee finding the optimal solution, or a number (also accepted as a string) by which the greedy solution is multiplied to obtain a bound for the weight of a matching (heuristic) returns a pair consisting of the list of k-edges in the found matching and the total weight of this matching''' if algo == 'greedy': return match_greedy(hedges) if algo == 'btrack': return match_btrack(hedges) if algo == 'bandb' : return match_bandb (hedges, bound) return None # select function according to algorithm #----------------------------------------------------------------------- def match_ok (hedges): '''match_ok (hedges) Check a computed matching, that is, test whether all hyperedges have pairwise disjoint node sets. hedges list of weighted hyperedges (pairs of a node set and a weight) returns True if hyperedges are pairwise disjoint, False otherwise''' for i in range(len(hedges)-1): r = hedges[i][0] # traverse all hyperedge pairs for h in hedges[i+1:]: # if two hyperedges overlap, no matching if r & h[0]: return False return True # otherwise hyperedges are a matching #----------------------------------------------------------------------- def kpmsupp (trains, width=0.003, wgtfn=idwgt, algo='greedy', bound='num'): '''kpmsupp (trains, width=0.003, wgtfn=idwgt, algo='greedy', bound='num') Compute support of an item set with k-partite matching (support is the weight of a maximum weight k-partite matching). trains list of pairs (item id, list of points/times) width width of time window or influence region (width >= 0) wgtfn weighting function (for example, idwgt or binwgt) algo algorithm to use; one of 'greedy', 'btrack', or 'bandb' bound bounding function indicator for the branch and bound search; either 'num' or 'wgt' for bounding functions that guarantee finding the optimal solution, or a number (also accepted as a string) by which the greedy solution is multiplied to obtain a bound for the weight of a matching (heuristic) returns the support of the item set''' if len(trains) <= 1: # if there is only one train, return len(trains[0][1])# simply return its length hedges = get_hedges(trains, width) ccomps = get_ccomps(hedges) # get connected comps. of hyperedges ccomps = [wgt_hedges(c, width, wgtfn) for c in ccomps] match = [h for c in ccomps for h in matching(c, algo, bound)[0]] return sum([h[1] for h in match]) # find a max. weight matching #----------------------------------------------------------------------- # Main CoCoNAD Functions #----------------------------------------------------------------------- def report (trains, elim, supp, data): '''report (trains, elim, supp, data) Check and report a (closed/maximal) frequent item set. trains list of trains of collected items, each train a pair (item id, list of points/times) elim list of trains of eliminated items, each train a pair (item id, list of points/times) (needed to check whether the item set is closed or maximal) supp support of the item set to report data recursion data, that is, a 9-tuple/list with: target type of frequent item set to find (in ['a','c','m']) width width of time window or influence region (width >= 0) supp minimum support (number of synchronous spiking events) suppfn support function (e.g. support or kpmsupp) zmin minimum size of an item set zmax maximum size of an item set result list/dictionary in which the result is stored report reporting mode (pattern information/pattern spectrum) border support border (minimum per item set size)''' if data[0] != 'a': # if closed or maximal item sets s = supp if data[0] == 'c' else data[2] for e in reversed([e for e in elim if len(e[1]) > s-1]): if data[3](trains+[e], data[1]) >= s: return # check for a closed/maximal item set n = len(trains) # report the (closed/maximal) item set if n < data[4] or n > data[5]: return if data[8] and n < len(data[8]) and supp <= data[8][n]: return # check against sizes and support border if data[7] == 'a': # if to collect the found pattern iset = tuple([i for i,t in trains]) data[6].append((iset,supp)) elif data[7] in '+-': # if map from sizes to support ranges if n not in data[6]: data[6][n] = [supp,supp] elif supp < data[6][n][0]: data[6][n][0] = supp elif supp > data[6][n][1]: data[6][n][1] = supp else: # if to collect a pattern spectrum s = (n,supp) # as a map from signatures to counts if s in data[6]: data[6][s] += 1 # update the else: data[6][s] = 1 # pattern spectrum #----------------------------------------------------------------------- def seqreport (trains, exts, supp, data): '''seqreport (trains, exts, supp, data) Check and report a (closed/maximal) frequent item sequence. trains list of trains of collected items, each train a pair (item id, list of points/times) exts list of trains of possible extension items, each train a pair (item id, list of points/times) (needed to check whether the item set is closed or maximal) supp support of the item set to report data recursion data, that is, a 9-tuple/list with: target type of frequent item set to find (in ['a','c','m']) width width of time window or influence region (width >= 0) supp minimum support (number of synchronous spiking events) suppfn support function (e.g. support or kpmsupp) zmin minimum size of an item set zmax maximum size of an item set result list/dictionary in which the result is stored report reporting mode (pattern information/pattern spectrum) border support border (minimum per item set size)''' if data[0] != 'a': # if closed or maximal item sets s = supp if data[0] == 'c' else data[2] for e in [e for e in exts if len(e[1]) >= s]: for i in range(len(trains)): if data[3](trains[:i]+[e]+trains[i:], data[1]) >= s: return # check for a closed/maximal item set n = len(trains) # report the (closed/maximal) item set if n < data[4] or n > data[5]: return if data[8] and n < len(data[8]) and supp <= data[8][n]: return # check against sizes and support border if data[7] == 'a': # if to collect the found pattern iset = tuple([i for i,t in trains]) data[6].append((iset,supp)) elif data[7] in '+-': # if map from sizes to support ranges if n not in data[6]: data[6][n] = [supp,supp] elif supp < data[6][n][0]: data[6][n][0] = supp elif supp > data[6][n][1]: data[6][n][1] = supp else: # if to collect a pattern spectrum s = (n,supp) # as a map from signatures to counts if s in rep: data[6][s] += 1 # update the else: data[6][s] = 1 # pattern spectrum #----------------------------------------------------------------------- def ovlreport (items, supp, ivals, elim, data): '''ovlreport (items, supp, ivals, elim, data) Check and report a (closed/maximal) frequent item set. item list of collected items supp support of the item set to report ivals list of intervals [start,end] resulting from intersecting the interval lists for the spike trains of the collected items elim list of trains of eliminated items, each train a pair (support, list of intervals) (needed to check whether the item set is closed or maximal) data recursion data, that is, a 9-tuple/list with: target type of frequent item set to find (in ['a','c','m']) width width of time window or influence region (width >= 0) supp minimum support (number of synchronous spiking events) suppfn support function (e.g. support or kpmsupp) zmin minimum size of an item set zmax maximum size of an item set result list in which the result is stored report reporting mode (pattern information/pattern spectrum) border support border (minimum per item set size)''' if data[0] != 'a': # if closed or maximal item sets s = supp if data[0] == 'c' else data[2] for e in reversed([e for t,e in elim if t >= s]): if sum([b-a for a,b in isect(ivals,e)]) /data[1] >= s: return # check for a closed/maximal item set n = len(items) # report the (closed/maximal) item set if n < data[4] or n > data[5]: return if data[8] and n < len(data[8]) and supp <= data[8][n]: return # check against sizes and support border if data[7] == 'a': # if to collect the found pattern data[6].append((items,supp)) elif data[7] in '+-': # if map from sizes to support ranges if n not in data[6]: data[6][n] = [supp,supp] elif supp < data[6][n][0]: data[6][n][0] = supp elif supp > data[6][n][1]: data[6][n][1] = supp else: # if to collect a pattern spectrum s = (n,supp) # as a map from signatures to counts if s in rep: data[6][s] += 1 # update the else: data[6][s] = 1 # pattern spectrum #----------------------------------------------------------------------- def pexcube (trains, exts, elim, supp, data): '''pexcube (trains, exts, elim, supp, data) Recursively traverse a hypercube spanned by perfect extensions. trains list of trains of collected items, each train a pair (item id, list of points/times) exts list of trains of extension items, each train a pair (item id, list of points/times) elim list of trains of eliminated items, each train a pair (item id, list of points/times) supp support of current item set (for closed/maximal check) data recursion data, that is, a 10-tuple/list with: target type of frequent item set to find (in ['a','c','m']) width width of time window or influence region (width >= 0) supp minimum support (number of synchronous spiking events) suppfn support function (e.g. support or kpmsupp) zmin minimum size of an item set zmax maximum size of an item set result list/dictionary in which the result is stored report reporting mode (pattern information/pattern spectrum) border support border (minimum per item set size) limit lower bound for support of trains union exts''' if supp <= data[9]: # if the support limit is reached report(trains+exts, elim, supp, data) return # report the full item set and abort k = len(elim); m = 0 # init. maximum extension support for i in range(len(exts)): # traverse the extension items z = trains+[exts[i]] # add train of extension item s = data[3](z, data[1]) # compute the support of the set if s < data[2]: continue# skip infrequent extensions and if s > m: m = s # find maximum extension support pexcube(z, exts[i+1:], elim, s, data) elim.append(exts[i]) # find frequent item sets recursively del elim[k:] # delete added eliminated items if data[0] == 'a' or m < (supp if data[0] == 'c' else data[2]): report(trains, elim, supp, data) # report current item set #----------------------------------------------------------------------- def rec_pex (trains, pexs, exts, elim, ivals, supp, data): '''rec_pex (trains, pexs, exts, elim, ivals, supp, data) Recursion of COntinuous time ClOsed Neuron Assembly Detection (version with perfect extension pruning). trains list of trains of collected items, each train a pair (item id, list of points/times) pexs list of trains of perfect extensions, each train a pair (item id, list of points/times) exts list of trains of extension items, each train a pair (item id, list of points/times) elim list of trains of eliminated items, each train a pair (item id, list of points/times) ivals list of intervals [start,end] for point/time filtering (intersection of intervals induced by the trains) supp support of current item set (for closed/maximal check) data recursion data, that is, a 9-tuple/list with: target type of frequent item set to find (in ['a','c','m']) width width of time window or influence region (width >= 0) supp minimum support (number of synchronous spiking events) suppfn support function (e.g. support or kpmsupp) zmin minimum size of an item set zmax maximum size of an item set result list/dictionary in which the result is stored report reporting mode (pattern information/pattern spectrum) border support border (minimum per item set size)''' exts.sort(key = lambda t: len(t[1])) p = len(pexs) # sort extension items by point count k = len(elim); m = 0 # and note old train array lengths for i in range(len(exts)): # traverse the extension items v = pfisect(ivals, exts[i][1], data[1]) z = [(n,ivfilt(t,v)) for n,t in trains+[exts[i]]] if [t for n,t in z if len(t) < data[2]]: continue # check for sufficiently long trains s = data[3](z, data[1]) # compute the support of the ext. set if s < data[2]: continue# skip infrequent and perfect extensions if s >= supp: pexs.append(exts[i]); continue if s > m: m = s # find maximum extension support x = [(n,ivfilt(t,v)) for n,t in exts[i+1:]] x = [t for t in x if len(t[1]) >= data[2]] e = [(n,ivfilt(t,v)) for n,t in elim] e = [t for t in e if len(t[1]) >= data[2]] rec_pex(z, pexs, x, e, v, s, data) elim.append(exts[i]) # find frequent item sets recursively z = trains +pexs # get the item set to report if pexs: # if there are perfect extensions, s = data[3](z, data[1]) # compute support of the full item set if s < supp: # if perfect extensions cause a problem data[9] = s # fix by traversing spanned hypercube x = [(n,ivfilt(t, ivals)) for n,t in pexs] x = [t for t in x if len(t[1]) >= data[2]] if m < s: del elim[k:] pexcube(trains, x, elim, supp, data) del elim[k:] # delete added eliminated items del pexs[p:] # and added perfect extensions return # and abort the function del elim[k:] # delete added eliminated items del pexs[p:] # and added perfect extensions if data[0] == 'c' or m < data[2]: report(z, elim, supp, data) # report current item set #----------------------------------------------------------------------- def rec_nopex (trains, exts, elim, ivals, supp, data): '''rec_nopex (trains, exts, elim, ivals, supp, data) Recursion of COntinuous time ClOsed Neuron Assembly Detection (version without perfect extension pruning). trains list of trains of collected items, each train a pair (item id, list of points/times) exts list of trains of extension items, each train a pair (item id, list of points/times) elim list of trains of eliminated items, each train a pair (item id, list of points/times) ivals list of intervals [start,end] for point/time filtering (intersection of intervals induced by the trains) supp support of current item set (for closed/maximal check) data recursion data, that is, a 9-tuple/list with: target type of frequent item set to find (in ['a','c','m']) width width of time window or influence region (width >= 0) supp minimum support (number of synchronous spiking events) suppfn support function (e.g. support or kpmsupp) zmin minimum size of an item set zmax maximum size of an item set result list/dictionary in which the result is stored report reporting mode (pattern information/pattern spectrum) border support border (minimum per item set size)''' exts.sort(key = lambda t: len(t[1])) k = len(elim); m = 0 # sort extension items by time count for i in range(len(exts)): # traverse the extension items v = pfisect(ivals, exts[i][1], data[1]) z = [(n,ivfilt(t,v)) for n,t in trains+[exts[i]]] if [t for n,t in z if len(t) < data[2]]: continue # check for sufficiently long trains s = data[3](z, data[1]) # compute the support of the ext. set if s < data[2]: continue# skip infrequent extensions if s > m: m = s # find maximum extension support x = [(n,ivfilt(t,v)) for n,t in exts[i+1:]] x = [t for t in x if len(t[1]) >= data[2]] if data[0] == 'a': e = elim else: # filter extensions/eliminated trains e = [(n,ivfilt(t,v)) for n,t in elim] e = [t for t in e if len(t[1]) >= data[2]] rec_nopex(z, x, e, v, s, data) elim.append(exts[i]) # find frequent item sets recursively del elim[k:] # delete added eliminated items if data[0] == 'a' or m < (supp if data[0] == 'c' else data[2]): report(trains, elim, supp, data) # report current item set #----------------------------------------------------------------------- def rec_seq (trains, exts, ivals, supp, data): '''rec_seq (trains, exts, ivals, supp, data) Recursion of COntinuous time ClOsed Neuron Assembly Detection for frequent item sequences. trains list of trains of collected items, each train a pair (item id, list of points/times) exts list of trains of extension items, each train a pair (item id, list of points/times) ivals list of intervals [start,end] for point/time filtering (intersection of intervals induced by the trains) supp support of current item set (for closed/maximal check) data recursion data, that is, a 9-tuple/list with: target type of frequent item set to find (in ['a','c','m']) width width of time window or influence region (width >= 0) supp minimum support (number of synchronous spiking events) suppfn support function (e.g. seqsupp) zmin minimum size of an item set zmax maximum size of an item set result list/dictionary in which the result is stored report reporting mode (pattern information/pattern spectrum) border support border (minimum per item set size)''' m = 0 # init. the maximum extension support for i in range(len(exts)): # traverse the extension items v = pfisect(ivals, exts[i][1], data[1]) z = [(n,ivfilt(t,v)) for n,t in trains+[exts[i]]] if [t for n,t in z if len(t) < data[2]]: continue # check for sufficiently long trains s = data[3](z, data[1]) # compute the support of the ext. set if s < data[2]: continue# skip infrequent extensions if s > m: m = s # find maximum extension support x = [(n,ivfilt(t,v)) for n,t in exts[:i]+exts[i+1:]] x = [t for t in x if len(t[1]) >= data[2]] rec_seq(z,x,v,s,data) # find frequent item sets recursively if data[0] == 'a' or m < (supp if data[0] == 'c' else data[2]): seqreport(trains, exts, supp, data) # report current item set #----------------------------------------------------------------------- def rec_ovl (items, exts, elim, ivals, supp, data): '''rec_ovl (items, exts, elim, ivals, supp, data) Recursion of COntinuous time ClOsed Neuron Assembly Detection for frequent item sequences, using overlap-based support. items list of collected items exts list of trains of extension items, each train a pair (item id, list of intervals) elim list of trains of eliminated items, each train a pair (support, list of intervals) ivals list of intervals [start,end] resulting from intersecting the interval lists for the spike trains of the collected items supp support of current item set (for closed/maximal check) data recursion data, that is, a 9-tuple/list with: target type of frequent item set to find (in ['a','c','m']) width width of time window or influence region (width >= 0) supp minimum support (number of synchronous spiking events) suppfn support function (e.g. seqsupp) zmin minimum size of an item set zmax maximum size of an item set result list/dictionary in which the result is stored report reporting mode (pattern information/pattern spectrum) border support border (minimum per item set size)''' exts.sort(key = lambda t: len(t[1])) k = len(elim); m = 0 # init. the maximum extension support for i in range(len(exts)): # traverse the extension items v = exts[i][1] # compute overlap based support s = sum([b-a for a,b in v]) /data[1] if s < data[2]: continue# skip infrequent extensions if s > m: m = s # find maximum extension support x = [(n,isect(t,v)) for n,t in exts[i+1:]] if data == 'a': e = elim else: # filter extensions/eliminated trains e = [(sum([b-a for a,b in e])/data[1],e) for e in [isect(t,v) for u,t in elim]] e = [a for a in e if a[0] >= data[2]] rec_ovl(items+[exts[i][0]], x, e, v, s, data) elim.append((s,v)) # find freq. item sets recursively del elim[k:] # delete added eliminated items if data[0] == 'a' or m < (supp if data[0] == 'c' else data[2]): ovlreport(items, supp, ivals, elim, data) # report item set #----------------------------------------------------------------------- def coconad (trains, target='c', width=0.003, supp=2.0, zmin=2, zmax=None, report='a', algo=support, mode='', border=None): '''coconad (trains, target='c', width=0.003, supp=2.0, zmin=2, zmax=None, report='a', algo=support, mode='', border=None) COntinuous time ClOsed Neuron Assembly Detection. trains list of trains of items to process, each train a pair (item id, sequence of points/times) target type of frequent item sets to find (default: c) a/s sets/all all frequent item sets c closed closed frequent item sets m maximal maximal frequent item sets as allseqs all frequent item sequences cs allseqs closed frequent item sequences ms allseqs maximal frequent item sequences supp minimum support of an item set (default: 2) width width of time window/influence region (default: 0.003) zmin minimum size of an item set (default: 2) zmax maximum size of an item set (default: no limit) report values to report with an item set (default: a) a absolute item set support (number of coincidences) # pattern spectrum as a map from signatures to counts = pattern spectrum as a list of triplets (size,supp,frq) + pattern spectrum as a map from sizes to support ranges - pattern spectrum as a list of triplets (size,supp,frq) algo algorithm indicator/support function (default: support) algorithm indicators (to be passed as a string) b/f/r binary support (from leading times) o influence map overlap support (graded) i interval list intersection support (graded) support functions (to be passed as function) support binary support (from leading times) kpmsupp k-partite matching support (graded) ovlsupp influence map overlap support (graded) mode operation mode indicators/flags (default: None) x do not use perfect extension pruning border support border (minimum support per item set size) returns the patterns as a list of pairs (tuple of item ids, support) or a pattern spectrum as a specified with the report argument''' if target in ['cls','clsd','closed']: target = 'c' elif target in ['max','maxi','maximal']: target = 'm' elif target in ['a','all','set','sets','frequent', 'allsets','freqsets']: target = 'ax' elif target in ['seq','seqs','sequences', 'allseq','allseqs']: target = 'as' elif target in ['clsseq','clsseqs','clsdseq','clsdseqs', 'closedseq','closedseqs']: target = 'cs' elif target in ['maxseq','maxseqs','maxiseq','maxiseqs', 'maximalseq','maximalseqs']: target = 'ms' if target == 's' or target == 'ss': target = 'a'+target[1:] if target not in ['a', 'c', 'm', 'as','cs','ms', 'ao','co','mo','aq','cq','mq']: target = 'c' seq = (len(target) > 1) # evaluate the target indicator if seq: target = target[:1] # and set the sequence flag width = float(width) # ensure a floating point width if supp < 0: supp = -supp # ensure a positive support if zmin <= 0: zmin = 1 # check and adapt the size range if not zmax or zmax <= 0: zmax = len(trains) trains = [(n,sorted(t)) for n,t in trains if len(t) >= supp] result = dict() if report in '#=-' else [] if algo in ['b','f','r']: suppfn = support elif algo in ['o']: suppfn = ovlsupp elif algo in ['i']: suppfn = 'isect' else: suppfn = algo pex = ('x' not in mode) # get perfect extension flag if seq: suppfn = seqsupp # build the recursion data if report in '#=+-' and suppfn not in [support,seqsupp]: report = '-' if report == '-' else '+' data = [target,width,supp,suppfn,zmin,zmax,result,report,border,0] if seq: # if to find frequent sequences rec_seq ([], trains, [[-oo,oo]], oo, data) elif suppfn == 'isect': # if influence overlap based support trains = [(i,train2ivals(t,width)) for i,t in trains] rec_ovl ([], trains, [], [[-oo,oo]], oo, data) elif pex and target != 'a': # if to use perfect extension pruning rec_pex ([], [], trains, [], [[-oo,oo]], oo, data) else: # if to do a plain search rec_nopex([], trains, [], [[-oo,oo]], oo, data) if report == '=': # if triplet format (sig. to counts) return [(s[0],s[1],result[s]) for s in result] if report == '-': # if triplet format (size to range) return [(s,result[s][0],result[s][1]) for s in result] return result # return the found patterns/spectrum #----------------------------------------------------------------------- # File Reading Function #----------------------------------------------------------------------- def getrec (f, recseps, fldseps, blanks, comment): '''Read a record from a file and split it into fields. f file to read from recseps record separators fldseps field separators blanks blank characters comment comment characters returns a list of the fields read or None at end of file''' c = f.read(1).decode() # read the next character while c != '' and c in comment: while c != '' and c not in recseps: c = f.read(1).decode() # skip comment until record separator c = f.read(1).decode() # consume the record separator if c == '': return None # if at end of file, abort rec = [] # initialize the record sep = fldseps +recseps # combine all separators while 1: # field read loop while c != '' and c in blanks: c = f.read(1).decode() # skip leading blanks fld = '' # initialize the field while c != '' and c not in sep: fld += c # append character to the field c = f.read(1).decode() # and read the next character fld = fld.strip(blanks) # remove leading/trailing blanks */ if len(fld) > 0: # if the field is not empty, rec.append(fld) # append the field to the record if c == '' or c in recseps: break c = f.read(1).decode() # consume field separator return rec # return the record #----------------------------------------------------------------------- # Main Program #----------------------------------------------------------------------- def error (msg): stderr.write(msg); exit() # print an error message and abort #----------------------------------------------------------------------- if __name__ == '__main__': desc = 'COntinuous time ClOsed Neuron Assembly Detection' version = 'version 1.21 (2016.03.22) ' \ + '(c) 2012-2016 Christian Borgelt' args = [] # list of program arguments/options opts = {'t': ['target', 'c' ], # target type 'o': ['seq', False], # flag for item sequences 'w': ['width', 0.003], # width of window/region 's': ['supp', 2 ], # minimum support (absolute) 'm': ['zmin', 2 ], # minimum number of items 'n': ['zmax', -1 ], # maximum number of items 'a': ['beg', -oo ], # minimum of point/time span 'z': ['end', oo ], # minimum of point/time span 'e': ['stype', 'b' ], # support type (b, g or o) 'A': ['algo', 'g' ], # support algorithm (g,t,b,i) 'd': ['eps', 0.0 ], # weight function parameter 'g': ['bound', 'n' ], # bounding function (n or w) 'x': ['pex', True ], # perfect extension flag 'q': ['pssep', ' ' ], # pattern spectrum separator 'h': ['ishdr', '' ], # record header for output 'k': ['isep', ' ' ], # item separator for output 'v': ['outfmt',' (%g)'],# output format for support 'l': ['train', False], # one train per line/record 'y': ['item', True ], # item in line, item first 'r': ['recseps', '\\n' ], # record separators 'f': ['fldseps', ' ,\\t' ], # field separators 'b': ['blanks', ' \\t\\r' ], # blank characters 'C': ['comment', '#' ], # comment characters 'P': ['pspfn', '' ] } # name of pattern spectrum file if len(argv) <= 1: # if no arguments are given opts = dict([(o,opts[o][1]) for o in opts]) print('usage: coconad.py [options] infile [outfile]') print(desc); print(version) print('-t# target type ' +'(default: %s)' % opts['t']) print(' (s: frequent, c: closed, m: maximal item sets)') print('-o find item sequences (respect item order) ' +'(default: sets)'); print('-w# width of time window/maximum distance ' +'(default: %g)' % opts['w']) print('-s# minimum support of an item set ' +'(default: %g)' % opts['s']) print('-m# minimum number of items per item set ' +'(default: %d)' % opts['m']) print('-n# maximum number of items per item set ' +'(default: no limit)') print('-a# start of point/time range ' +'(default: auto)') print('-z# end of point/time range ' +'(default: auto)') print('-e# support type (synchrony type) ' +'(default: %s)' % opts['e']) print(' (b: binary synchrony, g: graded synchrony,') print(' o: influence map overlap based synchrony)') print('-d# parameter for graded synchrony ' +'(default: %g)' % opts['d']) print(' (negative: binary, ' +'positive: epsilon-insensitive)') print('-A# algorithm for graded support ' +'(default: %c)' % opts['A']) print(' (g: greedy, t: backtracking, ' + 'b: branch and bound,') print(' i: interval list intersection') print('-g# bounding function for graded support ' +'(default: %c)' % opts['g']) print(' (n: edge count, w: edge weight, ' +'number: greedy)') print('-x do not prune with perfect extensions ' +'(default: prune)'); print('-P# write a pattern spectrum to a file') print('-q# column separator for pattern spectrum ' +'(default: "%s")' % opts['q']) print('-h# record header for output ' +'(default: "%s")' % opts['h']) print('-k# item separator for output ' +'(default: "%s")' % opts['k']) print('-v# output format for pattern support ' +'(default: "%s")' % opts['v']) print('-l one train per line (item + points) ' +'(default: id/time pairs)'); print('-y no item (with -l) or item after point ' +'(default: item first)'); print('-r# record separators ' +'(default: "%s")' % opts['r']) print('-f# field separators ' +'(default: "%s")' % opts['f']) print('-b# blank characters ' +'(default: "%s")' % opts['b']) print('-C# comment characters ' +'(default: "%s")' % opts['C']) print('infile file to read trains from ' +'[required]') print('outfile file to write found item sets to ' +'[optional]') exit() # print usage message and abort stderr.write('coconad.py') # print a startup message stderr.write(' - ' +desc +'\n' +version +'\n') # --- evaluate program options and arguments --- for a in argv[1:]: # traverse the program arguments if a[0] != '-': # check for a fixed argument if len(args) < 2: args.append(a); continue else: error('too many arguments\n') if a[1] not in opts: # check for an option error('unknown option: %s\n' % a[:2]) o = opts[a[1]] # get the corresponding option v = a[2:] # and get the option argument if isinstance(o[1],True.__class__): o[1] = not o[1] elif isinstance(o[1], (0).__class__): o[1] = int(v) elif isinstance(o[1], 0.0.__class__): o[1] = float(v) else: o[1] = v if len(args) < 1: error('missing arguments\n') opts = dict([opts[o] for o in opts]) target = opts['target'] # target type for CoCoNAD seq = opts['seq'] # flag for item sequences width = opts['width'] # time window width/maximum distance supp = opts['supp'] # minimum support for CoCoNAD zmin = opts['zmin'] # minimum size of item sets zmax = opts['zmax'] # maximum size of item sets beg = opts['beg'] # beginning of allowed span of spikes end = opts['end'] # end of allowed span of spikes pex = opts['pex'] # perfect extension flag pssep = opts['pssep'] # pattern spectrum separator ishdr = opts['ishdr'] # record header for output isep = opts['isep'] # item separator for output outfmt = opts['outfmt'] # output format for support tpr = opts['train'] # flag for one train per line item = opts['item'] # and flag for item first recseps = opts['recseps'] # record separators fldseps = opts['fldseps'] # field separators blanks = opts['blanks'] # blank characters comment = opts['comment'] # comment characters pspfn = opts['pspfn'] # pattern spectrum file name if zmin <= 0: error('invalid minimum size %d\n' % zmin) if supp <= 0: error('invalid minimum support %d\n' % supp) if width < 0: error('invalid window width %g\n' % width) if seq: # add sequence flag to target if len(target) <= 1: target += 's' elif 'seq' not in target: target += 'seqs' seq |= 'seq' in target or (len(target) == 2 and target[1] in 'soq') x = [recseps,fldseps,blanks,comment] if version_info[0] >= 3: # decode ASCII escape sequences x = [bytes(s, 'utf-8').decode('unicode_escape') for s in x] else: # decode ASCII escape sequences x = [s.decode('string_escape') for s in x] recseps,fldseps,blanks,comment = x stype = opts['stype'] # get the support type algo = opts['algo'] # and the support algorithm if seq: stype = 'b' # sequences only with binary support if stype == 'b': # if to use a binary support suppfn = seqsupp if seq else support elif stype == 'o': # if influence map overlap support suppfn = ovlsupp if algo != 'i' else 'isect' else: # if to use a graded support, eps = opts['eps'] # choose the weighting function if eps < 0: wgtfn = binwgt elif eps > 0: wgtfn = lambda x: epswgt(x, eps) else: wgtfn = idwgt if algo == 't': algo = 'btrack' elif algo == 'b': algo = 'bandb' else: algo = 'greedy' bound = opts['bound'] # choose the bounding heuristic if bound == 'n': bound = 'num' elif bound == 'w': bound = 'wgt' else: bound = float(bound) suppfn = lambda t,w: kpmsupp(t, w, wgtfn, algo, bound) # --- read the data set to analyze --- t = time() # start timer, print log message stderr.write('reading %s ... ' % args[0]) trains = dict(); k = 1 # initialize the train dictionary with open(args[0], 'rb') as inp: while 1: # record read loop rec = getrec(inp, recseps, fldseps, blanks, comment) if rec == None: break if not rec: continue# read the next record if tpr: # if one train per line/record if item: n = rec[0]; rec = rec[1:] else: n = k+1# get the item/neuron identifier trains[n] = [x for x in [float(x) for x in line] if x >= beg and x <= end] k += 1 # collect the trains else: # if item/ point pairs if item: n,s = rec else: s,n = rec s = float(s) # collect points per item/neuron if s < beg or s > end: continue if n in trains: trains[n].append(s) else: trains[n] = [s] trains = trains.items() # convert to pairs (item, list of times) points = sum([len(s) for n,s in trains]) stderr.write('[%d item(s), %d point(s)]' % (len(trains), points)) stderr.write(' done [%.2fs].\n' % (time()-t)) # --- run CoCoNAD algorithm --- t = time() # start timer, print log message if seq: stderr.write('mining frequent item sequences ... ') else: stderr.write('mining frequent item sets ... ') rep = 'a' if len(args) > 1 else '=' pats = coconad(trains, target, width, supp, zmin, zmax, rep, suppfn, '' if pex else 'x') if rep == '=': stderr.write('[%d signature(s)] ' % len(pats)) else: stderr.write('[%d pattern(s)] ' % len(pats)) stderr.write('done [%.2fs].\n' % (time()-t)) # --- save generated pattern spectrum --- if pspfn: # if file name for pattern spectrum if rep == '=': # if it has been collected already, psp = pats # simply get the pattern spectrum elif stype == 'b': # if a full pattern set was collected psp = dict() # collect pattern spectrum from it for p,c in pats: # traverse the patterns and s = (len(p),c) # construct the pattern signature if s in psp: psp[s] += 1 # then update the else: psp[s] = 1 # pattern spectrum else: # if a full pattern set was collected psp = dict() # collect pattern spectrum from it for p,c in pats: # traverse the patterns and n = len(p) # construct the pattern signature if n not in psp: psp[n] = [c,c] elif c < psp[n][0]: psp[n][0] = c elif c > psp[n][1]: psp[n][1] = c t = time() # start timer, print log message stderr.write('writing %s ... ' % pspfn) with open(pspfn, 'w') as out: if stype == 'b': # print (size, support, frequency) for s in sorted([(s[0],s[1],psp[s]) for s in psp]): out.write(('%d'+pssep+'%d'+pssep+'%.16g\n') % s) else: # print (size, min. supp., max. supp) for s in sorted([(s,psp[s][0],psp[s][1]) for s in psp]): out.write(('%d'+pssep+'%.16g'+pssep+'%.16g\n') % s) stderr.write('[%d signature(s)] ' % len(psp)) stderr.write('done [%.2fs].\n' % (time()-t)) # --- write output file --- if len(args) < 2: exit() # if no output file is given, abort t = time() # write the found patterns stderr.write('writing ' +args[1] +' ... ') with open(args[1], 'w') as out: for p,c in pats: # traverse and print the patterns for n in p[:-1]: out.write(n +isep) out.write(p[-1]) # write the items of the pattern out.write(outfmt % c) out.write('\n') # write the support information stderr.write('[%d pattern(s)]' % len(pats)) stderr.write(' done [%.2fs].\n' % (time()-t))