diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index e094ce5..462060b 100644 --- a/chandra_aca/aca_image.py +++ b/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. diff --git a/chandra_aca/tests/test_aca_image.py b/chandra_aca/tests/test_aca_image.py index 7c7f31e..c55ebd1 100644 --- a/chandra_aca/tests/test_aca_image.py +++ b/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]))