Add test mode for comparison to C implementation

chandra_aca/aca_image.py

@@ -369,8 +369,9 @@ def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0, seed=None):
 
         self.flicker_mean_time = flicker_mean_time
         self.flicker_scale = flicker_scale
+        self.test_idx = 1 if seed == -1 else 0
 
-        if seed is not None:
+        if seed is not None and seed != -1:
             np.random.seed(seed)
             _numba_random_seed(seed)
 
@@ -395,16 +396,21 @@ def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0, seed=None):
                                                 self.flicker_vals0) - 1
 
         # Create an array of time (secs) until next flicker for each pixel
-        # This is drawing from an exponential distribution.  For the initial
-        # time assume the flickering is randomly phased within that interval.
-        phase = np.random.uniform(0.0, 1.0, size=self.flicker_n_vals)
-        rand_unifs = np.random.uniform(0.0, 1.0, size=self.flicker_n_vals)
-        t_flicker = -np.log(1.0 - rand_unifs) * flicker_mean_time
+        if seed == -1:
+            # Special case of testing, use a constant flicker_mean_time initially
+            t_flicker = np.ones(shape=(self.flicker_n_vals,)) * flicker_mean_time
+            phase = 1.0
+        else:
+            # This is drawing from an exponential distribution.  For the initial
+            # time assume the flickering is randomly phased within that interval.
+            phase = np.random.uniform(0.0, 1.0, size=self.flicker_n_vals)
+            rand_unifs = np.random.uniform(0.0, 1.0, size=self.flicker_n_vals)
+            t_flicker = -np.log(1.0 - rand_unifs) * flicker_mean_time
+
         self.flicker_times = t_flicker * phase
 
     def flicker_update(self, dt, use_numba=True):
-        """
-        Propagate the image forward by ``dt`` seconds and update any pixels
+        """Propagate the image forward by ``dt`` seconds and update any pixels
         that have flickered during that interval.
 
         This has the option to use one of two implementations.  The default is
@@ -418,7 +424,9 @@ def flicker_update(self, dt, use_numba=True):
             self.flicker_init()
 
         if use_numba:
-            _flicker_update_numba(dt, len(self.flicker_vals),
+            _flicker_update_numba(dt,
+                                  len(self.flicker_vals),
+                                  self.test_idx,
                                   self.flicker_vals0,
                                   self.flicker_vals,
                                   self.flicker_mask_vals,
@@ -428,6 +436,8 @@ def flicker_update(self, dt, use_numba=True):
                                   self.flicker_cdfs,
                                   self.flicker_scale,
                                   self.flicker_mean_time)
+            if self.test_idx > 0:
+                self.test_idx += 1
         else:
             self._flicker_update_vectorized(dt)
 
@@ -475,7 +485,7 @@ def _numba_random_seed(seed):
 
 
 @numba.jit(nopython=True)
-def _flicker_update_numba(dt, nvals,
+def _flicker_update_numba(dt, nvals, test_idx,
                           flicker_vals0,
                           flicker_vals,
                           flicker_mask_vals,
@@ -501,10 +511,16 @@ def _flicker_update_numba(dt, nvals,
         if flicker_times[ii] > 0:
             continue
 
-        # Random uniform used for (1) distribution of flickering amplitude
-        # via the CDFs and (2) distribution of time to next flicker.
-        rand_ampl = np.random.uniform(0.0, 1.0)
-        rand_time = np.random.uniform(0.0, 1.0)
+        if test_idx > 0:
+            # Deterministic and reproducible but bouncy sequence that is reproducible in C
+            # (which has a different random number generator).
+            rand_ampl = np.abs(np.sin(float(ii + test_idx)))
+            rand_time = np.abs(np.cos(float(ii + test_idx)))
+        else:
+            # Random uniform used for (1) distribution of flickering amplitude
+            # via the CDFs and (2) distribution of time to next flicker.
+            rand_ampl = np.random.uniform(0.0, 1.0)
+            rand_time = np.random.uniform(0.0, 1.0)
 
         # Determine the new value after flickering and set in array view.
         # First get the right CDF from the list of CDFs based on the pixel value.

chandra_aca/tests/test_aca_image.py

@@ -362,3 +362,39 @@ def test_flicker_no_seed():
         b.flicker_update(100.0)
 
     assert np.any(a != b)
+
+
+def test_flicker_test_sequence():
+    """Test the deterministic testing sequence that allows for
+    cross-implementation comparison.  This uses the seed=-1 pathway.
+    """
+    a = ACAImage(np.linspace(0, 800, 9).reshape(3, 3))
+    a.flicker_init(seed=-1, flicker_mean_time=15)
+    dt = 10
+    assert np.all(np.round(a) == [[0, 100, 200],
+                                  [300, 400, 500],
+                                  [600, 700, 800]])
+
+    a.flicker_update(dt)
+    assert np.all(np.round(a) == [[0, 100, 200],
+                                  [300, 400, 500],
+                                  [600, 700, 800]])
+
+    a.flicker_update(dt)
+    assert np.all(np.round(a) == [[0, 80, 229],
+                                  [447, 374, 531],
+                                  [952, 693, 815]])
+    a.flicker_update(dt)
+    assert np.all(np.round(a) == [[0, 80, 229],
+                                  [278, 374, 531],
+                                  [592, 693, 815]])
+
+    a.flicker_update(dt)
+    assert np.all(np.round(a) == [[0, 80, 185],
+                                  [278, 374, 531],
+                                  [592, 693, 815]])
+
+    assert np.allclose(a.flicker_times,
+                       np.array([15., 49.06630191, 48.34713118, 38.34713118,
+                                 28.34713118, 1.03039723, 26.30875397, 16.30875397,
+                                 7.40192939]))