@@ -146,3 +146,144 @@ def chi_squared(self) -> float:
146146 return ag .util .fit .chi_squared_from (
147147 chi_squared_map = self .chi_squared_map .array ,
148148 )
149+
150+
151+ class FitFluxesSolved (AbstractFitPoint ):
152+ """
153+ Fits the fluxes of a point source dataset with the source-plane flux solved for analytically (in flux space,
154+ magnification-first), following Lombardi 2024 (arXiv:2406.15280) §6.1, rather than read from a free `flux`
155+ model parameter.
156+
157+ With image-plane magnifications `µᵢ` (`magnifications_at_positions`), observed fluxes `f̂ᵢ` and noise `σᵢ`:
158+
159+ `F* = (Σᵢ µᵢ f̂ᵢ/σᵢ²) / (Σᵢ µᵢ²/σᵢ²)` (`solved_flux`)
160+
161+ with model fluxes `µᵢF*` (`model_data`), a standard chi-squared and noise normalization, and the likelihood
162+ analytically marginalized over `F*` (flat prior):
163+
164+ `log_likelihood = -0.5*(χ² + noise_norm) - 0.5*log((Σᵢ µᵢ²/σᵢ²)/(2π))`
165+
166+ The paper's magnitude-space form is not used here: the flux noise maps in this fit are flux-space Gaussians,
167+ and converting to magnitude space would change the error model, not just its parametrization.
168+
169+ Works with any profile that has **no** `flux` attribute (`ag.ps.Point` or `ag.ps.PointSolved`); a profile
170+ with a `flux` attribute (`ag.ps.PointFlux`) raises, since its flux prior would otherwise be sampled by the
171+ non-linear search but silently ignored by the analytic solve. Use `FitFluxes` for a free-flux fit.
172+ """
173+
174+ def __init__ (
175+ self ,
176+ name : str ,
177+ data : aa .ArrayIrregular ,
178+ noise_map : aa .ArrayIrregular ,
179+ positions : aa .Grid2DIrregular ,
180+ tracer : Tracer ,
181+ profile : Optional [ag .ps .Point ] = None ,
182+ xp = np ,
183+ ):
184+ """
185+ Parameters
186+ ----------
187+ name
188+ The name of the point source dataset which is paired to a `Point` profile.
189+ data
190+ The observed fluxes of the point source.
191+ noise_map
192+ The noise-map of the fluxes which are used to compute the log likelihood.
193+ positions
194+ The image-plane positions of the point source where the fluxes and magnifications are calculated.
195+ tracer
196+ The tracer of galaxies whose point source profile is used to fit the fluxes.
197+ profile
198+ Manually input the profile of the point source, used instead of one extracted from the tracer.
199+ """
200+ self .positions = positions
201+
202+ super ().__init__ (
203+ name = name ,
204+ data = data ,
205+ noise_map = noise_map ,
206+ tracer = tracer ,
207+ solver = None ,
208+ profile = profile ,
209+ xp = xp ,
210+ )
211+
212+ if hasattr (self .profile , "flux" ):
213+ raise exc .PointExtractionException (
214+ f"For the point-source named { name }
215+ f"{ self .profile .__class__ .__name__ }
216+ f"for the source flux analytically (F*), so a free `flux` prior would be sampled by the "
217+ f"non-linear search but silently ignored. Use `FitFluxes` with `ag.ps.PointFlux` for a "
218+ f"free-flux fit, or use a profile with no `flux` attribute (e.g. `ag.ps.Point` / "
219+ f"`ag.ps.PointSolved`) with `FitFluxesSolved`."
220+ )
221+
222+ @property
223+ def flux_precision_sum (self ) -> float :
224+ """
225+ `Σᵢ µᵢ²/σᵢ²` — the precision of the solved flux `F*`, and the marginalization normalization.
226+ """
227+ mu = self .magnifications_at_positions .array
228+ sigma_squared = self .noise_map .array ** 2.0
229+ return self ._xp .sum (mu ** 2.0 / sigma_squared )
230+
231+ @property
232+ def solved_flux (self ) -> float :
233+ """
234+ `F* = (Σᵢ µᵢ f̂ᵢ/σᵢ²) / (Σᵢ µᵢ²/σᵢ²)`.
235+ """
236+ mu = self .magnifications_at_positions .array
237+ f_hat = self .data .array
238+ sigma_squared = self .noise_map .array ** 2.0
239+ numerator = self ._xp .sum (mu * f_hat / sigma_squared )
240+ return numerator / self .flux_precision_sum
241+
242+ @property
243+ def model_data (self ) -> aa .ArrayIrregular :
244+ """
245+ The model fluxes `µᵢF*`.
246+ """
247+ return aa .ArrayIrregular (
248+ values = self .magnifications_at_positions .array * self .solved_flux
249+ )
250+
251+ @property
252+ def model_fluxes (self ) -> aa .ArrayIrregular :
253+ return self .model_data
254+
255+ @property
256+ def residual_map (self ) -> aa .ArrayIrregular :
257+ """
258+ Returns the difference between the observed and model fluxes of the point source.
259+ """
260+ residual_map = super ().residual_map
261+
262+ return aa .ArrayIrregular (values = residual_map )
263+
264+ @property
265+ def chi_squared (self ) -> float :
266+ """
267+ Returns the chi-squared of the fit of the point source fluxes.
268+ """
269+ return ag .util .fit .chi_squared_from (
270+ chi_squared_map = self .chi_squared_map .array ,
271+ )
272+
273+ @property
274+ def marginalization_term (self ) -> float :
275+ """
276+ The analytic-marginalization contribution to the log likelihood from integrating out the (flat-prior)
277+ source flux: `-0.5 * log((Σᵢ µᵢ²/σᵢ²)/(2π))`.
278+ """
279+ return - 0.5 * self ._xp .log (self .flux_precision_sum / (2.0 * np .pi ))
280+
281+ @property
282+ def log_likelihood (self ) -> float :
283+ """
284+ `log_likelihood = -0.5*(χ² + noise_norm) - 0.5*log((Σᵢ µᵢ²/σᵢ²)/(2π))`.
285+ """
286+ return (
287+ - 0.5 * (self .chi_squared + self .noise_normalization )
288+ + self .marginalization_term
289+ )
0 commit comments