GH-1404: Handle monitor delay errors

Previously we used consective photodiode edges of the same type
to determine whether an event is a valid transition. This led
sometimes to missing out on the final valid event, and resulted
in an error ("stimulus transitions" and "photodiode events" did
not have the same length).

The algorithm has been updated to use consecutive events as
described below:

From @dougo:
1. Every transition from rise/fall or fall/rise represents a
change in the state of the sync square. Call these ‘photodiode events’
2. During the experiment, the sync square changes state every N
frames (where N = 60 by default). Let’s call these ‘sentinel frames’
3. A transition marks a ‘sentinel frame’ if the next photodiode
transition is roughly 1.0s after the previous OR 1.0s before the
next (for N = 60 frames at a 60 Hz display rate). Call these ‘valid transitions’

He also notes that the time delay between each ‘valid transition’
and each ‘sentinel frame’ should be roughly constant (mean
between 20 and 40 ms, std of just a few ms for a given experiment)

Update the logging to show the min, max, and std delay.

No longer use a default value, but throw an error if there are
data issues or the value computed is outside expected range.

allensdk/internal/brain_observatory/time_sync.py

@@ -9,14 +9,11 @@
     cv2 = None
 
 TRANSITION_FRAME_INTERVAL = 60
-SHORT_PHOTODIODE_MIN = 0.1 # seconds
-SHORT_PHOTODIODE_MAX = 0.5 # seconds
-REG_PHOTODIODE_MIN = 1.9 # seconds
-REG_PHOTODIODE_MAX = 2.1 # seconds
-PHOTODIODE_ANOMALY_THRESHOLD = 0.5 # seconds
-LONG_STIM_THRESHOLD = 0.2 # seconds
-ASSUMED_DELAY = 0.0351 # seconds
-MAX_MONITOR_DELAY = 0.07 # seconds
+REG_PHOTODIODE_INTERVAL = 1.0     # seconds
+REG_PHOTODIODE_STD = 0.05    # seconds
+PHOTODIODE_ANOMALY_THRESHOLD = 0.5     # seconds
+LONG_STIM_THRESHOLD = 0.2     # seconds
+MAX_MONITOR_DELAY = 0.07     # seconds
 
 VERSION_1_KEYS = {
     "photodiode": "stim_photodiode",
@@ -53,58 +50,50 @@ def get_keys(sync_dset):
 
 def monitor_delay(sync_dset, stim_times, photodiode_key,
                   transition_frame_interval=TRANSITION_FRAME_INTERVAL,
-                  max_monitor_delay=MAX_MONITOR_DELAY,
-                  assumed_delay=ASSUMED_DELAY):
+                  max_monitor_delay=MAX_MONITOR_DELAY):
     """Calculate monitor delay."""
-    try:
-        transitions = stim_times[::transition_frame_interval]
-        photodiode_events = get_real_photodiode_events(sync_dset, photodiode_key)
-        transition_events = photodiode_events[0:len(transitions)]
-
-        delay = np.mean(transition_events-transitions)
-        logging.info("Calculated monitor delay: %s", delay)
-
-        if delay < 0 or delay > max_monitor_delay:
-            delay = assumed_delay
-            logging.warning("Setting delay to assumed value: %s", delay)
-    except (IndexError, ValueError) as e:
-        logging.error(e)
-        delay = assumed_delay
-        logging.warning("Bad photodiode signal, setting delay to assumed "
-                        "value: %s", delay)
-
+    transitions = stim_times[::transition_frame_interval]
+    photodiode_events = get_real_photodiode_events(sync_dset, photodiode_key)
+    transition_events = photodiode_events[0:len(transitions)]
+
+    delays = transition_events - transitions
+    delay = np.mean(delays)
+    logging.info(f"Calculated monitor delay: {delay}. \n "
+                 f"Max monitor delay: {np.max(delays)}. \n "
+                 f"Min monitor delay: {np.min(delays)}.\n "
+                 f"Std monitor delay: {np.std(delays)}.")
+
+    if delay < 0 or delay > max_monitor_delay:
+        raise ValueError("Delay ({delay}s) falls outside expected value "
+                         "range.")
     return delay
 
 
 def get_photodiode_events(sync_dset, photodiode_key):
     """Returns the photodiode events with the start/stop indicators and
     the window init flash stripped off.
     """
-    all_events = sync_dset.get_events_by_line(photodiode_key)
-    pdr = sync_dset.get_rising_edges(photodiode_key)
-    pdf = sync_dset.get_falling_edges(photodiode_key)
-
-    all_events_sec = all_events/sync_dset.sample_freq
-    pdr_sec = pdr/sync_dset.sample_freq
-    pdf_sec = pdf/sync_dset.sample_freq
-
-    pdf_diff = np.ediff1d(pdf_sec, to_end=0)
-    pdr_diff = np.ediff1d(pdr_sec, to_end=0)
-
-    reg_pd_falling = pdf_sec[(pdf_diff >= REG_PHOTODIODE_MIN) &
-                             (pdf_diff <= REG_PHOTODIODE_MAX)]
-
-    short_pd_rising = pdr_sec[(pdr_diff >= SHORT_PHOTODIODE_MIN) &
-                              (pdr_diff <= SHORT_PHOTODIODE_MAX)]
-
-    first_falling = reg_pd_falling[0]
-    last_falling = reg_pd_falling[-1]
-
-    end_indicators = short_pd_rising[short_pd_rising > last_falling]
-    first_end_indicator = end_indicators[0]
+    all_events = sync_dset.get_events_by_line(photodiode_key, units="seconds")
+    all_events_diff = np.ediff1d(all_events, to_begin=0, to_end=0)
+    all_events_diff_prev = all_events_diff[:-1]
+    all_events_diff_next = all_events_diff[1:]
+
+    min_interval = REG_PHOTODIODE_INTERVAL - REG_PHOTODIODE_STD
+    max_interval = REG_PHOTODIODE_INTERVAL + REG_PHOTODIODE_STD
+    valid_events = np.where(
+        ((all_events_diff_next >= min_interval) &
+         (all_events_diff_next <= max_interval)) |
+        ((all_events_diff_prev >= min_interval) &
+         (all_events_diff_prev <= max_interval)))[0]
+
+    if len(valid_events):
+        first_event_index = valid_events[0]
+        last_event_index = valid_events[-1]
+        pd_events = all_events[first_event_index:last_event_index+1]
+    else:
+        raise ValueError("No valid photodiode events found. Please check "
+                         "the input data for errors. ")
 
-    pd_events =  all_events_sec[(all_events_sec >= first_falling) &
-                                (all_events_sec < first_end_indicator)]
     return pd_events
 
 
@@ -283,7 +272,7 @@ def corrected_stim_timestamps(self):
 
         photodiode_key = self._keys["photodiode"]
         delay = monitor_delay(self.dataset, timestamps, photodiode_key)
-        
+
         return timestamps + delay, delta, delay
 
     @property

allensdk/test/internal/brain_observatory/test_time_sync.py

@@ -7,6 +7,11 @@
 from mock import patch
 from allensdk.internal.brain_observatory import time_sync as ts
 from allensdk.internal.pipeline_modules import run_ophys_time_sync
+from allensdk.brain_observatory.sync_dataset import Dataset
+
+
+ASSUMED_DELAY = 0.0351
+
 
 data_file = resource_filename(__name__, "time_sync_test_data.json")
 test_data = json.load(open(data_file, "r"))
@@ -20,6 +25,24 @@
 MIN_BOUND = .03
 MAX_BOUND = .04
 
+
+class MockSyncDataset(Dataset):
+    """
+    Mock the Dataset class so it doesn't load an h5 file upon
+    initialization.
+    """
+    def __init__(self, data):
+        self.dfile = data
+
+
+def mock_get_real_photodiode_events(data, key):
+    return data
+
+
+def mock_get_events_by_line(line, units="seconds"):
+    return line
+
+
 def calculate_stimulus_alignment(stim_time, valid_twop_vsync_fall):
     stimulus_alignment = np.empty(len(stim_time))
 
@@ -312,7 +335,7 @@ def test_get_corrected_stim_times(stim_data_length, start_delay):
             aligner = ts.OphysTimeAligner("test")
 
     aligner.stim_data_length = stim_data_length
-    with patch.object(ts, "monitor_delay", return_value=ts.ASSUMED_DELAY):
+    with patch.object(ts, "monitor_delay", return_value=ASSUMED_DELAY):
         with patch.object(ts.Dataset, "get_falling_edges",
                           return_value=true_falling) as mock_falling:
             with patch.object(ts.Dataset, "get_rising_edges",
@@ -329,15 +352,15 @@ def test_get_corrected_stim_times(stim_data_length, start_delay):
         assert mock_log.call_count == 2
         assert len(times) == len(true_falling) - 1
         assert delta == len(true_falling) - 1 - stim_data_length
-        assert np.all(times == true_falling[1:] + ts.ASSUMED_DELAY)
+        assert np.all(times == true_falling[1:] + ASSUMED_DELAY)
     elif stim_data_length != len(true_falling):
         mock_rising.assert_called_once()
         mock_log.assert_called_once()
         assert delta == len(true_falling) - stim_data_length
-        assert np.all(times == true_falling + ts.ASSUMED_DELAY)
+        assert np.all(times == true_falling + ASSUMED_DELAY)
     else:
         assert mock_rising.call_count == 0
-        assert np.all(times == true_falling + ts.ASSUMED_DELAY)
+        assert np.all(times == true_falling + ASSUMED_DELAY)
         assert mock_log.call_count == 0
         assert delta == 0
 
@@ -413,23 +436,66 @@ def test_module(input_json):
                            equal_nan=True)
 
 
-@pytest.mark.skipif(data_skip, reason="No sync or data")
-def test_monitor_delay(scientifica_input):
-    sync_file = scientifica_input.pop("sync_file")
-    dset = ts.Dataset(sync_file)
-    stim_times = dset.get_falling_edges("stim_vsync")
-
-    with patch("numpy.mean", side_effect=ValueError()) as mock_mean:
-        delay = ts.monitor_delay(dset, stim_times, "stim_photodiode",
-                                 assumed_delay=20)
-        assert delay == 20
-        mock_mean.assert_called_once()
-
-    with patch.object(ts, "get_real_photodiode_events",
-                      side_effect=IndexError()) as mock_events:
-        delay = ts.monitor_delay(dset, stim_times, "stim_photodiode",
-                                 assumed_delay=30)
-        assert delay == 30
+@pytest.mark.parametrize(
+    "sync_dset,stim_times,transition_interval,expected",
+    [
+        (np.array([1.0, 2.0, 3.0, 4.0, 5.0]),
+         np.array([0.99, 1.99, 2.99, 3.99, 4.99]), 1, 0.01),
+        (np.array([1.0, 2.0, 3.0, 4.0]),
+         np.array([0.95, 2.0, 2.95, 4.0]), 1, 0.025),
+        (np.array([1.0]), np.array([1.0]), 1, 0.0)
+    ],
+)
+def test_monitor_delay(sync_dset, stim_times, transition_interval, expected,
+                       monkeypatch):
+    monkeypatch.setattr(ts, "get_real_photodiode_events",
+                        mock_get_real_photodiode_events)
+    pytest.approx(expected, ts.monitor_delay(sync_dset, stim_times, "key",
+                                             transition_interval))
+
+
+@pytest.mark.parametrize(
+    "sync_dset,stim_times,transition_interval",
+    [
+        (np.array([1.0, 2.0, 3.0]), np.array([]), 1,),    # negative
+        (np.array([1.0, 2.0, 3.0, 4.0]), np.array([1.1, 2.1, 3.1, 4.1]), 1,),  # too big
+    ],
+)
+def test_monitor_delay_raises_error(
+        sync_dset, stim_times, transition_interval,
+        monkeypatch):
+    monkeypatch.setattr(ts, "get_real_photodiode_events",
+                        mock_get_real_photodiode_events)
+    with pytest.raises(ValueError):
+        ts.monitor_delay(sync_dset, stim_times, "key", transition_interval)
+
+
+@pytest.mark.parametrize(
+    "sync_dset,expected",
+    [
+        ([0.25, 0.5, 0.75, 1., 2., 3., 5., 5.75], [1., 2., 3.]),
+        ([1., 2., 3., 4.], [1., 2., 3., 4.]),
+    ],
+)
+def test_get_photodiode_events(sync_dset, expected, monkeypatch):
+    ds = MockSyncDataset(sync_dset)
+    monkeypatch.setattr(ds, "get_events_by_line", mock_get_events_by_line)
+    assert expected == ts.get_photodiode_events(ds, sync_dset)
+
+
+@pytest.mark.parametrize(
+    "sync_dset,expected",
+    [
+        ([], []),
+        ([0.25, 0.25], []),
+    ]
+)
+def test_photodiode_events_error_if_none_found(sync_dset, expected,
+                                               monkeypatch):
+    ds = MockSyncDataset(sync_dset)
+    monkeypatch.setattr(ds, "get_events_by_line", mock_get_events_by_line)
+    with pytest.raises(ValueError):
+        assert expected == ts.get_photodiode_events(ds, sync_dset)
 
 
 @pytest.mark.parametrize("deserialized_pkl,expected", [