GH-1404: Handle monitor delay errors

In order to compute monitor delay, a sync square monitored
by a photodiode changes state every 60 frames ('sentinel
event'). Each transition from a rise/fall to a fall/rise
('photodiode event') represents a change in the state of
the sync square.
Previously only consecutive falling edges of the photodiode
were used to determine whether the photodiode events triggered
by sentinel frames had started. Only a falling edge was allowed
to be the final event in a photodiode stream, even though it
is possible that the last valid event could be a rising edge.
This led sometimes to missing valid events, and resulted
in an error (because the number of sentinel frames should match
the number of photodiode events).
The algorithm has been updated to examine consecutive events
instead of relying on a specific edge type. A transition marks
a sentinel frame if the next photodiode transition is roughly
1.0s after the previous or 1.0s before the next. Due to noisy
data we can't just look for all contiguous events of this type.
Instad we mark the beginning and end of the photodiode event
stream as follows:
1. The first valid photodiode event is the first event where
the time between the next two consecutive events are ~1.0s
2. The last valid photodiode event is the final event where
the last event where the time between the previous two
consecutive events are ~1.0s

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

@@ -1,3 +1,5 @@
+from collections import deque
+
 import numpy as np
 import h5py
 from allensdk.brain_observatory.sync_dataset import Dataset
@@ -9,14 +11,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.0215 # 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 +52,101 @@ 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.
+def _find_last_n(arr, n, cond):
     """
-    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)]
+    Find the final index where the prior `n` values in an array meet
+    the condition `cond` (inclusive).
+    Parameters
+    ==========
+    arr: numpy.1darray
+    n: int
+    cond: Callable that returns True if condition is met, False
+    otherwise. Should be able to be applied to the array elements  
+    without any additional arguments.
+    """
+    reversed_ix = _find_n(arr[::-1], n, cond)
+    if reversed_ix is not None:
+        reversed_ix = len(arr) - reversed_ix - 1
+    return reversed_ix
 
-    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]
+def _find_n(arr, n, cond):
+    """
+    Find the index where the next `n` values in an array meet the 
+    condition `cond` (inclusive).
+    Parameters
+    ==========
+    arr: numpy.1darray
+    n: int
+    cond: Callable that returns True if condition is met, False
+    otherwise. Should be able to be applied to the array elements  
+    without any additional arguments.
+    """
+    if len(arr) < n:
+        return None
+    queue = deque(np.apply_along_axis(cond, 0, arr[:n]), maxlen=n)
+    i = 0
+    while queue.count(True) < n:
+        try:
+            i += 1
+            queue.append(cond(arr[i+n-1]))
+        except IndexError:
+            return None
+    return i
 
-    end_indicators = short_pd_rising[short_pd_rising > last_falling]
-    first_end_indicator = end_indicators[0]
 
-    pd_events =  all_events_sec[(all_events_sec >= first_falling) &
-                                (all_events_sec < first_end_indicator)]
+def get_photodiode_events(sync_dset, photodiode_key):
+    """Returns the photodiode events with the start/stop indicators and
+    the window init flash stripped off. These transitions occur roughly
+    ~1.0s apart, since the sync square changes state every N frames
+    (where N = 60, and frame rate is 60 Hz). Because there are no
+    markers for when the first transition of this type started, we
+    estimate based on the event intervals. For the first valid event,
+    find the first two events that both meet the following criteria:
+        The next event occurs ~1.0s later
+    First the last valid event, find the first two events that both meet
+    the following criteria:
+        The last valid event occured ~1.0s before
+    """
+    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
+    if not len(all_events):
+        raise ValueError("No photodiode events found. Please check "
+                         "the input data for errors. ")
+    first_valid_index = _find_n(
+        all_events_diff_next, 2,
+        lambda x: (x >= min_interval) & (x <= max_interval))
+    last_valid_index = _find_last_n(
+        all_events_diff_prev, 2,
+        lambda x: (x >= min_interval) & (x <= max_interval))
+    if first_valid_index is None:
+        raise ValueError("Can't find valid start event")
+    if last_valid_index is None:
+        raise ValueError("Can't find valid end event")
+    pd_events = all_events[first_valid_index:last_valid_index+1]
     return pd_events
 
 
@@ -283,7 +325,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,99 @@ 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([0.9, 1.9, 2.9]), 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(
+    "arr,cond,n,expected",
+    [
+        (np.array([1, 1, 1, 2]), lambda x: x < 2, 3, 0),
+        (np.array([2, 1, 1, 1]), lambda x: x < 2, 3, 1),
+        (np.array([1, 2, 2, 1, 1]), lambda x: x >= 1, 2, 0),
+        (np.array([]), lambda x: x < 1, 1, None),
+        (np.array([1, 2, 3]), lambda x: x < 3, 4, None),
+        (np.array([1, 2, 2, 3, 2]), lambda x: x == 2, 3, None),
+    ]
+)
+def test_find_n(arr, cond, n, expected):
+    assert expected == ts._find_n(arr, n, cond)
+
+
+@pytest.mark.parametrize(
+    "arr,cond,n,expected",
+    [
+        (np.array([1, 1, 1, 2]), lambda x: x < 2, 3, 2),
+        (np.array([2, 1, 1, 1]), lambda x: x < 2, 3, 3),
+        (np.array([1, 2, 2, 1, 1]), lambda x: x >= 1, 2, 4),
+        (np.array([]), lambda x: x < 1, 1, None),
+        (np.array([1, 2, 3]), lambda x: x < 3, 4, None),
+        (np.array([1, 2, 2, 3, 2]), lambda x: x == 2, 3, None),
+    ]
+)
+def test_find_last_n(arr, cond, n, expected):
+    assert expected == ts._find_last_n(arr, n, cond)
+
+
+@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.]),
+        ([0.25, 1., 2., 2.1, 2.2, 3., 4., 5.], [3., 4., 5.]),   # false alarm start
+        ([0.25, 1., 2., 3., 4., 4.5, 5.1, 6.1], [1., 2., 3., 4.]),  # false alarm end
+    ],
+)
+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)
+    np.testing.assert_array_equal(
+        expected, ts.get_photodiode_events(ds, sync_dset))
+
+
+@pytest.mark.parametrize(
+    "sync_dset,",
+    [
+        ([]),
+        ([0.25, 0.25]),
+        ([1., 2.]),
+    ]
+)
+def test_photodiode_events_error_if_none_found(sync_dset, monkeypatch):
+    ds = MockSyncDataset(sync_dset)
+    monkeypatch.setattr(ds, "get_events_by_line", mock_get_events_by_line)
+    with pytest.raises(ValueError):
+        ts.get_photodiode_events(ds, sync_dset)
 
 
 @pytest.mark.parametrize("deserialized_pkl,expected", [