Fp fixes (#520)

* WIP daq wrong polarity

* photometry sync algorithm: compute cost function per frame shifts

* bump version number

ibllib/__init__.py

@@ -1,5 +1,5 @@
 """Library implementing the International Brain Laboratory data pipeline."""
-__version__ = "2.16.0"
+__version__ = "2.16.1"
 import warnings
 
 from iblutil.util import get_logger

ibllib/io/extractors/fibrephotometry.py

@@ -67,32 +67,52 @@ def sync_photometry_to_daq(vdaq, fs, df_photometry, chmap=DAQ_CHMAP, v_threshold
  :param v_threshold:
  :return:
  """
+ # here we take the flag that is the most common
  daq_frames, tag_daq_frames = read_daq_timestamps(vdaq=vdaq, v_threshold=v_threshold)
  nf = np.minimum(tag_daq_frames.size, df_photometry['Input0'].size)
- ipeak = np.argmax(np.correlate(tag_daq_frames[:nf].astype(np.int8), df_photometry['Input0'].values[:nf], mode='full'))
- # if the frame shift is negative, it means that the photometry frames are early
- frame_shift = ipeak - nf + 1
 
+ # we compute the framecounter for the DAQ, and match the bpod up state frame by frame for different shifts
+ # the shift that minimizes the mismatch is usually good
  df = np.median(np.diff(df_photometry['Timestamp']))
- fcn_fp2daq = scipy.interpolate.interp1d(df_photometry['Timestamp'][:nf], daq_frames[:nf] / fs)
- drift_ppm = (np.polyfit(daq_frames[:nf] / fs, df_photometry['Timestamp'][:nf], 1)[0] - 1) * 1e6
- _logger.info(f"drift photometry to DAQ PPM: {drift_ppm}")
-
+ fc = np.cumsum(np.round(np.diff(daq_frames) / fs / df).astype(np.int32)) - 1 # this is a daq frame counter
+ fc = fc[fc < (nf - 1)]
+ max_shift = 300
+ error = np.zeros(max_shift * 2 + 1)
+ shifts = np.arange(-max_shift, max_shift + 1)
+ for i, shift in enumerate(shifts):
+ rolled_fp = np.roll(df_photometry['Input0'].values[fc], shift)
+ error[i] = np.sum(np.abs(rolled_fp - tag_daq_frames[:fc.size]))
+ # a negative shift means that the DAQ is ahead of the photometry and that the DAQ misses frame at the beginning
+ frame_shift = shifts[np.argmax(-error)]
+ if np.sign(frame_shift) == -1:
+ ifp = fc[np.abs(frame_shift):]
+ elif np.sign(frame_shift) == 0:
+ ifp = fc
+ elif np.sign(frame_shift) == 1:
+ ifp = fc[:-np.abs(frame_shift)]
+ t_photometry = df_photometry['Timestamp'].values[ifp]
+ t_daq = daq_frames[:ifp.size] / fs
+ # import matplotlib.pyplot as plt
+ # plt.plot(shifts, -error)
+ fcn_fp2daq = scipy.interpolate.interp1d(t_photometry, t_daq, fill_value='extrapolate')
+ drift_ppm = (np.polyfit(t_daq, t_photometry, 1)[0] - 1) * 1e6
+ if drift_ppm > 120:
+ _logger.warning(f"drift photometry to DAQ PPM: {drift_ppm}")
+ else:
+ _logger.info(f"drift photometry to DAQ PPM: {drift_ppm}")
  # here is a bunch of safeguards
- assert np.all(np.abs(np.diff(daq_frames) - df * fs) < 1) # check that there are no missed frames on daq
  assert np.unique(np.diff(df_photometry['FrameCounter'])).size == 1 # checks that there are no missed frames on photo
- assert frame_shift == 0 # it's always the end frames that are missing
- assert drift_ppm < 20
-
+ assert np.abs(frame_shift) <= 5 # it's always the end frames that are missing
+ assert np.abs(drift_ppm) < 60
  ts_daq = fcn_fp2daq(df_photometry['Timestamp'].values) # those are the timestamps in daq time
-
  return ts_daq, fcn_fp2daq, drift_ppm
 
 
 def read_daq_voltage(daq_file, chmap=DAQ_CHMAP):
  channel_names = [c.name for c in load_raw_daq_tdms(daq_file)['Analog'].channels()]
  assert all([v in channel_names for v in chmap.values()]), "Missing channel"
  vdaq, fs = load_channels_tdms(daq_file, chmap=chmap, return_fs=True)
+ vdaq = {k: v - np.median(v) for k, v in vdaq.items()}
  return vdaq, fs
 
 
@@ -105,6 +125,10 @@ def read_daq_timestamps(vdaq, v_threshold=V_THRESHOLD):
  :return:
  """
  daq_frames = rises(vdaq['photometry'], step=v_threshold, analog=True)
+ if daq_frames.size == 0:
+ daq_frames = rises(-vdaq['photometry'], step=v_threshold, analog=True)
+ _logger.warning(f'No photometry pulses detected, attempting to reverse voltage and detect again,'
+ f'found {daq_frames.size} in reverse voltage. CHECK YOUR FP WIRING TO THE DAQ !!')
  tagged_frames = vdaq['bpod'][daq_frames] > v_threshold
  return daq_frames, tagged_frames
 

release_notes.md

@@ -1,3 +1,8 @@
+## Release 2.16.1
+### Release Notes 2.16.1 2022-09-28
+### bugfixes
+- photometry extraction: recover from corrupt DAQ signal and reversed polarity of voltage pulses
+
 ## Release 2.16
 ### Release Notes 2.16.0 2022-09-27
 ### features