Merge branch 'develop' into bb_passive_fix

Author: mayofaulkner

brainbox/io/one.py

@@ -70,11 +70,11 @@ def load_channel_locations(eid, one=None, probe=None, aligned=False):
             counts = [0]
         else:
             tracing = [(insertions.get('json', {'temp': 0}).get('extended_qc', {'temp': 0}).
-                       get('tracing_exists', False))]
+                        get('tracing_exists', False))]
             resolved = [(insertions.get('json', {'temp': 0}).get('extended_qc', {'temp': 0}).
-                        get('alignment_resolved', False))]
+                         get('alignment_resolved', False))]
             counts = [(insertions.get('json', {'temp': 0}).get('extended_qc', {'temp': 0}).
-                      get('alignment_count', 0))]
+                       get('alignment_count', 0))]
         probe_id = [insertions['id']]
     # No specific probe specified, load any that is available
     # Need to catch for the case where we have two of the same probe insertions
@@ -420,7 +420,7 @@ def load_wheel_reaction_times(eid, one=None):
 
 
 def load_trials_df(eid, one=None, maxlen=None, t_before=0., t_after=0., ret_wheel=False,
-                   ret_abswheel=False, wheel_binsize=0.02, addtl_types=()):
+                   ret_abswheel=False, ext_DLC=False, wheel_binsize=0.02, addtl_types=[]):
     """
     TODO Test this with new ONE
     Generate a pandas dataframe of per-trial timing information about a given session.
@@ -451,6 +451,8 @@ def load_trials_df(eid, one=None, maxlen=None, t_before=0., t_after=0., ret_whee
         Whether to return the time-resampled wheel velocity trace, by default False
     ret_abswheel : bool, optional
         Whether to return the time-resampled absolute wheel velocity trace, by default False
+    ext_DLC : bool, optional
+        Whether to extract DLC data, by default False
     wheel_binsize : float, optional
         Time bins to resample wheel velocity to, by default 0.02
     addtl_types : list, optional

brainbox/modeling/linear.py

@@ -59,16 +59,38 @@ def __init__(self, design_matrix, spk_times, spk_clu,
     def _fit(self, dm, binned, cells=None):
         """
         Fitting primitive that brainbox.NeuralModel.fit method will call
+
+        Parameters
+        ----------
+        dm : np.ndarray
+            Design matrix to use for fitting
+        binned : np.ndarray
+            Array of binned spike times. Must share first dimension with dm
+        cells : iterable with .shape attribute, optional
+            List of cells which are being fit. Use to generate index for output
+            coefficients and intercepts, must share shape with second dimension
+            of binned. When None will default to a list of all cells in the model object,
+            by default None
+
+        Returns
+        -------
+        coefs, pd.Series
+            Series containing fit coefficients for cells
+        intercepts, pd.Series
+            Series containing intercepts for fits.
         """
         if cells is None:
             cells = self.clu_ids.flatten()
+        if cells.shape[0] != binned.shape[1]:
+            raise ValueError('Length of cells does not match shape of binned')
+
         coefs = pd.Series(index=cells, name='coefficients', dtype=object)
         intercepts = pd.Series(index=cells, name='intercepts')
 
         lm = self.estimator.fit(dm, binned)
         weight, intercept = lm.coef_, lm.intercept_
         for cell in cells:
-            cell_idx = np.argwhere(self.clu_ids == cell)[0, 0]
+            cell_idx = np.argwhere(cells == cell)[0, 0]
             coefs.at[cell] = weight[cell_idx, :]
             intercepts.at[cell] = intercept[cell_idx]
         return coefs, intercepts
@@ -84,7 +106,6 @@ def score(self):
         """
         if not hasattr(self, 'coefs'):
             raise AttributeError('Model has not been fit yet.')
-
         testmask = np.isin(self.design.trlabels, self.testinds).flatten()
         dm, binned = self.design[testmask, :], self.binnedspikes[testmask]
 

brainbox/modeling/neural_model.py

@@ -103,7 +103,7 @@ def __init__(self, design_matrix, spk_times, spk_clu,
         self.design = design_matrix
         self.spikes = spks
         self.clu = clu
-        self.clu_ids = np.argwhere(np.sum(trialspiking, axis=0) > mintrials)
+        self.clu_ids = np.argwhere(np.sum(trialspiking, axis=0) > mintrials).flatten()
         self.traininds = traininds
         self.testinds = testinds
         self.stepwise = stepwise

brainbox/modeling/poisson.py

@@ -38,16 +38,19 @@ def _fit(self, dm, binned, cells=None, noncovwarn=False):
         alpha : float
             Regularization strength, applied as multiplicative constant on ridge regularization.
         cells : list
-            List of cells which should be fit. If None is passed, will default to fitting all cells
-            in clu_ids
+            List of cells labels for columns in binned. Will default to all cells in model if None
+            is passed. Must be of the same length as columns in binned. By default None.
         """
         if cells is None:
             cells = self.clu_ids.flatten()
+        if cells.shape[0] != binned.shape[1]:
+            raise ValueError('Length of cells does not match shape of binned')
+
         coefs = pd.Series(index=cells, name='coefficients', dtype=object)
         intercepts = pd.Series(index=cells, name='intercepts')
         nonconverged = []
         for cell in tqdm(cells, 'Fitting units:', leave=False):
-            cell_idx = np.argwhere(self.clu_ids == cell)[0, 0]
+            cell_idx = np.argwhere(cells == cell)[0, 0]
             cellbinned = binned[:, cell_idx]
             with catch_warnings(record=True) as w:
                 fitobj = PoissonRegressor(alpha=self.alpha,

brainbox/modeling/utils.py

@@ -101,7 +101,7 @@ def fit(self, progress=False):
                 new_feature_idx, nf_score = self._get_best_new_feature(current_mask, cells)
                 for cell in cells:
                     maskdf.at[cell, self.features[new_feature_idx.loc[cell]]] = True
-                    seqdf.loc[cell, i] = self.features[new_feature_idx]
+                    seqdf.loc[cell, i] = self.features[new_feature_idx.loc[cell]]
                     scoredf.loc[cell, i] = nf_score.loc[cell]
         self.support_ = maskdf
         self.sequences_ = seqdf
@@ -110,7 +110,8 @@ def fit(self, progress=False):
     def _get_best_new_feature(self, mask, cells):
         mask = np.array(mask)
         candidate_features = np.flatnonzero(~mask)
-        my = self.model.binnedspikes[self.train]
+        cell_idxs = np.argwhere(np.isin(self.model.clu_ids, cells)).flatten()
+        my = self.model.binnedspikes[np.ix_(self.train, cell_idxs)]
         scores = pd.DataFrame(index=cells, columns=candidate_features, dtype=float)
         for feature_idx in candidate_features:
             candidate_mask = mask.copy()

examples/one/tutorial_script.py

@@ -60,20 +60,37 @@
           }
 
 
+def data_for_keys(keys, data):
+    """Check keys exist in 'data' dict and contain values other than None"""
+    return data is not None and all(k in data and data.get(k, None) is not None for k in keys)
+
+
 def get_ibl_sync_map(ef, version):
     """
     Gets default channel map for the version/binary file type combination
     :param ef: ibllib.io.spikeglx.glob_ephys_file dictionary with field 'ap' or 'nidq'
     :return: channel map dictionary
     """
+    # Determine default channel map
     if version == '3A':
         default_chmap = CHMAPS['3A']['ap']
     elif version == '3B':
         if ef.get('nidq', None):
             default_chmap = CHMAPS['3B']['nidq']
         elif ef.get('ap', None):
             default_chmap = CHMAPS['3B']['ap']
-    return spikeglx.get_sync_map(ef['path']) or default_chmap
+    # Try to load channel map from file
+    chmap = spikeglx.get_sync_map(ef['path'])
+    # If chmap provided but not with all keys, fill up with default values
+    if not chmap:
+        return default_chmap
+    else:
+        if data_for_keys(default_chmap.keys(), chmap):
+            return chmap
+        else:
+            _logger.warning("Keys missing from provided channel map, "
+                            "setting missing keys from default channel map")
+            return {**default_chmap, **chmap}
 
 
 def _sync_to_alf(raw_ephys_apfile, output_path=None, save=False, parts=''):
@@ -242,7 +259,8 @@ def _assign_events_audio(audio_t, audio_polarities, return_indices=False):
     # make sure that there are no 2 consecutive fall or consecutive rise events
     assert(np.all(np.abs(np.diff(audio_polarities)) == 2))
     # take only even time differences: ie. from rising to falling fronts
-    dt = np.diff(audio_t)[::2]
+    i0 = 0 if audio_polarities[0] == 1 else 1
+    dt = np.diff(audio_t)[i0::2]
     # detect ready tone by length below 110 ms
     i_ready_tone_in = np.r_[np.where(dt <= 0.11)[0] * 2]
     t_ready_tone_in = audio_t[i_ready_tone_in]

ibllib/io/extractors/ephys_fpga.py

@@ -125,7 +125,19 @@ def _get_spacer_times(spacer_template, jitter, ttl_signal, t_quiet):
     # adjust indices for
     # - `np.where` call above
     # - length of spacer_model
-    idxs_spacer_middle += 2 - int((np.floor(len(spacer_model) / 2)))
+    spacer_around = int((np.floor(len(spacer_model) / 2)))
+    idxs_spacer_middle += 2 - spacer_around
+
+    # for each spacer make sure the times are monotonically increasing before
+    # and monotonically decreasing afterwards
+    is_valid = np.zeros((idxs_spacer_middle.size), dtype=bool)
+    for i, t in enumerate(idxs_spacer_middle):
+        before = all(np.diff(dttl[t - spacer_around:t]) > 0)
+        after = all(np.diff(dttl[t + 1:t + 1 + spacer_around]) < 0)
+        is_valid[i] = np.bitwise_and(before, after)
+
+    idxs_spacer_middle = idxs_spacer_middle[is_valid]
+
     # pull out spacer times (middle)
     ts_spacer_middle = ttl_signal[idxs_spacer_middle]
     # put beginning/end of spacer times into an array

ibllib/io/extractors/ephys_passive.py

@@ -93,7 +93,7 @@ class SpikeSorting(tasks.Task):
         "Documents/PYTHON/iblscripts/deploy/serverpc/kilosort2/run_pykilosort.sh"
     )
     SPIKE_SORTER_NAME = 'pykilosort'
-    PYKILOSORT_REPO = '~/Documents/PYTHON/SPIKE_SORTING/pykilosort'
+    PYKILOSORT_REPO = Path.home().joinpath('Documents/PYTHON/SPIKE_SORTING/pykilosort')
 
     @staticmethod
     def _sample2v(ap_file):
@@ -113,7 +113,7 @@ def _fetch_pykilosort_version(repo_path):
                         version = line.split('=')[-1].strip().replace('"', '').replace("'", '')
         except Exception:
             pass
-        return version
+        return f"pykilosort_{version}"
 
     @staticmethod
     def _fetch_ks2_commit_hash(repo_path):
@@ -137,7 +137,7 @@ def _run_pykilosort(self, ap_file):
         session_path/spike_sorters/{self.SPIKE_SORTER_NAME}/probeXX folder
         :return: path of the folder containing ks2 spike sorting output
         """
-
+        self.version = self._fetch_pykilosort_version(self.PYKILOSORT_REPO)
         label = ap_file.parts[-2]  # this is usually the probe name
         if ap_file.parent.joinpath(f"spike_sorting_{self.SPIKE_SORTER_NAME}.log").exists():
             _logger.info(f"Already ran: spike_sorting_{self.SPIKE_SORTER_NAME}.log"
@@ -194,7 +194,6 @@ def _run_pykilosort(self, ap_file):
 
         shutil.copytree(temp_dir.joinpath('output'), sorter_dir, dirs_exist_ok=True)
         shutil.rmtree(temp_dir, ignore_errors=True)
-        self.version = self._fetch_ks2_commit_hash(self.PYKILOSORT_REPO)
         return sorter_dir
 
     def _run(self, overwrite=False):
@@ -225,7 +224,7 @@ def _run(self, overwrite=False):
                 logfile = ks2_dir.joinpath(f"spike_sorting_{self.SPIKE_SORTER_NAME}.log")
                 if logfile.exists():
                     shutil.copyfile(logfile, probe_out_path.joinpath(
-                        f"ibl_log.info_{self.SPIKE_SORTER_NAME}.log"))
+                        f"_ibl_log.info_{self.SPIKE_SORTER_NAME}.log"))
                 out, _ = spikes.sync_spike_sorting(ap_file=ap_file, out_path=probe_out_path)
                 out_files.extend(out)
                 # convert ks2_output into tar file and also register

ibllib/pipes/ephys_preprocessing.py

@@ -274,6 +274,9 @@ def get_active_wheel_period(wheel, duration_range=(3., 20.), display=False):
         edges = np.c_[on, off]
         indices, _ = np.where(np.logical_and(
             np.diff(edges) > duration_range[0], np.diff(edges) < duration_range[1]))
+        if len(indices) == 0:
+            _log.warning('No period of wheel movement found for motion alignment.')
+            return None
         # Pick movement somewhere in the middle
         i = indices[int(indices.size / 2)]
         if display:
@@ -304,7 +307,7 @@ def ensure_required_data(self):
             # Assert 3A probe model; if so download all probe data
             det = self.one.get_details(self.eid, full=True)
             probe_model = next(x['model'] for x in det['probe_insertion'])
-            assert probe_model == '3A', 'raw ephys data not missing'
+            assert probe_model == '3A', 'raw ephys data missing'
             collections += ('raw_ephys_data/probe00', 'raw_ephys_data/probe01')
             assert_unique = False
         for dstype in dtypes: