pynapple draft · MaxInterval and logISI
Burst detection, R beside Python
Each section puts a piece of the R reference on the left and the Python that answers to it on the right. Coloured bands mark the lines that do the same job. Hover or tap a band, or one of the labels above the code, to light up its partner on the other side.
The Python was written from the published description of each method, with the R read for constants and edge cases, so it is a counterpart and not a line-by-line translation. Where the two are structured differently, the notes under each pair say why the result is the same.
| R | Python | Section |
|---|---|---|
mi.find.bursts, phase 1 | _find_runs | Find runs |
mi.find.bursts, phase 2 | _merge_close_runs | Merge close bursts |
mi.find.bursts, phase 3 | _drop_small_runs, _runs_to_interval_set | Drop small bursts |
logisi.break.calc | _log_isi_histogram | Histogram |
get.peaks | _find_peaks | Peaks |
find.thresh | _log_isi_valley | Threshold |
logisi.pasq.method | detect_bursts_log_isi | Which threshold |
add.brs | _extend_to_related_spikes | Burst-related spikes |
logisi.find.burst | _find_runs and the two phases after it, reused | Where they disagree |
MaxInterval, phase 1: find runs
Both versions walk the spike train once with the same two-state machine. R keeps an in.burst flag and writes rows into a preallocated matrix; Python keeps start (None means not in a burst) and appends (first, last) index pairs. Indices go from 1-based to 0-based.
| Figure 3 | R, mi.par | Python | Default |
|---|---|---|---|
| max ISI at start of burst | beg.isi | max_isi_start | 0.17 s |
| max ISI in burst | end.isi | max_isi_burst | 0.3 s |
| min inter-burst interval | min.ibi | min_inter_burst_interval | 0.2 s |
| min burst duration | min.durn | min_burst_duration | 0.01 s |
| min number of spikes in burst | min.spikes | min_spikes_in_burst | 3 |
R reads these from a global list, mi.par. Python takes them as arguments; the defaults are Table 1 of Cotterill et al. 2016.
Rrepos/sjemea/R/maxinterval.Rlines 52–95
52 n = 2
53 in.burst = FALSE
54
55 while ( n <= nspikes) {
56
57 next.isi = spikes[n] - spikes[n-1]
58 if (in.burst) {
59 if (next.isi > end.isi) {
60 ## end of burst
61 end = n-1; in.burst = FALSE
62
63
64 ibi = spikes[beg] - last.end; last.end = spikes[end]
65 res = c(beg, end, ibi)
66 burst = burst + 1
67 if (burst > max.bursts) {
68 print("too many bursts!!!")
69 browser()
70 }
71 bursts[burst,] <- res
72 }
73 } else {
74 ## not yet in burst.
75 if (next.isi < beg.isi) {
76 ## Found the start of a new burst.
77 beg = n-1; in.burst = TRUE
78 }
79 }
80 n = n+1
81 }
82
83 ## At the end of the burst, check if we were in a burst when the
84 ## train finished.
85 if (in.burst) {
86 end = nspikes
87 ibi = spikes[beg] - last.end
88 res = c(beg, end, ibi)
89 burst = burst + 1
90 if (burst > max.bursts) {
91 print("too many bursts!!!")
92 browser()
93 }
94 bursts[burst,] <- res
95 }
Pythonnotes/drafts/pynapple_burst_detection.pylines 33–51
33def _find_runs(times, max_isi_start, max_isi_burst):
34 """
35 Phase 1: runs of spikes. A run starts at the first of two spikes closer than
36 `max_isi_start` and ends at the last spike before an interval longer than
37 `max_isi_burst`. Returns (first, last) spike indices, both inclusive.
38 """
39 runs = []
40 start = None
41 for i in range(1, len(times)):
42 isi = times[i] - times[i - 1]
43 if start is None:
44 if isi < max_isi_start:
45 start = i - 1
46 elif isi > max_isi_burst:
47 runs.append((start, i - 1))
48 start = None
49 if start is not None:
50 runs.append((start, len(times) - 1))
51 return runs
- The comparisons are the same: a burst starts on an interval strictly below the start threshold and ends on one strictly above the in-burst threshold.
- R stores each burst's gap to the previous one (
IBI) here for phase 2. Python does not; phase 2 reads the gap from the spike times. - The
max.burstspreallocation and itsbrowser()guard have no counterpart; a Python list grows.
MaxInterval, phase 2: merge close bursts
R marks every burst whose stored gap is below min.ibi, walks those in reverse to copy each end onto the burst before it, then deletes the marked rows. Python builds the merged list front to back: a run that begins too soon after the previous run's last spike extends that run.
Rrepos/sjemea/R/maxinterval.Rlines 112–131
112 ## Phase 2 -- merging of bursts. Here we see if any pair of bursts
113 ## have an IBI less than MIN.IBI; if so, we then merge the bursts.
114 ## We specifically need to check when say three bursts are merged
115 ## into one.
116
117
118 ibis = bursts[,"IBI"]
119 merge.bursts = which(ibis < min.ibi)
120
121 if (any(merge.bursts)) {
122 ## Merge bursts efficiently. Work backwards through the list, and
123 ## then delete the merged lines afterwards. This works when we
124 ## have say 3+ consecutive bursts that merge into one.
125
126 for (burst in rev(merge.bursts)) {
127 bursts[burst-1, "end"] = bursts[burst, "end"]
128 bursts[burst, "end"] = NA #not needed, but helpful.
129 }
130 bursts = bursts[-merge.bursts,,drop=FALSE] #delete the unwanted info.
131 }
Pythonnotes/drafts/pynapple_burst_detection.pylines 54–65
54def _merge_close_runs(times, runs, min_inter_burst_interval):
55 """
56 Phase 2: two consecutive runs separated by less than `min_inter_burst_interval`
57 (last spike of one to first spike of the next) become one run.
58 """
59 merged = []
60 for first, last in runs:
61 if merged and times[first] - times[merged[-1][1]] < min_inter_burst_interval:
62 merged[-1] = (merged[-1][0], last)
63 else:
64 merged.append((first, last))
65 return merged
- The gap is first spike of this run minus last spike of the previous original run in both versions, so a chain of three or more bursts merges the same way. R needs the reverse walk to get that; the forward build gets it for free.
MaxInterval, phase 3: drop small bursts
R computes a length and a duration column and deletes the rows that fail either test. Python keeps the runs that pass both.
Rrepos/sjemea/R/maxinterval.Rlines 139–174
139 ## Phase 3 -- remove small bursts: less than min duration (MIN.DURN), or
140 ## having too few spikes (less than MIN.SPIKES).
141 ## In this phase we have the possibility of deleting all spikes.
142
143 ## LEN = number of spikes in a burst.
144 ## DURN = duration of burst.
145 len = bursts[,"end"] - bursts[,"beg"] + 1
146 durn = spikes[bursts[,"end"]] - spikes[bursts[,"beg"]]
147 bursts = cbind(bursts, len, durn)
148
149 rejects = which ( (durn < min.durn) | ( len < min.spikes) )
150
151 if (any(rejects)) {
152 bursts = bursts[-rejects,,drop=FALSE]
153 }
154
155 if (nrow(bursts) == 0) {
156 ## All the bursts were removed during phase 3.
157 bursts = no.bursts
158 } else {
159 ## Compute mean ISIS
160 len = bursts[,"end"] - bursts[,"beg"] + 1
161 durn = spikes[bursts[,"end"]] - spikes[bursts[,"beg"]]
162 mean.isis = durn/(len-1)
163
164 ## Recompute IBI (only needed if phase 3 deleted some cells).
165 if (nrow(bursts)>1) {
166 ibi2 = c(NA, calc.ibi(spikes, bursts))
167 } else {
168 ibi2 = NA
169 }
170 bursts[,"IBI"] = ibi2
171
172 SI = rep(1, length(mean.isis ))
173 bursts = cbind(bursts, mean.isis, SI)
174 }
Pythonnotes/drafts/pynapple_burst_detection.pylines 68–85
68def _drop_small_runs(times, runs, min_burst_duration, min_spikes_in_burst):
69 """
70 Phase 3: a run shorter than `min_burst_duration` or with fewer than
71 `min_spikes_in_burst` spikes is not a burst.
72 """
73 return [
74 (first, last)
75 for first, last in runs
76 if times[last] - times[first] >= min_burst_duration
77 and last - first + 1 >= min_spikes_in_burst
78 ]
79
80
81def _runs_to_interval_set(times, runs):
82 starts = [times[first] for first, _ in runs]
83 ends = [times[last] for _, last in runs]
84 n_spikes = [last - first + 1 for first, last in runs]
85 return nap.IntervalSet(start=starts, end=ends, metadata={"n_spikes": n_spikes})
- R returns a matrix with
beg,end,IBI,len,durn,mean.isisandSI. Python returns anIntervalSetfrom first to last spike with one metadata column,n_spikes(R'slen). Duration and inter-burst intervals follow from the intervals themselves. - The same three Python functions serve logISI and the instantaneous-rate detector. In R,
logisi.find.burstis a second copy of this whole function.
logISI: the histogram
The threshold is read off a histogram of the interspike intervals on a logarithmic axis. Both versions work in milliseconds, leave out intervals under 1 ms, and normalise the counts.
CIBP is the intra-burst peak, Cmin the minimum before each later peak.Rrepos/burstanalysis/Burst_detection_methods/logisi_pasq_method.Rlines 112–126
112#Calculates cutoff for burst detection
113logisi.break.calc<-function(st, cutoff){
114 isi<-diff(st)*1000
115 max.isi<-ceiling(log10(max(isi)))
116 isi<-isi[isi>=1]
117 br<-logspace(0, max.isi, 10*max.isi)
118 h<-hist(isi, breaks=br, plot=FALSE)
119 h$density<-h$counts/sum(h$counts)
120 h$density<-lowess(h$density, f=0.05)$y
121 thr<-find.thresh(h, cutoff*1000)
122 if(!is.null(thr)){
123 thr<-thr/1000
124 }
125 thr
126}
Pythonnotes/drafts/pynapple_burst_detection.pylines 210–223
210def _log_isi_histogram(isi_ms):
211 """
212 Histogram of the interspike intervals (in milliseconds) on a logarithmic
213 axis: 10 bin edges per decade from 1 ms to the decade above the longest
214 interval, intervals under 1 ms left out, counts normalised to sum to one.
215 Returns (left edges in ms, normalised counts).
216 """
217 decades = max(1, int(np.ceil(np.log10(isi_ms.max()))))
218 edges = np.logspace(0, decades, 10 * decades)
219 isi_ms = isi_ms[isi_ms >= 1]
220 # Bin k holds edges[k] < isi <= edges[k + 1]; 1 ms itself goes to bin 0.
221 bins = np.clip(np.searchsorted(edges, isi_ms, side="left") - 1, 0, None)
222 counts = np.bincount(bins, minlength=len(edges) - 1)
223 return edges[:-1], counts / counts.sum()
- R's
histcloses bins on the right,(a, b]; NumPy'shistogramcloses them on the left. Python therefore bins withsearchsortedto get the same edges, so an interval sitting exactly on an edge lands in the same bin. - The smoothing line has no Python counterpart. With ten bins per decade the histogram has at most 60 bins for any train whose longest interval is under 1000 s. At
f = 0.05R'slowessthen uses a two-point neighbourhood and returns its input unchanged. Checked in R for 20 to 60 bins: zero difference. - The conversion to milliseconds (
diff(st) * 1000, R line 114) sits one function up in Python, in_log_isi_valley.
logISI: finding the peaks
This is where the two differ most line by line. R's get.peaks is a general peak finder with a minimum peak distance Pd, a height margin Th, a peak limit Np and a skip-ahead after each peak. It is only ever called with the defaults. With Pd = 2 and Th = 0 it reduces to one rule, which Python states directly.
Rrepos/burstanalysis/Burst_detection_methods/logisi_pasq_method.Rlines 38–75
38#Finds peaks in logISI histogram
39get.peaks<-function(h, Pd=2, Th=0, Np=NULL){
40 m<-0
41 L<-length(h$density)
42 j<-0
43 Np<-ifelse(is.null(Np), L, Np)
44 pks<-NULL
45 locs<-NULL
46 void.th<-0.7
47 while((j<L)&&(m<Np)){
48 j<-j+1
49 endL<-max(1,j-Pd)
50 if (m>0 && j<min(c(locs[m]+Pd, L-1))){
51 j<-min(c(locs[m]+Pd, L-1))
52 endL<-j-Pd
53 }
54 endR<-min(L, j+Pd)
55 temp<-h$density[endL:endR]
56 aa<-which(j==endL:endR)
57 temp[aa]<--Inf
58 if (Pd>1){
59 idx1<-max(1, aa-2)
60 idx2<-min(aa+2, length(temp))
61 idx3<-max(1, aa-1)
62 idx4<-min(aa+1, length(temp))
63 if (sum((h$density[j]>(temp[c(1:idx1, idx2:length(temp))]+Th))==FALSE)==0 && sum((h$density[j]>(temp[idx3:idx4]))==FALSE)==0 && j!=1 && j!=L){
64 m<-m+1
65 pks[m]<-h$density[j]
66 locs[m]<-j
67 } } else if (sum((h$density[j]>(temp+Th))==FALSE)==0 ) {
68 m<-m+1
69 pks[m]<-h$density[j]
70 locs[m]<-j
71 }
72
73 }
74 ret<-data.frame(pks=pks, locs=locs)
75}
Pythonnotes/drafts/pynapple_burst_detection.pylines 226–236
226def _find_peaks(density):
227 """
228 Bins strictly higher than every bin within two bins of them, excluding the
229 first and last bin.
230 """
231 peaks = []
232 for k in range(1, len(density) - 1):
233 neighbours = np.delete(density[max(0, k - 2) : k + 3], min(k, 2))
234 if np.all(density[k] > neighbours):
235 peaks.append(k)
236 return peaks
- The skip-ahead (R lines 50–53) changes nothing: a bin within two bins of a peak cannot itself be strictly higher than that peak.
- R tests the edge condition inside the comparison (
j != 1 && j != L); Python leaves the edge bins out of the loop. - The
Pd <= 1branch (lines 67–71) andNpare never reached with the defaults and are not carried over.
logISI: from peaks to threshold
Take the largest peak below the cutoff as the intra-burst peak. For each later peak, find the lowest bin between the two and score the separation with the void parameter. The threshold is the interval at the first minimum that scores 0.7 or more.
Rrepos/burstanalysis/Burst_detection_methods/logisi_pasq_method.Rlines 77–110
77#Function to find cutoff threshold.
78find.thresh<-function(h, ISITh=100){
79 void.th<-0.7
80 gp<-get.peaks(h)
81 num.peaks<-length(gp$pks)
82 pkx<-h$breaks[gp$locs]
83 intra.index<-which(pkx<ISITh)
84 if(length(intra.index)>=1){
85 max.intra<-max(gp$pks[intra.index])
86 max.idx<-which.max(gp$pks[intra.index])
87 } else {
88 return(-1000)
89 }
90 x1<-pkx[max.idx]
91 y1<-max.intra
92 locs1<-gp$locs[max.idx]
93 num.peaks.after.burst<-num.peaks-max.idx
94 if (num.peaks.after.burst==0){
95 return(NULL)
96 } else {
97 gp2<-gp[(max.idx+1):num.peaks,]
98 ymin<-sapply(gp2$locs, function(x) min(h$density[locs1:x]))
99 xmin<-sapply(gp2$locs, function(x) which.min(h$density[locs1:x]))+locs1-1
100 voidParameter<-1-(ymin/sqrt(y1*gp2$pks))
101 }
102 indxvoid<-suppressWarnings(min(which(voidParameter>=void.th)))
103 if(is.infinite(indxvoid)) {
104 flags<-c(1,0)
105 return(NULL)
106 } else {
107 ISImax<-h$breaks[xmin[indxvoid]]
108 return(ISImax)
109 }
110}
Pythonnotes/drafts/pynapple_burst_detection.pylines 239–264
239def _log_isi_valley(isi_s, max_cutoff):
240 """
241 The threshold of Pasquale et al. (2010): the ISI at the first valley that
242 separates the intra-burst peak of the log-ISI histogram from a later peak
243 with a void parameter of at least 0.7.
244
245 Returns (has_intraburst_peak, threshold in seconds or None).
246 """
247 isi_ms = isi_s * 1000
248 if len(isi_ms) == 0 or not np.any(isi_ms >= 1):
249 return False, None
250 edges, density = _log_isi_histogram(isi_ms)
251 peaks = _find_peaks(density)
252 intraburst = [k for k in peaks if edges[k] < max_cutoff * 1000]
253 if not intraburst:
254 return False, None
255 peak = max(intraburst, key=lambda k: density[k])
256 for later in peaks:
257 if later <= peak:
258 continue
259 between = density[peak : later + 1]
260 valley = peak + int(np.argmin(between))
261 void = 1 - between.min() / np.sqrt(density[peak] * density[later])
262 if void >= 0.7:
263 return True, float(edges[valley] / 1000)
264 return True, None
- R signals its three outcomes with three kinds of value:
-1000(no intra-burst peak),NULL(no clear valley) and a number. Python returns a pair,(has_intraburst_peak, threshold or None). - R computes the void parameter for every later peak and then takes the first that passes; Python stops at the first that passes. Same result.
- The reference uses
>=0.7 where the chapter writes>; Python follows the reference so the two agree.
logISI: which threshold to detect with
Four cases, in both versions. R builds a parameter list for each and calls logisi.find.burst; Python calls the three MaxInterval phases with the start and in-burst thresholds set to the same value.
Rrepos/burstanalysis/Burst_detection_methods/logisi_pasq_method.Rlines 4–36
4logisi.pasq.method<-function(spike.train, cutoff=0.1){
5 cutoff<-ifelse(is.null(cutoff), 0.1, cutoff)
6 if (length(spike.train)>3) {
7 isi.low <- logisi.break.calc(spike.train, cutoff) #Calculates threshold as isi.low
8 if (is.null(isi.low) || isi.low>=1 ){
9 logisi.par <- list(min.ibi=0, min.durn=0, min.spikes=3,
10 isi.low=cutoff) #If no value for isi.low found, or isi.low above 1 second, find bursts using threshold equal to cutoff (default 100ms)
11 result<-logisi.find.burst(spike.train, logisi.par)
12 } else if (isi.low<0) {
13 result<-NA
14 } else if (isi.low>cutoff & isi.low <1) {
15 logisi.par <- list(min.ibi=isi.low, min.durn=0, min.spikes=3,
16 isi.low=cutoff) #If isi.low >cutoff, find bursts using threshold equal to cutoff (default 100ms)
17 bursts<-logisi.find.burst(spike.train, logisi.par)
18 if (!is.na(bursts)[1]){
19 logisi.par2 <- list(min.ibi=0, min.durn=0, min.spikes=3,
20 isi.low=isi.low) #If bursts have been found, add burst related spikes using threshold of isi.low
21 brs<-logisi.find.burst(spike.train, logisi.par2)
22 result<-add.brs(bursts, brs, spike.train)
23 } else {
24 result<-bursts
25 }
26 } else {
27 logisi.par <- list(min.ibi=0, min.durn=0, min.spikes=3,
28 isi.low=isi.low) #If isi.low<cutoff, find bursts using a threshold equal to isi.low
29 result<-logisi.find.burst(spike.train, logisi.par)
30 }
31
32 } else {
33 result<-NA
34 }
35 result
36}
Pythonnotes/drafts/pynapple_burst_detection.pylines 424–443
424 has_peak, threshold = _log_isi_valley(_spike_isis(spikes), max_cutoff)
425 no_valley = threshold is None or threshold >= 1.0
426
427 def detect(times):
428 if not has_peak or (no_valley and not fallback):
429 return _runs_to_interval_set(times, [])
430 if no_valley:
431 runs = _find_runs(times, max_cutoff, max_cutoff)
432 runs = _drop_small_runs(times, runs, 0, min_spikes_in_burst)
433 elif threshold > max_cutoff:
434 cores = _find_runs(times, max_cutoff, max_cutoff)
435 cores = _merge_close_runs(times, cores, threshold)
436 cores = _drop_small_runs(times, cores, 0, min_spikes_in_burst)
437 runs = _extend_to_related_spikes(times, cores, threshold)
438 else:
439 runs = _find_runs(times, threshold, threshold)
440 runs = _drop_small_runs(times, runs, 0, min_spikes_in_burst)
441 return _runs_to_interval_set(times, runs)
442
443 return _detect_per_epoch(spikes, detect)
fallback=Falseis Python only and not in the published method: a train without a clear valley then gets no bursts instead of the fixed cutoff.- In the third case R passes
min.ibi = isi.low, so burst cores closer than the threshold merge before the three-spike minimum is applied. Python does the same with_merge_close_runs(times, cores, threshold). Two two-spike cores close together therefore survive as one burst in both. - R requires more than three spikes before doing anything (line 6). Python has no such guard; it needs only enough spikes to make a burst.
logISI: burst-related spikes
When the threshold lies above the cutoff, each core is widened to the spikes within the threshold of its ends. R detects a second set of bursts at the wider threshold (brs), searches it for the one overlapping each core with a nested loop and two early exits, then removes duplicates by comparing neighbours. Python looks up the run that contains the core and skips it if already taken.
Rrepos/burstanalysis/Burst_detection_methods/logisi_pasq_method.Rlines 131–179
131###Function to add burst related spikes to edges of bursts
132add.brs<-function(bursts, brs, spike.train){
133is.between<-function(x,a,b){ betw<-0
134 if(x>=a &x<=b)
135 {
136 betw<-1
137 }
138 betw
139}
140
141burst.adj<-data.frame(beg=rep(0, dim(bursts)[1]), end=rep(0, dim(bursts)[1]) )
142for (i in 1:dim(bursts)[1]) {
143 for (j in 1:dim(brs)[1]) {
144 if(is.between(bursts[i,1], brs[j,1], brs[j,2]) | is.between(bursts[i,2], brs[j,1], brs[j,2]))
145 {
146 burst.adj$beg[i]<-min(bursts[i,1], brs[j,1])
147 burst.adj$end[i]<-max(bursts[i,2], brs[j,2])
148 break
149 } else {
150 burst.adj$beg[i]<-bursts[i,1]
151 burst.adj$end[i]<-bursts[i,2]
152 }
153 if(brs[j,2]>bursts[i,2]) {
154 break
155 }
156 }
157}
158
159diff.begs<-diff(burst.adj[,"beg"])
160rep.bursts.begs<-which(diff.begs==0)
161if (any(rep.bursts.begs)) {
162 burst.adj<-burst.adj[-rep.bursts.begs,]
163}
164diff.ends<-diff(burst.adj[,"end"])
165rep.bursts.end<-which(diff.ends==0)+1
166if (any(rep.bursts.end)) {
167 burst.adj<-burst.adj[-rep.bursts.end,]
168}
169start.times<-spike.train[burst.adj$beg]
170end.times<- spike.train[burst.adj$end]
171durn<-end.times-start.times
172len<-burst.adj$end-burst.adj$beg+1
173mean.isis<-durn/(len-1)
174N.burst<-dim(burst.adj)[1]
175IBI<-c(NA, start.times[-1]-end.times[-N.burst])
176result<-cbind(beg=burst.adj$beg, end=burst.adj$end, IBI=IBI, len=len, durn=durn, mean.isis=mean.isis, SI=rep(1, N.burst))
177result
178
179}
Pythonnotes/drafts/pynapple_burst_detection.pylines 446–457
446def _extend_to_related_spikes(times, cores, threshold):
447 """
448 Extend each core to the run of spikes with consecutive intervals no longer
449 than `threshold` that contains it. Cores inside the same run become one burst.
450 """
451 runs = _find_runs(times, threshold, threshold)
452 extended = []
453 for first, last in cores:
454 container = next((r for r in runs if r[0] <= first and last <= r[1]), (first, last))
455 if container not in extended:
456 extended.append(container)
457 return extended
- Every core lies inside exactly one wider run, because its intervals are below the cutoff and the cutoff is below the threshold. So R's overlap test and Python's containment test pick the same run.
- R lines 169–177 rebuild the output columns (duration, mean ISI, IBI); in Python that is
_runs_to_interval_set, shared with every detector.
The one place they disagree
logisi.find.burst is a copy of the MaxInterval loop with two changes. One is a tolerance of 10-10 on the comparisons. The other is the loop bound: n < nspikes where MaxInterval has n <= nspikes. The last interval is never examined, and a burst still open at that point is closed on the last spike whatever that interval was. Python uses the one _find_runs for both methods and examines every interval.
Rrepos/burstanalysis/Burst_detection_methods/logisi_pasq_method.Rlines 226–270
226 eps<-10^(-10)
227 n = 2
228 in.burst = FALSE
229
230 while ( n < nspikes) {
231
232 next.isi = spikes[n] - spikes[n-1]
233 if (in.burst) {
234 if (next.isi - isi.low>eps) {
235 ## end of burst
236 end = n-1; in.burst = FALSE
237
238
239 ibi = spikes[beg] - last.end; last.end = spikes[end]
240 res = c(beg, end, ibi)
241 burst = burst + 1
242 if (burst > max.bursts) {
243 print("too many bursts!!!")
244 browser()
245 }
246 bursts[burst,] <- res
247 }
248 } else {
249 ## not yet in burst.
250 if (next.isi - isi.low <=eps) {
251 ## Found the start of a new burst.
252 beg = n-1; in.burst = TRUE
253 }
254 }
255 n = n+1
256 }
257
258 ## At the end of the burst, check if we were in a burst when the
259 ## train finished.
260 if (in.burst) {
261 end = nspikes
262 ibi = spikes[beg] - last.end
263 res = c(beg, end, ibi)
264 burst = burst + 1
265 if (burst > max.bursts) {
266 print("too many bursts!!!")
267 browser()
268 }
269 bursts[burst,] <- res
270 }
Pythonnotes/drafts/pynapple_burst_detection.pylines 33–51
33def _find_runs(times, max_isi_start, max_isi_burst):
34 """
35 Phase 1: runs of spikes. A run starts at the first of two spikes closer than
36 `max_isi_start` and ends at the last spike before an interval longer than
37 `max_isi_burst`. Returns (first, last) spike indices, both inclusive.
38 """
39 runs = []
40 start = None
41 for i in range(1, len(times)):
42 isi = times[i] - times[i - 1]
43 if start is None:
44 if isi < max_isi_start:
45 start = i - 1
46 elif isi > max_isi_burst:
47 runs.append((start, i - 1))
48 start = None
49 if start is not None:
50 runs.append((start, len(times) - 1))
51 return runs
- The comparison script reproduces this bound by wrapping
_find_runs, and counts a mismatch as explained only when the wrapped Python equals the R output exactly.
Agreement on 500 synthetic trains
Burst boundaries compared as spike indices: 125 trains each of bursty with Poisson background, slow bursts around the 100 ms cutoff, Poisson, and gamma renewal; 300 s each, 1,052,872 spikes in total. Run with python tests/_test_burst_r_reference.py.
| Exact | Explained by the loop bound | Unexplained | |
|---|---|---|---|
| MaxInterval bursts | 500 / 500 | – | 0 |
| logISI threshold | 500 / 500 | – | 0 |
| logISI bursts | 448 / 500 | 52 | 0 |
logISI cases covered: 192 trains with the threshold at or below the cutoff, 65 above it (burst-related spikes), 167 falling back to the cutoff, 76 with no intra-burst peak.