Fixes the second bug in the spike holes issue #594

CUDA/mexMPnu8.cu

@@ -324,9 +324,11 @@ __global__ void  bestFilter(const double *Params, const float *data,
               if  (Cf > Cnextbest + 1e-6)
                     Cnextbest = Cf;
       }
-      err[tid0] 	= Cbest;
-      eloss[tid0] 	= Cbest - Cnextbest;
-      ftype[tid0] 	= ibest;
+      if (tid0>=nt0){
+          err[tid0] 	= Cbest;
+          eloss[tid0] 	= Cbest - Cnextbest;
+          ftype[tid0] 	= ibest;
+      }
 
       tid0 += blockDim.x * gridDim.x;
   }
@@ -351,7 +353,7 @@ __global__ void  bestFilterUpdate(const double *Params, const float *data,
 
   if (ind<counter[0]){
       t = st[ind]-nt0 + tid;
-      if (t>=0 && t<NT){
+      if (t>=nt0 && t<NT){
           Cbest = 0.0f;
           for (i=0; i<Nfilt;i++){
               a = 1+ lam;

configFiles/configFile384.m

@@ -41,8 +41,11 @@
 ops.GPU                 = 1; % has to be 1, no CPU version yet, sorry
 % ops.Nfilt               = 1024; % max number of clusters
 ops.nfilt_factor        = 4; % max number of clusters per good channel (even temporary ones)
-ops.ntbuff              = 64;    % samples of symmetrical buffer for whitening and spike detection
-ops.NT                  = 64*1024+ ops.ntbuff; % must be multiple of 32 + ntbuff. This is the batch size (try decreasing if out of memory). 
+
+% changing NT and ntbuff may introduce the spike holes problem issue #594. 
+ops.ntbuff              = 64;    % MUST BE MULTIPLE OF 32. Samples of symmetrical buffer for whitening and spike detection. 
+ops.NT                  = 64*1024 - ops.ntbuff; % MUST BE MULTIPLE OF 1024 - ntbuff. This is the batch size (try decreasing if out of memory). 
+
 ops.whiteningRange      = 32; % number of channels to use for whitening each channel
 ops.nSkipCov            = 25; % compute whitening matrix from every N-th batch
 ops.scaleproc           = 200;   % int16 scaling of whitened data

mainLoop/trackAndSort.m

@@ -56,7 +56,11 @@
 % spike threshold for finding missed spikes in residuals
 ops.spkTh = -6; % why am I overwriting this here?
 
-batchstart = 0:NT:NT*nBatches;
+% This was updated in 2.5.2 (fixes for spike holes issue #594)
+ntpad = ops.ntbuff; %16 * ceil(2*nt0 / 32);
+batchstart = ntpad:NT:NT*nBatches;
+
+%batchstart = 0:NT:NT*nBatches;
 
 % find the closest NchanNear channels, and the masks for those channels
 [iC, mask, C2C] = getClosestChannels(rez, sigmaMask, NchanNear);
@@ -120,12 +124,28 @@
 for ibatch = 1:niter    
     k = iorder(ibatch); % k is the index of the batch in absolute terms
 
+    % new version (fixing spikes hole issue #594) -------------------------
     % loading a single batch (same as everywhere)
-    offset = 2 * ops.Nchan*batchstart(k);
+    offset = 2 * ops.Nchan*(batchstart(k) - ntpad);
     fseek(fid, offset, 'bof');
-    % dat = fread(fid, [ops.Nchan NT + ops.ntbuff], '*int16');
+    %dat = fread(fid, [ops.Nchan NT + ops.ntbuff], '*int16'); 
     % spike holes bug issue #594, location 1/2
-    dat = fread(fid, [ops.Nchan NT], '*int16');
+    dat = fread(fid, [ops.Nchan NT+2*ntpad], '*int16');
+    if size(dat,2)<NT+2*ntpad
+        NT_current = size(dat,2);
+        ndiff = NT+2*ntpad - NT_current;
+        dat(:, NT_current+1:NT+2*ntpad) = repmat(dat(:, NT_current), 1, ndiff);
+    end
+
+    % ---------------------------------------------------------------------
+
+     % old version
+%     % loading a single batch (same as everywhere)
+%     offset = 2 * ops.Nchan*batchstart(k);
+%     fseek(fid, offset, 'bof');
+%     % dat = fread(fid, [ops.Nchan NT + ops.ntbuff], '*int16');
+%     % spike holes bug issue #594, location 1/2
+%     dat = fread(fid, [ops.Nchan NT], '*int16');
 
     dat = dat';
     dataRAW = single(gpuArray(dat))/ ops.scaleproc;
@@ -156,6 +176,18 @@
         mexMPnu8(Params, dataRAW, single(U), single(W), single(mu), iC-1, iW-1, UtU, iList-1, ...
         wPCA);
 
+
+    % exclude spike times in the padding (to fix issue #594)
+    ix = (st0>ntpad) & (st0<=NT+ntpad);
+    %ix = [1:numel(st0)] > nmax;
+
+    st0 = st0(ix);
+    id0 = id0(ix);
+    x0 = x0(ix);
+    featW = featW(:, ix);
+    featPC = featPC(:,:,ix);
+    vexp = vexp(ix); 
+
     % \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
     % \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
 
@@ -189,11 +221,14 @@
     % we carefully assign the correct absolute times to spikes found in this batch
     % toff = nt0min + t0 + NT*(k-1);
     % spike holes bug issue #594, location 2/2    
-    ioffset         = ops.ntbuff;
-    if k==1
-        ioffset         = 0; % the first batch is special (no pre-buffer)
-    end
-    toff = nt0min + t0 -ioffset + (NT-ops.ntbuff)*(k-1);
+%     ioffset         = ops.ntbuff;
+%     if k==1
+%         ioffset         = 0; % the first batch is special (no pre-buffer)
+%     end
+%     toff = nt0min + t0 -ioffset + (NT-ops.ntbuff)*(k-1);
+
+    % new version in patch 2.5.2
+    toff = t0 + nt0min -ntpad + batchstart(k); 
 
     st = toff + double(st0);