Skip to content

Commit 84d859f

Browse files
Jammy2211claude
authored and
Jammy2211
committed
feat: over/under-prediction policies for point-source pairing likelihoods
FitPositionsImagePairRepeat (model-fit default) gains an unmatched_model_policy class attribute: 'magnification_filter' (default — extra model images below magnification_threshold=0.1 are exempt, the demagnified-central convention; brighter extras add distance-to-nearest-observed residuals normalized by mean noise), 'penalize' (no exemption), 'ignore' (historical behaviour, explicit). Under-prediction gains a finite no_image_residual floor (1e4) when the solver returns no images. n_unmatched_model_positions diagnostic added. residual_map vectorized, fixed-shape/NaN-mask JAX-safe. FitPositionsImagePair (Hungarian): unmatched observed positions (n_model < n_obs) now contribute distance-to-nearest-model residuals in data order — previously silently dropped, rewarding under-prediction (removes the long- standing all-caps do-not-use warning). 5 new regression tests (exact penalty values, policy switching, demagnified exemption via IsothermalSph, no-image floors, Hungarian anti-reward). Full suite 358 pass. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PUuWXiS23FvmfQPLvMNjeM
1 parent 645cc8f commit 84d859f

5 files changed

Lines changed: 485 additions & 160 deletions

File tree

‎autolens/point/fit/positions/image/abstract.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@
1515
one observed position, useful for highly magnified systems where some images may be
1616
too close to separate.
1717
"""
18+
1819
from abc import ABC
1920
import numpy as np
2021
from typing import Optional

‎autolens/point/fit/positions/image/pair.py‎

Lines changed: 108 additions & 94 deletions
Original file line numberDiff line numberDiff line change
@@ -1,94 +1,108 @@
1-
import numpy as np
2-
from scipy.optimize import linear_sum_assignment
3-
4-
import autoarray as aa
5-
6-
from autolens.point.fit.positions.image.abstract import AbstractFitPositionsImagePair
7-
8-
9-
class FitPositionsImagePair(AbstractFitPositionsImagePair):
10-
"""
11-
Fits the positions of a point source dataset using a `Tracer` object with an image-plane chi-squared where every
12-
model position of the point-source is paired with its closest observed position, without allowing for repeated
13-
pairings of the same observed position to model positions.
14-
15-
By not allowing for repeated pairings, this can produce behaviour such as a model position not being paired to
16-
its closest observed position, but instead being paired to a further observed position, if doing so
17-
means that the overall distances of pairings are reduced.
18-
19-
THIS FIT CURRENTLY GIVES UNRELIABLE RESULTS, BECAUSE IT GOES TO SOLUTIONS WHERE THE NUMBER OF MODEL POSITIONS
20-
IS BELOW THE NUMBER OF DATA POSITIONS, REDUCING THE CHI-SQUARED TO LOW VALUES. PYAUTOLENS SHOULD BE UPDATED TO
21-
PENALIZE THIS BEHAVIOUR BEFORE THIS FIT CAN BE USED. THIS REISDUAL MAP PROPERTY MAY ALSO NEED TO BE EXTENDED
22-
TO ACCOUNT FOR NOISE.
23-
24-
The fit performs the following steps:
25-
26-
1) Determine the source-plane centre of the point source, which could be a free model parameter or computed
27-
as the barycenter of ray-traced positions in the source-plane, using name pairing (see below).
28-
29-
2) Determine the image-plane model positions using the `PointSolver` and the source-plane centre of the point
30-
source (e.g. ray tracing triangles to and from the image and source planes), including accounting for
31-
multi-plane ray-tracing.
32-
33-
3) Pair each model position with the observed position, not allowing for repeated pairings of the same
34-
observed position to model positions, to compute the `residual_map`. This may result in some observed
35-
positions not being paired to their closest model position, if doing so reduces the overall distances of
36-
pairings.
37-
38-
5) Compute the chi-squared of each position as the square of the residual divided by the RMS noise-map value.
39-
40-
6) Sum the chi-squared values to compute the overall log likelihood of the fit.
41-
42-
Point source fitting uses name pairing, whereby the `name` of the `Point` object is paired to the name of the
43-
point source dataset to ensure that point source datasets are fitted to the correct point source.
44-
45-
This fit object is used in the `FitPointDataset` to perform position based fitting of a `PointDataset`,
46-
which may also fit other components of the point dataset like fluxes or time delays.
47-
48-
When performing a `model-fit`via an `AnalysisPoint` object the `figure_of_merit` of this object
49-
is called and returned in the `log_likelihood_function`.
50-
51-
Parameters
52-
----------
53-
name
54-
The name of the point source dataset which is paired to a `Point` profile.
55-
data
56-
The positions of the point source in the image-plane which are fitted.
57-
noise_map
58-
The noise-map of the positions which are used to compute the log likelihood of the positions.
59-
tracer
60-
The tracer of galaxies whose point source profile are used to fit the positions.
61-
solver
62-
Solves the lens equation in order to determine the image-plane positions of a point source by ray-tracing
63-
triangles to and from the source-plane.
64-
profile
65-
Manually input the profile of the point source, which is used instead of the one extracted from the
66-
tracer via name pairing if that profile is not found.
67-
"""
68-
69-
@property
70-
def residual_map(self) -> aa.ArrayIrregular:
71-
residual_map = []
72-
73-
cost_matrix = np.linalg.norm(
74-
np.array(
75-
self.data,
76-
)[:, np.newaxis]
77-
- np.array(
78-
self.model_data.array,
79-
),
80-
axis=2,
81-
)
82-
83-
data_indexes, model_indexes = linear_sum_assignment(cost_matrix)
84-
85-
for data_index, model_index in zip(data_indexes, model_indexes):
86-
distance = np.sqrt(
87-
self.square_distance(
88-
self.data[data_index], self.model_data.array[model_index]
89-
)
90-
)
91-
92-
residual_map.append(distance)
93-
94-
return aa.ArrayIrregular(values=residual_map)
1+
import numpy as np
2+
from scipy.optimize import linear_sum_assignment
3+
4+
import autoarray as aa
5+
6+
from autolens.point.fit.positions.image.abstract import AbstractFitPositionsImagePair
7+
8+
9+
class FitPositionsImagePair(AbstractFitPositionsImagePair):
10+
"""
11+
Fits the positions of a point source dataset using a `Tracer` object with an image-plane chi-squared where every
12+
model position of the point-source is paired with its closest observed position, without allowing for repeated
13+
pairings of the same observed position to model positions.
14+
15+
By not allowing for repeated pairings, this can produce behaviour such as a model position not being paired to
16+
its closest observed position, but instead being paired to a further observed position, if doing so
17+
means that the overall distances of pairings are reduced.
18+
19+
**Under-prediction penalty**: the Hungarian assignment pairs ``min(n_observed, n_model)`` positions, which
20+
historically meant that a model predicting *fewer* images than observed silently dropped the unmatched
21+
observed positions from the chi-squared — rewarding under-prediction (samplers drove n_model below
22+
n_observed to shrink the chi-squared). Unmatched observed positions now contribute their distance to the
23+
nearest model position as a residual, and if the model predicts no images at all every observed position
24+
contributes the ``no_image_residual`` floor. ``FitPositionsImagePairRepeat`` remains the model-fit default;
25+
it additionally offers over-prediction policies.
26+
27+
The fit performs the following steps:
28+
29+
1) Determine the source-plane centre of the point source, which could be a free model parameter or computed
30+
as the barycenter of ray-traced positions in the source-plane, using name pairing (see below).
31+
32+
2) Determine the image-plane model positions using the `PointSolver` and the source-plane centre of the point
33+
source (e.g. ray tracing triangles to and from the image and source planes), including accounting for
34+
multi-plane ray-tracing.
35+
36+
3) Pair each model position with the observed position, not allowing for repeated pairings of the same
37+
observed position to model positions, to compute the `residual_map`. This may result in some observed
38+
positions not being paired to their closest model position, if doing so reduces the overall distances of
39+
pairings.
40+
41+
5) Compute the chi-squared of each position as the square of the residual divided by the RMS noise-map value.
42+
43+
6) Sum the chi-squared values to compute the overall log likelihood of the fit.
44+
45+
Point source fitting uses name pairing, whereby the `name` of the `Point` object is paired to the name of the
46+
point source dataset to ensure that point source datasets are fitted to the correct point source.
47+
48+
This fit object is used in the `FitPointDataset` to perform position based fitting of a `PointDataset`,
49+
which may also fit other components of the point dataset like fluxes or time delays.
50+
51+
When performing a `model-fit`via an `AnalysisPoint` object the `figure_of_merit` of this object
52+
is called and returned in the `log_likelihood_function`.
53+
54+
Parameters
55+
----------
56+
name
57+
The name of the point source dataset which is paired to a `Point` profile.
58+
data
59+
The positions of the point source in the image-plane which are fitted.
60+
noise_map
61+
The noise-map of the positions which are used to compute the log likelihood of the positions.
62+
tracer
63+
The tracer of galaxies whose point source profile are used to fit the positions.
64+
solver
65+
Solves the lens equation in order to determine the image-plane positions of a point source by ray-tracing
66+
triangles to and from the source-plane.
67+
profile
68+
Manually input the profile of the point source, which is used instead of the one extracted from the
69+
tracer via name pairing if that profile is not found.
70+
"""
71+
72+
no_image_residual = 1.0e4
73+
74+
@property
75+
def residual_map(self) -> aa.ArrayIrregular:
76+
77+
model = np.asarray(self.model_data.array)
78+
model = model[np.isfinite(model).all(axis=1)]
79+
80+
if model.shape[0] == 0:
81+
return aa.ArrayIrregular(
82+
values=[float(self.no_image_residual)] * len(self.data)
83+
)
84+
85+
cost_matrix = np.linalg.norm(
86+
np.array(
87+
self.data,
88+
)[:, np.newaxis]
89+
- model,
90+
axis=2,
91+
)
92+
93+
data_indexes, model_indexes = linear_sum_assignment(cost_matrix)
94+
95+
# Residuals are ordered by observed position (matching the noise-map ordering): assigned
96+
# positions get their Hungarian pairing distance; positions the assignment could not pair
97+
# (under-prediction, n_model < n_observed) get their distance to the nearest model
98+
# position, so a model that cannot produce an observed image is penalized, not rewarded.
99+
residuals = np.min(cost_matrix, axis=1)
100+
101+
for data_index, model_index in zip(data_indexes, model_indexes):
102+
residuals[data_index] = np.sqrt(
103+
self.square_distance(self.data[data_index], model[model_index])
104+
)
105+
106+
residual_map = [float(r) for r in residuals]
107+
108+
return aa.ArrayIrregular(values=residual_map)

0 commit comments

Comments
 (0)