Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 78 additions & 1 deletion autolens/point/model/analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,12 @@
It also manages result output (``ResultPoint``) and on-the-fly visualisation
(``VisualizerPoint``).
"""

import logging
import os

import numpy as np
from autonerves.test_mode import is_test_mode

import autofit as af
import autogalaxy as ag
Expand All @@ -32,6 +37,8 @@
from autolens.point.model.visualizer import VisualizerPoint
from autolens.point.solver import PointSolver

logger = logging.getLogger(__name__)


class AnalysisPoint(AgAnalysis, AnalysisLens):
Visualizer = VisualizerPoint
Expand Down Expand Up @@ -84,7 +91,9 @@ def __init__(
time-delays.
solver
Solves the lens equation in order to determine the image-plane positions of a point source by ray-tracing
triangles to and from the source-plane.
triangles to and from the source-plane. At construction, a warning is logged if observed positions lie
within two initial triangle scales plus three position sigmas of a grid edge, or outside it. This
diagnostic does not establish image completeness across the model prior and never changes the grid.
fit_positions_cls
The class used to fit the positions of the point source dataset, which could be an image-plane or
source-plane chi-squared.
Expand Down Expand Up @@ -114,6 +123,74 @@ def __init__(
self.fit_flux_cls = fit_flux_cls
self.fit_time_delays_cls = fit_time_delays_cls
self.title_prefix = title_prefix
self._check_solver_extent()

def _check_solver_extent(self):
"""Diagnose observed-position coverage once, outside likelihood tracing.

The margin is two initial triangle scales plus three position sigmas.
This is a data-based hint, not a completeness test over the model prior.
Custom solvers without rectangular bounds and absent/invalid diagnostic
inputs are left to their existing validation paths.
"""
positions = getattr(self.dataset, "positions", None)
noise = getattr(self.dataset, "positions_noise_map", None)
if positions is None or noise is None or self.solver is None:
return
try:
bounds = np.asarray(
[
self.solver.y_min,
self.solver.y_max,
self.solver.x_min,
self.solver.x_max,
],
dtype=float,
)
scale = float(self.solver.scale)
positions = np.asarray(getattr(positions, "array", positions), dtype=float)
noise = np.asarray(getattr(noise, "array", noise), dtype=float)
if positions.ndim != 2 or positions.shape[1] != 2 or not len(positions):
return
noise = np.broadcast_to(noise, (len(positions),))
except (AttributeError, TypeError, ValueError):
return
if (
not np.all(np.isfinite(bounds))
or not np.isfinite(scale)
or scale <= 0
or not np.all(np.isfinite(positions))
or not np.all(np.isfinite(noise))
or np.any(noise < 0)
):
return
lower, upper = bounds[[0, 2]], bounds[[1, 3]]
if np.any(upper <= lower):
return
margin = (2 * scale + 3 * noise)[:, None]
if np.any(positions - lower <= margin) or np.any(upper - positions <= margin):
logger.warning(
"PointSolver extent %s is too close to, or excludes, observed "
"positions for dataset %r (margin: 2 * initial scale + 3 * "
"position sigma). Images near or beyond the grid edge may be "
"missed. Enlarge the solver grid and check the model prior: "
"covering observed positions does not guarantee image "
"completeness across the model prior.",
tuple(bounds.tolist()),
getattr(self.dataset, "name", None),
)
elif not is_test_mode() and os.environ.get("PYAUTO_SMALL_DATASETS") != "1":
midpoint = (lower + upper) / 2
required = np.max(np.abs(positions - midpoint) + margin, axis=0)
if np.all((upper - lower) / 2 > 3 * required):
logger.info(
"PointSolver extent for dataset %r is more than 3 times "
"the observed-position envelope on both axes, including "
"the safety margin. At fixed initial scale, initial grid "
"work grows with grid area. Consider a smaller grid only "
"after checking image completeness across the model prior.",
getattr(self.dataset, "name", None),
)

def log_likelihood_function(self, instance):
"""
Expand Down
139 changes: 139 additions & 0 deletions test_autolens/point/model/test_analysis_point_extent.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,139 @@
"""Construction-time extent diagnostics use only observed NumPy data."""

import logging
from types import SimpleNamespace

import numpy as np
import pytest

import autolens as al

LOGGER = "autolens.point.model.analysis"


@pytest.fixture(autouse=True)
def production_mode(monkeypatch, caplog):
monkeypatch.setenv("PYAUTO_TEST_MODE", "0")
monkeypatch.delenv("PYAUTO_SMALL_DATASETS", raising=False)
caplog.set_level(logging.INFO, logger=LOGGER)


def make_analysis(positions=((0, 0),), noise=0.1, bounds=(-2, 2, -2, 2), scale=0.2):
dataset = al.PointDataset(
name="point_0", positions=positions, positions_noise_map=noise
)
solver = al.PointSolver(
y_min=bounds[0],
y_max=bounds[1],
x_min=bounds[2],
x_max=bounds[3],
scale=scale,
pixel_scale_precision=0.001,
)
return al.AnalysisPoint(dataset=dataset, solver=solver, use_jax=False)


def records(caplog):
return [r for r in caplog.records if r.name == LOGGER]


def test_safe_extent(caplog):
make_analysis()
assert records(caplog) == []


@pytest.mark.parametrize("axis", [0, 1])
@pytest.mark.parametrize("sign", [-1, 1])
@pytest.mark.parametrize("coordinate", [1.4, 2.1])
def test_each_edge_near_and_outside(caplog, axis, sign, coordinate):
position = [0.0, 0.0]
position[axis] = sign * coordinate
make_analysis(positions=[position])
assert len(records(caplog)) == 1
assert records(caplog)[0].levelno == logging.WARNING
assert "does not guarantee image completeness" in records(caplog)[0].message


def test_offset_asymmetric_bounds(caplog):
make_analysis(positions=[(10, -4)], bounds=(8, 13, -7, -2))
assert records(caplog) == []
make_analysis(positions=[(12.5, -4)], bounds=(8, 13, -7, -2))
assert len(records(caplog)) == 1


@pytest.mark.parametrize("scale,noise", [(0.5, 0.1), (0.2, 0.3)])
def test_margin_accounts_for_scale_and_noise(caplog, scale, noise):
make_analysis(positions=[(0.8, 0)], scale=scale, noise=noise)
assert records(caplog)[0].levelno == logging.WARNING


def test_noise_is_paired_with_each_position(caplog):
make_analysis(positions=[(1.5, 0), (0, 0)], noise=[0.01, 0.4])
assert records(caplog) == []


@pytest.mark.parametrize(
"bounds,expected",
[
((-2.1, 2.1, -2.1, 2.1), False),
((-2.2, 2.2, -2.2, 2.2), True),
((-10, 10, -2, 2), False),
((8, 12, -20, 20), False),
],
)
def test_oversize_requires_both_axes_and_margin(caplog, bounds, expected):
positions = [(10, 0)] if bounds[0] > 0 else [(0, 0)]
make_analysis(positions=positions, bounds=bounds)
assert bool(records(caplog)) == expected
if expected:
assert records(caplog)[0].levelno == logging.INFO


@pytest.mark.parametrize(
"variable,value",
[
("PYAUTO_TEST_MODE", "1"),
("PYAUTO_TEST_MODE", "2"),
("PYAUTO_SMALL_DATASETS", "1"),
],
)
def test_modes_suppress_hint_but_keep_warning(monkeypatch, caplog, variable, value):
monkeypatch.setenv(variable, value)
make_analysis(bounds=(-10, 10, -10, 10))
assert records(caplog) == []
make_analysis(positions=[(2.1, 0)])
assert records(caplog)[0].levelno == logging.WARNING


@pytest.mark.parametrize(
"positions,noise",
[(None, None), ([], []), ([[np.nan, 0]], [0.1]), ([[0, 0]], None)],
)
def test_missing_or_unusable_positions(caplog, positions, noise):
solver = make_analysis().solver
dataset = SimpleNamespace(positions=positions, positions_noise_map=noise)
al.AnalysisPoint(dataset=dataset, solver=solver, use_jax=False)
assert records(caplog) == []


@pytest.mark.parametrize("solver", [None, SimpleNamespace()])
def test_missing_or_custom_solver(caplog, solver):
dataset = make_analysis().dataset
al.AnalysisPoint(dataset=dataset, solver=solver, use_jax=False)
assert records(caplog) == []


def test_each_analysis_checks_shared_solver_once_not_per_likelihood(
caplog, monkeypatch
):
analysis = make_analysis(positions=[(2.1, 0)])
other = al.AnalysisPoint(
dataset=analysis.dataset, solver=analysis.solver, use_jax=False
)
for item in (analysis, other):
monkeypatch.setattr(
item, "fit_from", lambda instance: SimpleNamespace(log_likelihood=1.0)
)
for _ in range(3):
assert item.log_likelihood_function(None) == 1.0
assert len(records(caplog)) == 2
Loading