spikeglx.Reader open=True

Author: k1o0

atlaselectrophysiology/extract_files.py

@@ -27,8 +27,7 @@ def rmsmap(fbin, spectra=True):
     :return: a dictionary with amplitudes in channeltime space, channelfrequency space, time
      and frequency scales
     """
-    if not isinstance(fbin, spikeglx.Reader):
-        sglx = spikeglx.Reader(fbin)
+    sglx = spikeglx.Reader(fbin, open=True)
     rms_win_length_samples = 2 ** np.ceil(np.log2(sglx.fs * RMS_WIN_LENGTH_SECS))
     # the window generator will generates window indices
     wingen = dsp.WindowGenerator(ns=sglx.ns, nswin=rms_win_length_samples, overlap=0)
@@ -58,6 +57,7 @@ def rmsmap(fbin, spectra=True):
         # print at least every 20 windows
         if (iw % min(20, max(int(np.floor(wingen.nwin / 75)), 1))) == 0:
             print_progress(iw, wingen.nwin)
+    sglx.close()
     return win
 
 

launch_phy/defined_metrics.py

@@ -84,30 +84,28 @@ def noise_cutoff(amps,quartile_length=.2, nbins=100, end_low = 2):
 
 def peak_to_peak_amp(ephys_file, samp_inds, nsamps):
 
-    #read raw ephys file
-    sr = spikeglx.Reader(ephys_file)
-
-    #take a subset (nsamps) of the spike samples
-    samples = np.random.choice(samp_inds,nsamps)
-    #initialize arrays
-    amps=np.zeros(len(samples))
-    wfs=np.zeros((384,len(samples)))
-    wfs_baseline=np.zeros((384,len(samples)))
-    cnt=0
-    for i in samples:
-        wf = sr.data[int(i)]
-        wf_baseline = wf[:-1]-np.median(wf[:-1]) #subtract median baseline
-        # plt.plot(wf_baseline)
-        wfs[:,cnt] = wf[:-1]
-        wfs_baseline[:,cnt] = wf_baseline
-        amps[cnt] = np.max(wf_baseline)-np.min(wf_baseline)
-        cnt+=1
-    amps = np.max(wfs_baseline,axis=0)-np.min(wfs_baseline,axis=0)
-    mean_amp = np.mean(amps)
-
-
+    # read raw ephys file
+    with spikeglx.Reader(ephys_file) as sr:
+        # take a subset (nsamps) of the spike samples
+        samples = np.random.choice(samp_inds,nsamps)
+        # initialize arrays
+        amps=np.zeros(len(samples))
+        wfs=np.zeros((384,len(samples)))
+        wfs_baseline=np.zeros((384,len(samples)))
+        cnt=0
+        for i in samples:
+            wf = sr.data[int(i)]
+            wf_baseline = wf[:-1]-np.median(wf[:-1])  # subtract median baseline
+            # plt.plot(wf_baseline)
+            wfs[:,cnt] = wf[:-1]
+            wfs_baseline[:,cnt] = wf_baseline
+            amps[cnt] = np.max(wf_baseline)-np.min(wf_baseline)
+            cnt+=1
+        amps = np.max(wfs_baseline,axis=0)-np.min(wfs_baseline,axis=0)
+        mean_amp = np.mean(amps)
     return mean_amp
 
+
 def max_acceptable_cont(FR, RP, rec_duration,acceptableCont, thresh ):
     trueContRate = np.linspace(0,FR,100)
     timeForViol = RP*2*FR*rec_duration
@@ -640,25 +638,25 @@ def ptp_over_noise(ephys_file, ts, ch, t=2.0, sr=30000, n_ch_probe=385, dtype='i
                                    np.min(wf[:, :, cur_ch], axis=1))
 
     # Compute MAD for `ch` in chunks.
-    s_reader = spikeglx.Reader(ephys_file)
-    file_m = s_reader.data  # the memmapped array
-    n_chunk_samples = 5e6  # number of samples per chunk
-    n_chunks = np.ceil(file_m.shape[0] / n_chunk_samples).astype('int')
-    # Get samples that make up each chunk. e.g. `chunk_sample[1] - chunk_sample[0]` are the
-    # samples that make up the first chunk.
-    chunk_sample = np.arange(0, file_m.shape[0], n_chunk_samples, dtype=int)
-    chunk_sample = np.append(chunk_sample, file_m.shape[0])
-    # Give time estimate for computing MAD.
-    t0 = time.perf_counter()
-    stats.median_absolute_deviation(file_m[chunk_sample[0]:chunk_sample[1], ch], axis=0)
-    dt = time.perf_counter() - t0
-    print('Performing MAD computation. Estimated time is {:.2f} mins.'
-          ' ({})'.format(dt * n_chunks / 60, time.ctime()))
-    # Compute MAD for each chunk, then take the median MAD of all chunks.
-    mad_chunks = np.zeros((n_chunks, ch.size), dtype=np.int16)
-    for chunk in range(n_chunks):
-        mad_chunks[chunk, :] = stats.median_absolute_deviation(
-            file_m[chunk_sample[chunk]:chunk_sample[chunk + 1], ch], axis=0, scale=1)
+    with spikeglx.Reader(ephys_file) as s_reader:
+        file_m = s_reader.data  # the memmapped array
+        n_chunk_samples = 5e6  # number of samples per chunk
+        n_chunks = np.ceil(file_m.shape[0] / n_chunk_samples).astype('int')
+        # Get samples that make up each chunk. e.g. `chunk_sample[1] - chunk_sample[0]` are the
+        # samples that make up the first chunk.
+        chunk_sample = np.arange(0, file_m.shape[0], n_chunk_samples, dtype=int)
+        chunk_sample = np.append(chunk_sample, file_m.shape[0])
+        # Give time estimate for computing MAD.
+        t0 = time.perf_counter()
+        stats.median_absolute_deviation(file_m[chunk_sample[0]:chunk_sample[1], ch], axis=0)
+        dt = time.perf_counter() - t0
+        print('Performing MAD computation. Estimated time is {:.2f} mins.'
+              ' ({})'.format(dt * n_chunks / 60, time.ctime()))
+        # Compute MAD for each chunk, then take the median MAD of all chunks.
+        mad_chunks = np.zeros((n_chunks, ch.size), dtype=np.int16)
+        for chunk in range(n_chunks):
+            mad_chunks[chunk, :] = stats.median_absolute_deviation(
+                file_m[chunk_sample[chunk]:chunk_sample[chunk + 1], ch], axis=0, scale=1)
     print('Done. ({})'.format(time.ctime()))
 
     # Return `mean_ptp` over `mad`

viewspikes/data.py

@@ -53,6 +53,15 @@ def get_spikes(dsets, one):
 
 
 def stream(pid, t0, one=None, cache=True, dsets=None):
+    """
+    NB: returned Reader object must be closed after use
+    :param pid: Probe UUID
+    :param t0:
+    :param one: An instance of ONE
+    :param cache:
+    :param dsets:
+    :return:
+    """
     tlen = 1
     assert one
     if cache:
@@ -62,7 +71,8 @@ def stream(pid, t0, one=None, cache=True, dsets=None):
         dsets = one.alyx.rest('datasets', 'list', probe_insertion=pid)
     if cache and samples_folder.joinpath(sample_file_name).exists():
         print(f'loading {sample_file_name} from cache')
-        sr = spikeglx.Reader(samples_folder.joinpath(sample_file_name).with_suffix('.bin'))
+        sr = spikeglx.Reader(samples_folder.joinpath(sample_file_name).with_suffix('.bin'),
+                             open=True)
         return sr, dsets
 
     dset_ch = next(dset for dset in dsets if dset['dataset_type'] == "ephysData.raw.ch" and '.ap.' in dset['name'])
@@ -85,6 +95,7 @@ def stream(pid, t0, one=None, cache=True, dsets=None):
         out_meta = samples_folder.joinpath(sample_file_name)
         shutil.copy(sr.file_meta_data, out_meta)
         with open(out_meta.with_suffix('.bin'), 'wb') as fp:
+            sr.open()
             sr._raw[:].tofile(fp)
 
     return sr, dsets

viewspikes/main.py

@@ -56,7 +56,7 @@
     file_sample = pid
     pid, t0 = file_sample.stem.split('_')
     t0 = float(t0)
-    sr = spikeglx.Reader(file_sample)
+    sr = spikeglx.Reader(file_sample, open=True)
     dsets = one.alyx.rest('datasets', 'list', probe_insertion=pid)
 else:
     sr, dsets = stream(pid, t0, one=one, samples_folder=folder_samples)
@@ -90,3 +90,4 @@
     overlay_spikes(eqc_butt, spikes, clusters, channels)
     overlay_spikes(eqc_dest, spikes, clusters, channels)
     overlay_spikes(eqc_ks2, spikes, clusters, channels)
+sr.close()