diff --git a/README.md b/README.md index 831d39f..dcb6c0e 100644 --- a/README.md +++ b/README.md @@ -81,13 +81,15 @@ projection is what makes a two-dimensional formulation valid: under linear relat the secondary crosses that plane in a straight line, so the three-dimensional question becomes a two-dimensional one about a region. -The probability is the mass of a bivariate Gaussian inside that region, and it is evaluated -two independent ways. Foster and Estes (1992) integrate in polar coordinates with adaptive -quadrature. Alfano (2005a) performs the inner integral analytically with the error function -and applies Simpson's rule to what is left. The two share nothing but the reduction to -principal axes, so agreement between them checks both. Chan's series (2008) and a Monte Carlo -estimator that samples the three-dimensional relative position provide two further checks of -a different kind. +The probability is the mass of a bivariate Gaussian inside that region, and it is evaluated three +independent ways. Foster and Estes (1992) integrate in polar coordinates with adaptive quadrature. +Alfano (2005a) performs the inner integral analytically with the error function and applies +Simpson's rule to what is left. Patera (2001) applies Green's theorem and integrates around the +boundary of the region instead of over its interior, so the region enters only as the curve that +bounds it and an ellipse costs a different curve rather than a different derivation. The three +share nothing but the reduction to principal axes, so agreement between them checks all three. +Chan's series (2008) and a Monte Carlo estimator that samples the three-dimensional relative +position provide two further checks of a different kind. ```python from conjunction_screening import generate_catalog, run_screening @@ -285,7 +287,7 @@ uv run ruff format --check . uv run mypy ``` -209 tests cover 95.75 percent of the 1696 statements in the package. Continuous integration +222 tests cover 95.90 percent of the 1754 statements in the package. Continuous integration runs that same command with `--cov-fail-under=93` on Ubuntu and on Windows, which is the measured figure rounded down and given two points of headroom, so that a platform difference in which branch a filter takes cannot fail a build on its own. @@ -303,14 +305,14 @@ equal-period circles crossing a quarter of a revolution out of phase. The Lipsch the orbit path filter depends on is checked against a finely sampled numerical derivative, because a filter whose bound is not a bound could discard a real conjunction. -Other invariants covered: the state transition matrix is symplectic; the time of closest -approach has zero relative range rate; miss distance is symmetric under swapping the two -objects; the encounter plane projection preserves the magnitude of a perpendicular relative -position; the covariance stays symmetric and positive semi-definite through every stage; -Foster and Alfano agree; Chan is exact for a circular covariance and departs monotonically as -the aspect ratio grows; the combined hard body contains the Minkowski sum of the two bodies in -every direction; the dilution curve rises then falls; and a Monte Carlo estimate agrees with -every analytic value. +Other invariants covered: the state transition matrix is symplectic; the time of closest approach +has zero relative range rate; miss distance is symmetric under swapping the two objects; the +encounter plane projection preserves the magnitude of a perpendicular relative position; the +covariance stays symmetric and positive semi-definite through every stage; Foster, Alfano, and +Patera agree on a disc, and Foster and Patera still agree on an outline that is not one; Chan is +exact for a circular covariance and departs monotonically as the aspect ratio grows; the combined +hard body contains the Minkowski sum of the two bodies in every direction; the dilution curve +rises then falls; and a Monte Carlo estimate agrees with every analytic value. Two rules govern the tolerances. Only values from a converged solve are pinned, and the regression module asserts that every pinned event converged, because the state of a @@ -365,6 +367,11 @@ Methods: Space Center, August 1992. Stable record: [Stanford SearchWorks 13354320](https://searchworks.stanford.edu/view/13354320). Source of the polar quadrature formulation of the two-dimensional probability of collision. +- Patera, R. P. "General Method for Calculating Satellite Collision Probability." Journal of + Guidance, Control, and Dynamics, Vol. 24, No. 4, 2001, pp. 716 to 722. + DOI [10.2514/2.4771](https://doi.org/10.2514/2.4771). Source of the reduction of the area + integral to a contour integral around the boundary of the hard body outline, which is the + analytic method that does not assume the outline is a circle. - Alfano, S. "A Numerical Implementation of Spherical Object Collision Probability." The Journal of the Astronautical Sciences, Vol. 53, No. 1, 2005, pp. 103 to 109. DOI [10.1007/BF03546397](https://doi.org/10.1007/BF03546397). Source of the reduction of diff --git a/docs/design-notes.md b/docs/design-notes.md index a7603b5..22159d4 100644 --- a/docs/design-notes.md +++ b/docs/design-notes.md @@ -111,11 +111,36 @@ each term is one minus a partial sum of a Poisson series, which cancels catastro the ratio of hard body radius to covariance scale is small, and that ratio is small in every real conjunction. +Patera, "General Method for Calculating Satellite Collision Probability", Journal of Guidance, +Control, and Dynamics 24(4), 2001, DOI 10.2514/2.4771, converts the area integral into a +contour integral around the boundary of the region with Green's theorem. In the plane scaled +by the two principal standard deviations the density is isotropic and has a vector potential +in closed form, `(1 - exp(-r^2 / 2)) / r^2` times the perpendicular of the radius vector +measured from the miss point, so the probability becomes an integral in one variable along a +closed curve. That integrand is periodic and analytic in the variable, which is the regime the +trapezoidal rule converges geometrically in, so nothing adaptive is needed: doubling from 128 +nodes meets the 1e-11 tolerance at 256 on every case in the test suite. Patera agrees with +Foster to 2.5e-15 relative over the eight comparison cases, the level Alfano reaches, and +unlike Alfano it keeps that agreement once the region stops being a disc. + +One numerical point decides whether the contour method is usable at screening depths. Its +kernel splits as `1 / r^2` minus `exp(-r^2 / 2) / r^2`, and the first part integrates to the +winding number of the outline about the miss point, which is one when the miss vector lies +inside the body and zero otherwise. Left inside the quadrature it cancels away: on a +probability of 1e-234 the terms it contributes are of order one and the sum returns 5e-21, +which is noise. Subtracted in closed form, what is left is the decaying part alone and the +value is right. The subtracted form divides by zero when the miss point lies exactly on the +outline, where the other form is well behaved, so both are kept and the one whose integrand +has the smaller peak is used. That comparison needs no threshold: the two forms differ by a +term that is known exactly, so the smaller integrand is by construction the one whose sum +cancels less. + The region being integrated over is a disc only when both objects are spheres. A hard body may be given as an ellipsoid instead, in which case the region is the shadow it casts along the -relative velocity and the polar quadrature takes a radial limit that varies with the angle. -Foster and Monte Carlo accept that region; Alfano and Chan reject it, for reasons recorded -under closed limitations below. +relative velocity: the polar quadrature takes a radial limit that varies with the angle, and +the contour integral runs around an ellipse in place of a circle. Foster, Patera, and Monte +Carlo accept that region; Alfano and Chan reject it, for reasons recorded under closed +limitations below. The Monte Carlo estimator samples the three-dimensional relative position from the combined covariance and projects each draw onto the plane normal to the relative velocity, rather than @@ -250,8 +275,9 @@ direction of approach. That is the whole content of the limitation. The probability integral then runs over that ellipse. Foster's polar quadrature needs one change: the upper limit of the radial integral becomes `R(theta)`, and the angle moves from -the inner to the outer variable. Monte Carlo needs one change: the hit test becomes a -quadratic form in the plane rather than a distance in three dimensions. +the inner to the outer variable. Patera's contour integral needs one change: the curve it runs +around becomes the ellipse. Monte Carlo needs one change: the hit test becomes a quadratic +form in the plane rather than a distance in three dimensions. What it cost. Four things, none of them hidden. @@ -260,10 +286,12 @@ expands a non-central chi-square tail about a circular region; in both derivatio enters before the density does, so neither extends by changing a limit. They raise on a non-circular cross section rather than substituting a disc of the same area, because a silently substituted disc would produce a number that looks like a cross check of the Foster -result and is not one. The consequence is that the strongest validation available, two -independent quadratures agreeing to 4.3e-16, exists only for spheres. For a non-spherical body -the cross check is Foster against Monte Carlo, which is four binomial standard errors wide -rather than at machine precision. +result and is not one. Two of the five methods are therefore unavailable for a shaped body, +and that is the cost. It is not a loss of cross validation: Patera's contour integral does not +use the circle either, so it follows an elliptical outline by changing the curve rather than +the formulation, and it agrees with Foster to 1.0e-15 relative on the elliptical cases in the +test suite, the level the two disc methods reach on a disc. Monte Carlo stays the check that +also covers the covariance square root and the projection, four binomial standard errors wide. The combined body is an outer approximation, so the probability it produces is an overestimate. The size of that is measured rather than asserted: for a 30 m by 3 m by 3 m body diff --git a/src/conjunction_screening/__init__.py b/src/conjunction_screening/__init__.py index 9dace6a..68534a5 100644 --- a/src/conjunction_screening/__init__.py +++ b/src/conjunction_screening/__init__.py @@ -21,6 +21,7 @@ CloseApproachSettings, FosterMethod, MonteCarloMethod, + PateraMethod, ProbabilityMethod, ProbabilityResult, find_close_approaches, @@ -51,6 +52,7 @@ "KeplerianElements", "MonteCarloMethod", "OrbitState", + "PateraMethod", "ProbabilityMethod", "ProbabilityResult", "ScreeningConfig", diff --git a/src/conjunction_screening/algorithm/__init__.py b/src/conjunction_screening/algorithm/__init__.py index a28bc9a..337124a 100644 --- a/src/conjunction_screening/algorithm/__init__.py +++ b/src/conjunction_screening/algorithm/__init__.py @@ -39,10 +39,12 @@ CHAN, FOSTER, MONTE_CARLO, + PATERA, AlfanoMethod, ChanMethod, FosterMethod, MonteCarloMethod, + PateraMethod, ProbabilityMethod, ProbabilityResult, ) @@ -61,6 +63,7 @@ "CHAN", "FOSTER", "MONTE_CARLO", + "PATERA", "PATH_FILTER", "PERIGEE_APOGEE_FILTER", "TIME_FILTER", @@ -74,6 +77,7 @@ "FosterMethod", "MaximumProbability", "MonteCarloMethod", + "PateraMethod", "PathSeparation", "ProbabilityMethod", "ProbabilityResult", diff --git a/src/conjunction_screening/algorithm/probability.py b/src/conjunction_screening/algorithm/probability.py index f272924..acac10a 100644 --- a/src/conjunction_screening/algorithm/probability.py +++ b/src/conjunction_screening/algorithm/probability.py @@ -1,6 +1,6 @@ """Two-dimensional probability of collision. -All four implementations answer the same question: given a bivariate Gaussian on +All five implementations answer the same question: given a bivariate Gaussian on the encounter plane with mean at the projected miss vector and covariance equal to the projected combined covariance, what is the mass of that density inside the region the combined hard body covers, which for two spheres is a disc of the @@ -10,10 +10,12 @@ exp(-0.5 * (((x - mx) / sx)^2 + ((y - my) / sy)^2)) dx dy written here in the principal axes of the covariance, where the density -separates. The four methods differ in how the integral is evaluated: +separates. The five methods differ in how the integral is evaluated: * ``FosterMethod`` integrates in polar coordinates over the disc with adaptive quadrature, which is the formulation of Foster and Estes (1992). +* ``PateraMethod`` applies Green's theorem and integrates around the boundary of + the region instead of over its interior, following Patera (2001). * ``AlfanoMethod`` performs the inner integral analytically with the error function and applies Simpson's rule to what is left, following Alfano (2005). * ``ChanMethod`` evaluates Chan's convergent series, which is exact for a @@ -21,15 +23,16 @@ * ``MonteCarloMethod`` samples the three-dimensional relative position error and counts how many straight-line trajectories pass within the hard body radius. It exercises the projection as well as the integral and is the reference the other - three are checked against. + four are checked against. -Foster and Alfano are independent formulations of the same integral, so their +Foster, Alfano, and Patera are independent formulations of the same integral, an +area quadrature, a reduction to one dimension, and a contour integral, so their agreement is a genuine cross validation rather than a restatement. A non-spherical hard body casts an elliptical shadow on the encounter plane -rather than a circular one. Foster and Monte Carlo take that region as it is. -Alfano and Chan reject it, because the circle enters their derivations before -the density does. +rather than a circular one. Foster, Patera, and Monte Carlo take that region as +it is. Alfano and Chan reject it, because the circle enters their derivations +before the density does. """ from __future__ import annotations @@ -41,6 +44,7 @@ from scipy.integrate import dblquad from scipy.special import erf, erfc, gammaln +from conjunction_screening.model.arrays import Matrix, Vector from conjunction_screening.model.encounter import ( EncounterGeometry, PrincipalForm, @@ -52,15 +56,18 @@ "CHAN", "FOSTER", "MONTE_CARLO", + "PATERA", "AlfanoMethod", "ChanMethod", "FosterMethod", "MonteCarloMethod", + "PateraMethod", "ProbabilityMethod", "ProbabilityResult", ] FOSTER: Final[str] = "foster" +PATERA: Final[str] = "patera" ALFANO: Final[str] = "alfano" CHAN: Final[str] = "chan" MONTE_CARLO: Final[str] = "monte-carlo" @@ -183,6 +190,169 @@ def density(distance: float, angle: float) -> float: ) +@dataclass(frozen=True, slots=True) +class PateraMethod: + """Patera's contour integral around the boundary of the hard body outline. + + Scaling each principal axis by its own standard deviation turns the density + into the standard isotropic Gaussian, in which the mass inside a region has a + vector potential: the field ``g(r) (-y, x)`` with + + g(r) = (1 - exp(-r^2 / 2)) / r^2 + + has curl ``exp(-r^2 / 2)``, where ``r`` is measured from the miss vector. + Green's theorem then replaces the integral over the region by one around its + boundary, + + Pc = (1 / (2 pi)) * contour integral of g(r) (x dy - y dx) + + taken anticlockwise. The disc never enters the derivation, so unlike Alfano + and Chan this follows an elliptical outline by changing the curve rather than + the formulation, and it is therefore the analytic cross check that survives a + non-spherical hard body. + + The outline is a smooth closed curve and the integrand is periodic and + analytic in the parameter along it, so the trapezoidal rule converges + geometrically rather than at a fixed order and no adaptive subdivision is + needed. + + Attributes: + relative_tolerance: Convergence threshold on the change between + successive node counts, relative to the newer value. + initial_nodes: Number of boundary points at the first level. + max_nodes: Cap on the node count before the solve is declared + non-converged. + """ + + relative_tolerance: float = 1e-11 + initial_nodes: int = 128 + max_nodes: int = 1 << 18 + + @property + def name(self) -> str: + """Identifier of this method.""" + return PATERA + + def probability(self, encounter: EncounterGeometry) -> ProbabilityResult: + """Integrate around the outline, doubling the node count until it settles.""" + form = principal_axis_form(encounter) + generator, centre = _patera_outline(form) + inside = float(np.linalg.norm(np.linalg.solve(generator, centre))) < 1.0 + winding = 1.0 if inside else 0.0 + nodes = max(self.initial_nodes, 8) + tail_form = _patera_prefers_the_tail_kernel(generator, centre, nodes) + previous = _patera_contour_sum(generator, centre, nodes, winding, tail_form) + while nodes < self.max_nodes: + nodes *= 2 + current = _patera_contour_sum(generator, centre, nodes, winding, tail_form) + change = abs(current - previous) + if change <= self.relative_tolerance * max(current, np.finfo(float).tiny): + return ProbabilityResult( + method=PATERA, + value=float(np.clip(current, 0.0, 1.0)), + error_estimate=change, + converged=True, + detail=f"contour rule converged at {nodes} nodes", + ) + previous = current + return ProbabilityResult( + method=PATERA, + value=float(np.clip(previous, 0.0, 1.0)), + error_estimate=float("nan"), + converged=False, + detail=f"contour rule did not converge by {self.max_nodes} nodes", + ) + + +def _patera_outline(form: PrincipalForm) -> tuple[Matrix, Vector]: + """Return the outline generator and the miss vector, both in units of the sigmas. + + Dividing each principal axis by its own standard deviation is what makes the + kernel a function of the distance alone. The outline becomes ``{L u : |u| = 1}`` + for any square root ``L`` of the scaled shape matrix, and the disc of the + combined radius is the case where that matrix is ``R^2 I``, so the sphere needs + no separate path here either. The sign of the second column is set so that + increasing the parameter traverses the curve anticlockwise, which is the + orientation Green's theorem is written for. + """ + section = form.cross_section + shape = np.eye(2, dtype=np.float64) * form.radius_m**2 if section is None else section.matrix + scaling = np.array([form.sigma_x_m, form.sigma_y_m], dtype=np.float64) + eigenvalues, eigenvectors = np.linalg.eigh(shape / np.outer(scaling, scaling)) + generator = eigenvectors * np.sqrt(np.clip(eigenvalues, 0.0, None)) + if float(np.linalg.det(generator)) < 0.0: + generator = generator * np.array([1.0, -1.0], dtype=np.float64) + centre = np.array([form.mean_x_m, form.mean_y_m], dtype=np.float64) / scaling + return np.asarray(generator, dtype=np.float64), centre + + +def _patera_boundary(generator: Matrix, centre: Vector, nodes: int) -> tuple[Vector, Vector]: + """Return the squared distance to the density centre and the swept area element. + + The swept element is ``x dy - y dx`` divided by the parameter step, which is + twice the area the radius vector covers per unit parameter. + """ + angles = np.arange(nodes, dtype=np.float64) * (2.0 * np.pi / nodes) + cosine, sine = np.cos(angles), np.sin(angles) + point = generator @ np.stack((cosine, sine)) - centre[:, None] + tangent = generator @ np.stack((-sine, cosine)) + squared_radius = np.asarray(point[0] ** 2 + point[1] ** 2, dtype=np.float64) + swept = np.asarray(point[0] * tangent[1] - point[1] * tangent[0], dtype=np.float64) + return squared_radius, swept + + +def _patera_kernel(squared_radius: Vector, tail_form: bool) -> Vector: + """Return the contour kernel at each boundary point. + + The kernel splits as ``1 / r^2`` minus ``exp(-r^2 / 2) / r^2``, and the first + part integrates to the winding number of the outline about the density centre, + which is one or zero and is known without integrating. Two algebraically equal + forms follow and they fail in opposite regimes. The regular form leaves that + part in the quadrature, where it cancels away numerically: on the deep tail + case in the test suite it returns 5e-21 for a probability of 1e-234. The tail + form subtracts it in closed form and is accurate there, but divides by zero + when the density centre lies on the outline, which the regular form handles. + """ + safe = np.where(squared_radius > 0.0, squared_radius, 1.0) + if tail_form: + return np.asarray( + np.where(squared_radius > 0.0, np.exp(-0.5 * safe) / safe, np.inf), dtype=np.float64 + ) + return np.asarray( + np.where(squared_radius > 0.0, -np.expm1(-0.5 * safe) / safe, 0.5), dtype=np.float64 + ) + + +def _patera_prefers_the_tail_kernel(generator: Matrix, centre: Vector, nodes: int) -> bool: + """Return whether the tail form of the kernel is the better conditioned one. + + The two forms differ by a term whose integral is the winding number, so they + approximate the same value and the choice between them is purely numerical. + The one whose integrand is smaller in magnitude is the one whose sum cancels + less, so the peak over the boundary decides it. That rule carries no threshold + and picks the tail form in the tail, where the regular one cancels the answer + away, and the regular form when the density centre sits on the outline, where + the tail one is singular. + """ + squared_radius, swept = _patera_boundary(generator, centre, nodes) + tail = float(np.max(np.abs(_patera_kernel(squared_radius, True) * swept))) + regular = float(np.max(np.abs(_patera_kernel(squared_radius, False) * swept))) + return tail < regular + + +def _patera_contour_sum( + generator: Matrix, centre: Vector, nodes: int, winding: float, tail_form: bool +) -> float: + """Apply the trapezoidal rule to the contour integral at a fixed node count. + + The integrand is periodic, so the trapezoidal rule is the mean over equally + spaced nodes and the endpoint weights of the general rule do not arise. + """ + squared_radius, swept = _patera_boundary(generator, centre, nodes) + total = float(np.mean(_patera_kernel(squared_radius, tail_form) * swept)) + return winding - total if tail_form else total + + def _require_circular_cross_section(form: PrincipalForm, method: str) -> None: """Raise unless the hard body cross section is a disc. diff --git a/tests/test_hardbody.py b/tests/test_hardbody.py index 93e19c6..b9e777f 100644 --- a/tests/test_hardbody.py +++ b/tests/test_hardbody.py @@ -30,6 +30,7 @@ ChanMethod, FosterMethod, MonteCarloMethod, + PateraMethod, ) from conjunction_screening.model.encounter import ( EncounterGeometry, @@ -341,13 +342,59 @@ def test_an_elongated_cross_section_collects_more_mass_pointing_at_the_density() assert towards > across +def test_patera_agrees_with_the_quadrature_over_an_ellipse() -> None: + """The analytic cross check that survives the shape change. + + Alfano and Chan use the circle before the density enters, so neither extends + to an elliptical outline. Patera's contour integral never uses it: the region + appears only as the curve bounding it, so an ellipse costs a different curve + and not a different formulation. That restores a second analytic evaluation + for the shaped body, at the same tolerance the disc case is checked to rather + than at the width of a sampling error. + """ + method = FosterMethod() + contour = PateraMethod() + base = planar_encounter( + miss_distance_m=150.0, sigma_x_m=300.0, sigma_y_m=200.0, hard_body_radius_m=10.0 + ) + for major, minor, orientation in ((20.0, 5.0, 0.3), (40.0, 2.5, 0.0), (12.0, 8.0, 1.2)): + encounter = base.with_cross_section(CrossSection.ellipse(major, minor, orientation)) + area = method.probability(encounter) + boundary = contour.probability(encounter) + assert area.converged + assert boundary.converged + assert boundary.value == pytest.approx(area.value, rel=_QUADRATURE_AGREEMENT) + + +def test_patera_reads_an_elliptical_outline_as_the_ellipse_and_not_its_area() -> None: + """Two outlines of equal area pointed differently give different answers. + + An equal-area disc substitution, which is what Chan would have to make, cannot + tell these two apart. The contour method separates them because the curve it + integrates around is the outline itself, and the ordering it produces is the + one the geometry requires: the outline reaching towards the density peak + collects more mass. + """ + method = PateraMethod() + base = planar_encounter( + miss_distance_m=400.0, sigma_x_m=300.0, sigma_y_m=300.0, hard_body_radius_m=10.0 + ) + towards = method.probability(base.with_cross_section(CrossSection.ellipse(40.0, 2.5))) + across = method.probability( + base.with_cross_section(CrossSection.ellipse(40.0, 2.5, 0.5 * np.pi)) + ) + disc = method.probability(base.with_cross_section(CrossSection.disc(10.0))) + assert towards.value > disc.value > across.value + + def test_monte_carlo_agrees_with_the_quadrature_over_an_ellipse() -> None: - """The only independent check available once Alfano and Chan are out. + """A check of the elliptical path that goes through the plane construction too. - Foster and Alfano cross validate each other on a disc, and neither of the two - series methods extends to an elliptical outline, so the elliptical path is - checked against sampling instead. Tolerance: four binomial standard errors - computed from the estimate itself. + Foster and Patera cross validate each other on the outline analytically, but + both start from the same principal axis reduction of an encounter that is + already planar. Sampling in three dimensions and projecting covers the + covariance square root and the projection as well. Tolerance: four binomial + standard errors computed from the estimate itself. """ encounter = planar_encounter( miss_distance_m=150.0, sigma_x_m=300.0, sigma_y_m=200.0, hard_body_radius_m=10.0 diff --git a/tests/test_probability.py b/tests/test_probability.py index 56d0da4..0b071e2 100644 --- a/tests/test_probability.py +++ b/tests/test_probability.py @@ -8,6 +8,10 @@ differences in the platform's exponential and error function implementations, each of which can differ by an ulp between operating systems. +Patera's contour rule is asked for the same 1e-11 and is compared on the same +band. It converges geometrically rather than at a fixed order, so a level that +meets the threshold is usually several digits past it. + Chan's series is exact only when the in-plane covariance is circular. Where it is circular the tolerance is again the quadrature tolerance; where it is not, the disagreement is the equal-area approximation error and is checked to be present @@ -30,10 +34,12 @@ CHAN, FOSTER, MONTE_CARLO, + PATERA, AlfanoMethod, ChanMethod, FosterMethod, MonteCarloMethod, + PateraMethod, ProbabilityMethod, ) from conjunction_screening.model.encounter import EncounterGeometry, planar_encounter @@ -81,13 +87,104 @@ def test_foster_and_alfano_agree(case: tuple[str, float, float, float, float, fl assert alfano.value == pytest.approx(foster.value, rel=_QUADRATURE_AGREEMENT) +@pytest.mark.parametrize("case", _CASES, ids=[case[0] for case in _CASES]) +def test_foster_and_patera_agree(case: tuple[str, float, float, float, float, float]) -> None: + """A third route to the same integral lands on the same answer. + + Foster covers the region with an adaptive area quadrature. Patera never + enters the region at all: Green's theorem turns the mass inside it into a + contour integral around its boundary, evaluated with the trapezoidal rule on + a periodic integrand. The two share the principal axis reduction and nothing + else, so this is the second independent check of Foster and the first that + survives a non-circular outline. + """ + encounter = _encounter(case) + foster = FosterMethod().probability(encounter) + patera = PateraMethod().probability(encounter) + assert foster.converged + assert patera.converged + assert patera.value == pytest.approx(foster.value, rel=_QUADRATURE_AGREEMENT) + + +def test_patera_matches_the_closed_form_for_a_centred_disc() -> None: + """A disc centred on the density has a closed-form mass, and it is not a quadrature. + + With a circular covariance and no miss distance the integral is the Rayleigh + distribution function, ``1 - exp(-R^2 / (2 s^2))``, in closed form. Checking + against it fixes the contour orientation, the area element, and the kernel + normalisation at once, without any other method being involved. + """ + method = PateraMethod() + for sigma_m, radius_m in ((250.0, 10.0), (100.0, 90.0), (40.0, 400.0)): + encounter = planar_encounter( + miss_distance_m=0.0, + sigma_x_m=sigma_m, + sigma_y_m=sigma_m, + hard_body_radius_m=radius_m, + ) + expected = -float(np.expm1(-0.5 * (radius_m / sigma_m) ** 2)) + result = method.probability(encounter) + assert result.converged + assert result.value == pytest.approx(expected, rel=1e-13) + + +def test_patera_resolves_a_probability_far_below_the_dismissal_threshold() -> None: + """The deep tail is where the two forms of the contour kernel differ. + + This geometry is the shape of the faintest event a screening run produces: a + miss of eighty standard deviations across the tight in-plane direction, whose + probability is two hundred orders below one. Leaving the winding term inside + the quadrature there sums quantities of order one that must cancel to that + value, and returns 5e-21 rather than the answer, so a value agreeing with + Foster to several significant figures is evidence the subtraction is being + done in closed form. + + Tolerance: 1e-6 relative. Foster's adaptive quadrature is itself working two + hundred orders below one here and does not hold the 1e-11 it reports at that + depth. What is being checked is agreement to several figures rather than a + disagreement by orders of magnitude. + """ + encounter = planar_encounter( + miss_distance_m=3_000.0, + sigma_x_m=4_449.0, + sigma_y_m=35.8, + hard_body_radius_m=6.4, + orientation_rad=0.4, + ) + foster = FosterMethod().probability(encounter) + patera = PateraMethod().probability(encounter) + assert patera.converged + assert 0.0 < patera.value < 1e-200 + assert patera.value == pytest.approx(foster.value, rel=1e-6) + + +def test_patera_reports_a_starved_node_count_as_not_converged() -> None: + """A contour resolved by eight points is not an answer, and it says so. + + The outline here is large enough against the covariance that the integrand + varies by orders of magnitude around it, so eight and sixteen nodes cannot + agree to the requested tolerance and the refinement runs out of levels. + """ + encounter = planar_encounter( + miss_distance_m=800.0, sigma_x_m=100.0, sigma_y_m=100.0, hard_body_radius_m=100.0 + ) + result = PateraMethod(initial_nodes=8, max_nodes=16).probability(encounter) + assert not result.converged + assert "did not converge" in result.detail + + @pytest.mark.parametrize("case", _CASES, ids=[case[0] for case in _CASES]) def test_every_method_returns_a_probability( case: tuple[str, float, float, float, float, float], ) -> None: """All results lie in the unit interval and report their own accuracy.""" encounter = _encounter(case) - methods: tuple[ProbabilityMethod, ...] = (FosterMethod(), AlfanoMethod(), ChanMethod()) + methods: tuple[ProbabilityMethod, ...] = ( + FosterMethod(), + PateraMethod(), + AlfanoMethod(), + ChanMethod(), + ) for method in methods: result = method.probability(encounter) assert 0.0 <= result.value <= 1.0 @@ -233,6 +330,7 @@ def test_monte_carlo_reports_a_starved_estimate_as_not_converged() -> None: def test_methods_expose_their_names() -> None: """Every method reports the identifier used to key comparison tables.""" assert FosterMethod().name == FOSTER + assert PateraMethod().name == PATERA assert AlfanoMethod().name == ALFANO assert ChanMethod().name == CHAN assert MonteCarloMethod().name == MONTE_CARLO