diff --git a/src/pyrecest/filters/daum_huang_particle_filter.py b/src/pyrecest/filters/daum_huang_particle_filter.py index 3c142ddcc..5aeedd167 100644 --- a/src/pyrecest/filters/daum_huang_particle_filter.py +++ b/src/pyrecest/filters/daum_huang_particle_filter.py @@ -155,6 +155,7 @@ def gaussian_flow_affine_increment( ): """Apply one exact affine Gaussian bridge increment for a linear model.""" jitter = _validate_nonnegative_float(jitter, "jitter") + delta_lambda = _validate_nonnegative_float(delta_lambda, "delta_lambda") X = _as_particle_matrix_np(particles) m0, P0, H, y, R = _validate_linear_bridge_np( mean, @@ -191,6 +192,7 @@ def gaussian_bridge_moments( ): """Return Gaussian moments after a likelihood-power increment.""" jitter = _validate_nonnegative_float(jitter, "jitter") + delta_lambda = _validate_nonnegative_float(delta_lambda, "delta_lambda") m0, P0, H, y, R = _validate_linear_bridge_np( mean, covariance, @@ -575,10 +577,12 @@ def _measurement_matrix(model): def _try_batch_measurement(function, X): try: - values = to_numpy(function(asarray(X))) + raw_values = np.asarray(to_numpy(function(asarray(X)))) except (TypeError, ValueError, NotImplementedError, IndexError): return None - values = np.asarray(values, dtype=float) + if np.iscomplexobj(raw_values): + raise ValueError("measurement value must be real-valued.") + values = np.asarray(raw_values, dtype=float) if values.ndim == 1 and X.shape[0] == 1: return values.reshape(1, -1) if values.ndim == 2 and values.shape[0] == X.shape[0]: @@ -588,10 +592,12 @@ def _try_batch_measurement(function, X): def _try_batch_jacobian(function, X): try: - values = to_numpy(function(asarray(X))) + raw_values = np.asarray(to_numpy(function(asarray(X)))) except (TypeError, ValueError, NotImplementedError, IndexError): return None - values = np.asarray(values, dtype=float) + if np.iscomplexobj(raw_values): + raise ValueError("measurement jacobian must be real-valued.") + values = np.asarray(raw_values, dtype=float) if values.ndim == 2 and X.shape[0] == 1: return values.reshape(1, values.shape[0], values.shape[1]) if values.ndim == 3 and values.shape[0] == X.shape[0]: @@ -650,7 +656,10 @@ def _weighted_mean_np(particles, weights): def _as_particle_matrix_np(value): if _contains_masked_value(value): raise ValueError("particles must not contain masked values.") - X = np.asarray(to_numpy(value), dtype=float) + raw = np.asarray(to_numpy(value)) + if np.iscomplexobj(raw): + raise ValueError("particles must be real-valued.") + X = np.asarray(raw, dtype=float) if X.ndim == 1: X = X[None, :] if X.ndim != 2: @@ -665,7 +674,10 @@ def _as_particle_matrix_np(value): def _as_weights_np(value, n_particles: int): if _contains_masked_value(value): raise ValueError("weights must not contain masked values.") - weights = np.asarray(to_numpy(value), dtype=float).reshape(-1) + raw = np.asarray(to_numpy(value)) + if np.iscomplexobj(raw): + raise ValueError("weights must be real-valued.") + weights = np.asarray(raw, dtype=float).reshape(-1) if weights.shape != (n_particles,): raise ValueError("weights must have one entry per particle.") if not np.all(np.isfinite(weights)): @@ -685,7 +697,10 @@ def _as_weights_np(value, n_particles: int): def _as_vector_np(value, name): if _contains_masked_value(value): raise ValueError(f"{name} must not contain masked values.") - vector = np.asarray(to_numpy(value), dtype=float) + raw = np.asarray(to_numpy(value)) + if np.iscomplexobj(raw): + raise ValueError(f"{name} must be real-valued.") + vector = np.asarray(raw, dtype=float) if vector.ndim == 0: vector = vector.reshape(1) vector = vector.reshape(-1) if vector.ndim == 1 else vector @@ -699,7 +714,10 @@ def _as_vector_np(value, name): def _as_matrix_np(value, name, *, scalar_dim: int | None = None): if _contains_masked_value(value): raise ValueError(f"{name} must not contain masked values.") - matrix = np.asarray(to_numpy(value), dtype=float) + raw = np.asarray(to_numpy(value)) + if np.iscomplexobj(raw): + raise ValueError(f"{name} must be real-valued.") + matrix = np.asarray(raw, dtype=float) if matrix.ndim == 0 and scalar_dim == 1: matrix = matrix.reshape(1, 1) if matrix.ndim != 2: @@ -715,7 +733,10 @@ def _lambda_deltas_np(n_steps, step_schedule): return np.full(n_steps, 1.0 / float(n_steps)) if _contains_masked_value(step_schedule): raise ValueError("step_schedule must not contain masked values.") - deltas = np.asarray(step_schedule, dtype=float).reshape(-1) + raw = np.asarray(step_schedule) + if np.iscomplexobj(raw): + raise ValueError("step_schedule must be real-valued.") + deltas = np.asarray(raw, dtype=float).reshape(-1) if deltas.size == 0: raise ValueError("step_schedule must not be empty.") if np.any(deltas <= 0.0): @@ -731,7 +752,7 @@ def _lambda_deltas_np(n_steps, step_schedule): def _regularize_cov_np(covariance, jitter): - covariance = _symmetrize_np(np.asarray(covariance, dtype=float)) + covariance = np.asarray(covariance, dtype=float) if covariance.ndim != 2 or covariance.shape[0] != covariance.shape[1]: raise ValueError("covariance must be square.") if covariance.shape[0] == 0: @@ -739,6 +760,17 @@ def _regularize_cov_np(covariance, jitter): if not np.all(np.isfinite(covariance)): raise ValueError("covariance must be finite.") + matrix_scale = max(float(np.max(np.abs(covariance))), 1.0) + symmetry_tolerance = ( + 10.0 + * np.finfo(float).eps + * max(covariance.shape[0], 1) + * matrix_scale + ) + if float(np.max(np.abs(covariance - covariance.T))) > symmetry_tolerance: + raise ValueError("covariance must be symmetric.") + covariance = _symmetrize_np(covariance) + eigenvalues = np.linalg.eigvalsh(covariance) spectral_scale = max(float(np.max(np.abs(eigenvalues))), 1.0) tolerance = ( diff --git a/tests/filters/test_daum_huang_particle_filter.py b/tests/filters/test_daum_huang_particle_filter.py index 53dc01f96..0864db7d4 100644 --- a/tests/filters/test_daum_huang_particle_filter.py +++ b/tests/filters/test_daum_huang_particle_filter.py @@ -71,6 +71,84 @@ def test_affine_increment_matches_bridge_moments(self): to_numpy(actual_covariance), to_numpy(expected_covariance), atol=1e-10 ) + def test_bridge_rejects_nonsymmetric_covariances(self): + mean = array([0.0, 0.0]) + identity_covariance = array([[1.0, 0.0], [0.0, 1.0]]) + nonsymmetric_covariance = array([[1.0, 0.5], [0.0, 1.0]]) + + with self.assertRaisesRegex(ValueError, "symmetric"): + gaussian_bridge_moments( + mean, + nonsymmetric_covariance, + array([[1.0, 0.0]]), + array([0.0]), + array([[1.0]]), + 1.0, + jitter=0.0, + ) + + with self.assertRaisesRegex(ValueError, "symmetric"): + gaussian_bridge_moments( + mean, + identity_covariance, + identity_covariance, + array([0.0, 0.0]), + nonsymmetric_covariance, + 1.0, + jitter=0.0, + ) + + def test_bridge_rejects_complex_numpy_inputs_without_truncation(self): + real_vector = array([0.0]) + real_matrix = array([[1.0]]) + + invalid_cases = [ + (np.array([1.0 + 2.0j]), real_matrix), + (real_vector, np.array([[1.0 + 2.0j]])), + ] + for mean, covariance in invalid_cases: + with self.subTest(mean=mean, covariance=covariance): + with self.assertRaisesRegex(ValueError, "real-valued"): + gaussian_bridge_moments( + mean, + covariance, + real_matrix, + real_vector, + real_matrix, + 1.0, + jitter=0.0, + ) + + def test_bridge_exponent_must_be_finite_nonnegative_scalar(self): + mean = array([0.0]) + covariance = array([[1.0]]) + measurement_matrix = array([[1.0]]) + measurement = array([0.0]) + measurement_noise = array([[1.0]]) + + for delta_lambda in (float("nan"), float("inf"), -0.1, True, [1.0]): + with self.subTest(delta_lambda=delta_lambda): + with self.assertRaisesRegex(ValueError, "delta_lambda"): + gaussian_bridge_moments( + mean, + covariance, + measurement_matrix, + measurement, + measurement_noise, + delta_lambda, + ) + + with self.assertRaisesRegex(ValueError, "delta_lambda"): + gaussian_flow_affine_increment( + array([[0.0], [1.0]]), + mean, + covariance, + measurement_matrix, + measurement, + measurement_noise, + float("nan"), + ) + def test_edh_filter_linear_update_matches_bridge_moments(self): particles = array( [