diff --git a/docs/conf.py b/docs/conf.py index 3a9b2bd1..c62812ac 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -107,7 +107,13 @@ } # Mock imports for modules that may not be available during doc build +# These packages can only be installed from GitHub and are not available +# on ReadTheDocs or in the standard docs build environment: +# ne2001: pip install git+https://github.com/FRBs/ne2001.git +# frb: pip install git+https://github.com/FRBs/FRB.git +# astropath: pip install git+https://github.com/FRBs/astropath.git autodoc_mock_imports = [ 'ne2001', 'frb', + 'astropath', ] diff --git a/docs/index.rst b/docs/index.rst index 08c536ad..6bebe91a 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -43,6 +43,7 @@ Getting Started architecture parameters api + optical .. toctree:: :maxdepth: 1 diff --git a/docs/optical.rst b/docs/optical.rst new file mode 100644 index 00000000..808897ae --- /dev/null +++ b/docs/optical.rst @@ -0,0 +1,435 @@ +.. _optical: + +================================ +Optical Host Galaxy Association +================================ + +This section describes the three modules that connect zdm's redshift-DM +predictions to the `PATH `_ (Probabilistic +Association of Transients to their Hosts) algorithm. Together they provide +physically motivated priors on FRB host galaxy apparent magnitudes, derived +from the zdm posterior p(z | DM\ :sub:`EG`). + +- :mod:`zdm.optical_params` — parameter dataclasses configuring each model +- :mod:`zdm.optical` — host magnitude models and the PATH interface wrapper +- :mod:`zdm.optical_numerics` — numerical evaluation, optimisation, and + statistics for fitting the models to CRAFT ICS optical data + +Overview +======== + +Standard PATH assigns host galaxy candidates a prior based only on galaxy +surface density and angular size. The zdm optical modules replace this with a +prior informed by p(z | DM\ :sub:`EG`): given an FRB's extragalactic DM, zdm +predicts a redshift distribution, which is convolved with a host galaxy +luminosity model to produce p(m\ :sub:`r` | DM\ :sub:`EG`). + +The modules are built around a two-layer design: + +.. code-block:: text + + ┌──────────────────────────────────────────────────────────────┐ + │ model_wrapper (optical.py) │ + │ Convolves p(m_r|z) with zdm p(z|DM_EG) → p(m_r|DM_EG) │ + │ Estimates P_U (undetected host prior) │ + │ Plugs into PATH via pathpriors.USR_raw_prior_Oi │ + └──────────────────────────────────────────────────────────────┘ + │ wraps + ┌──────────────────────────────────────────────────────────────┐ + │ Host magnitude models (optical.py) │ + │ │ + │ simple_host_model — parametric p(M_r) histogram │ + │ loudas_model — mass/SFR-weighted tables (Loudas) │ + │ marnoch_model — Gaussian fit to known hosts │ + │ (Marnoch et al. 2023) │ + └──────────────────────────────────────────────────────────────┘ + │ configured by + ┌──────────────────────────────────────────────────────────────┐ + │ Parameter dataclasses (optical_params.py) │ + │ │ + │ OpticalState ← SimpleParams, LoudasParams, │ + │ Apparent, Identification │ + └──────────────────────────────────────────────────────────────┘ + +Host Magnitude Models +===================== + +Three models are available, all implementing the same interface +``get_pmr_gz(mrbins, z)`` which returns p(m\ :sub:`r` | z) for a set of +apparent magnitude bin edges at a given redshift. + +simple_host_model +----------------- + +A parametric model describing the intrinsic absolute magnitude distribution +p(M\ :sub:`r`) as N amplitudes (default 10) at uniformly spaced points +between ``Absmin`` and ``Absmax``. The amplitudes are normalised to sum to +unity and interpolated onto a fine internal grid via one of four schemes +controlled by ``AbsModelID``: + +.. list-table:: + :header-rows: 1 + :widths: 15 85 + + * - ``AbsModelID`` + - Description + * - 0 + - Step-function histogram — each parameter value applies uniformly to + its bin + * - 1 + - Linear interpolation between parameter points *(default)* + * - 2 + - Cubic spline interpolation (negative values clamped to zero) + * - 3 + - Cubic spline in log-space (parameters are log\ :sub:`10` weights) + +Conversion from M\ :sub:`r` to m\ :sub:`r` is controlled by ``AppModelID``: + +.. list-table:: + :header-rows: 1 + :widths: 15 85 + + * - ``AppModelID`` + - Description + * - 0 + - Pure distance modulus, no k-correction *(default)* + * - 1 + - Distance modulus plus power-law k-correction + ``2.5 × k × log10(1 + z)`` + +loudas_model +------------ + +Uses precomputed p(m\ :sub:`r` | z) tables from Nick Loudas, constructed by +weighting galaxy luminosities by either stellar mass or star-formation rate. +The single free parameter ``fSFR`` interpolates between the two: + +.. math:: + + p(m_r | z) = (1 - f_{\rm SFR})\,p_{\rm mass}(m_r | z) + + f_{\rm SFR}\,p_{\rm SFR}(m_r | z) + +Interpolation between tabulated redshift bins is performed in log-z space +with a luminosity-distance shift applied to each tabulated distribution before +combining, ensuring correct apparent magnitude evolution at low redshift. + +marnoch_model +------------- + +A zero-parameter model based on Marnoch et al. 2023 (MNRAS 525, 994). Fits +a Gaussian to the r-band magnitude distribution of known CRAFT ICS FRB host +galaxies, with mean and standard deviation described as cubic splines of +redshift. No free parameters; the model is fixed by the observed host sample. + +The ``model_wrapper`` Class +=========================== + +:class:`~zdm.optical.model_wrapper` is a survey-independent wrapper around +any host model. Its key responsibilities are: + +1. **Precomputation**: at initialisation it calls ``model.get_pmr_gz`` for + every redshift value in the zdm grid to build a cached + p(m\ :sub:`r` | z) array. +2. **DM integration**: ``init_path_raw_prior_Oi(DM, grid)`` extracts + p(z | DM\ :sub:`EG`) from the grid and convolves it with the cached + array to produce p(m\ :sub:`r` | DM\ :sub:`EG`). +3. **P_U estimation**: ``estimate_unseen_prior()`` integrates the magnitude + prior against the detection probability curve + (logistic function centred on ``pU_mean`` with width ``pU_width``) to + obtain the prior probability that the true host is below the detection + limit. +4. **PATH interface**: after ``init_path_raw_prior_Oi`` is called, + ``pathpriors.USR_raw_prior_Oi`` is automatically pointed at + ``path_raw_prior_Oi``, so PATH uses the zdm-derived prior transparently. + +Typical Workflow +================ + +The following example shows how to obtain zdm-informed PATH posteriors for a +single CRAFT ICS FRB. + +.. code-block:: python + + from zdm import optical as opt + from zdm import optical_numerics as on + from zdm import loading, cosmology as cos, parameters + + # 1. Initialise zdm grid + state = parameters.State() + cos.set_cosmology(state) + cos.init_dist_measures() + ss, gs = loading.surveys_and_grids(survey_names=['CRAFT_ICS_1300']) + g, s = gs[0], ss[0] + + # 2. Choose a host magnitude model + model = opt.marnoch_model() # or simple_host_model / loudas_model + + # 3. Wrap it for the survey's redshift grid + wrapper = opt.model_wrapper(model, g.zvals) + + # 4. For a specific FRB, look up its DM_EG + frb = 'FRB20190608B' + imatch = opt.matchFRB(frb, s) + DMEG = s.DMEGs[imatch] + + # 5. Compute p(m_r | DM_EG) and estimate P_U + wrapper.init_path_raw_prior_Oi(DMEG, g) # also sets pathpriors.USR_raw_prior_Oi + PU = wrapper.estimate_unseen_prior() + + # 6. Run PATH with the zdm prior + P_O, P_Ox, P_Ux, mags, ptbl = on.run_path(frb, usemodel=True, P_U=PU) + +To process the full CRAFT ICS sample and compare models, use +:func:`~zdm.optical_numerics.calc_path_priors` directly. To fit model +parameters, pass :func:`~zdm.optical_numerics.function` as the objective to +``scipy.optimize.minimize`` — see ``zdm/scripts/Path/optimise_host_priors.py`` +for a complete example. + +Parameter Reference +=================== + +All host galaxy model parameters are held in dataclasses collected by +:class:`~zdm.optical_params.OpticalState`. The four constituent dataclasses +and their parameters are described below. + +SimpleParams +------------ + +Controls the :class:`~zdm.optical.simple_host_model`. + +.. list-table:: + :header-rows: 1 + :widths: 20 12 12 56 + + * - Parameter + - Default + - Units + - Description + * - ``Absmin`` + - −25 + - M\ :sub:`r` + - Minimum absolute magnitude of the host distribution + * - ``Absmax`` + - −15 + - M\ :sub:`r` + - Maximum absolute magnitude of the host distribution + * - ``NAbsBins`` + - 1000 + - — + - Number of internal absolute magnitude bins (fine grid for + computing p(m\ :sub:`r` | z)) + * - ``NModelBins`` + - 10 + - — + - Number of free parameter bins describing p(M\ :sub:`r`) + * - ``AbsPriorMeth`` + - 0 + - — + - Initial prior on absolute magnitudes: 0 = uniform + * - ``AbsModelID`` + - 1 + - — + - Interpolation scheme for p(M\ :sub:`r`): 0 = histogram, + 1 = linear, 2 = spline, 3 = log-spline + * - ``AppModelID`` + - 0 + - — + - Absolute-to-apparent conversion: 0 = distance modulus only, + 1 = with power-law k-correction + * - ``k`` + - 0.0 + - — + - k-correction power-law index (only used when ``AppModelID=1``) + +LoudasParams +------------ + +Controls the :class:`~zdm.optical.loudas_model`. + +.. list-table:: + :header-rows: 1 + :widths: 20 12 12 56 + + * - Parameter + - Default + - Units + - Description + * - ``fSFR`` + - 0.5 + - — + - Fraction of FRB hosts tracing star-formation rate (0 = pure + mass-weighted, 1 = pure SFR-weighted) + * - ``NzBins`` + - 10 + - — + - Number of redshift bins for histogram calculations + * - ``zmin`` + - 0.0 + - — + - Minimum redshift for p(m\ :sub:`r`) calculation + * - ``zmax`` + - 0.0 + - — + - Maximum redshift for p(m\ :sub:`r`) calculation + * - ``NMrBins`` + - 0 + - — + - Number of absolute magnitude bins + * - ``Mrmin`` + - 0.0 + - M\ :sub:`r` + - Minimum absolute magnitude + * - ``Mrmax`` + - 0.0 + - M\ :sub:`r` + - Maximum absolute magnitude + +Apparent +-------- + +Controls the apparent magnitude grid used by :class:`~zdm.optical.model_wrapper`. + +.. list-table:: + :header-rows: 1 + :widths: 20 12 12 56 + + * - Parameter + - Default + - Units + - Description + * - ``Appmin`` + - 10 + - m\ :sub:`r` + - Minimum apparent magnitude of the internal grid + * - ``Appmax`` + - 35 + - m\ :sub:`r` + - Maximum apparent magnitude of the internal grid + * - ``NAppBins`` + - 250 + - — + - Number of apparent magnitude bins + +Identification +-------------- + +Controls the survey detection completeness model used to compute P_U. +The detection probability is modelled as a logistic function: +p(detected | m\ :sub:`r`) = 1 − p(U | m\ :sub:`r`) where + +.. math:: + + p(U | m_r) = \frac{1}{1 + \exp\!\left(\frac{\mu - m_r}{w}\right)} + +with μ = ``pU_mean`` and w = ``pU_width``. + +.. list-table:: + :header-rows: 1 + :widths: 20 12 12 56 + + * - Parameter + - Default + - Units + - Description + * - ``pU_mean`` + - 26.385 + - m\ :sub:`r` + - Magnitude at which 50 % of host galaxies are undetected + (the survey's half-completeness limit). Default value is + calibrated to VLT/FORS2 R-band observations. + * - ``pU_width`` + - 0.279 + - m\ :sub:`r` + - Characteristic width of the completeness rolloff. Smaller + values give a sharper transition between detected and + undetected regimes. + +Optimisation and Statistics +============================ + +:mod:`zdm.optical_numerics` provides two goodness-of-fit statistics for +comparing the model-predicted apparent magnitude prior to observed PATH +posteriors across a sample of FRBs: + +**Maximum-likelihood statistic** (:func:`~zdm.optical_numerics.calculate_likelihood_statistic`) + +For each FRB, evaluates + +.. math:: + + \ln \mathcal{L}_i = \log_{10}\!\left(\sum_j \frac{P(O_j|x)}{s_i} + P_{U,i}^{\rm prior}\right) + +where the sum runs over candidate host galaxies and *s*\ :sub:`i` is a +rescale factor that undoes PATH's internal renormalisation. The total +statistic is Σ ln *ℒ*\ :sub:`i` over all FRBs. This is the recommended +objective for parameter fitting. + +**KS-like statistic** (:func:`~zdm.optical_numerics.calculate_ks_statistic`) + +Builds normalised cumulative distributions of the model prior and the +observed posteriors (weighted by P(O|x)) over apparent magnitude, then +returns the maximum absolute difference — analogous to the KS test statistic. +Smaller values indicate a better fit. + +Both statistics accept a ``POxcut`` argument to restrict the sample to FRBs +with a confidently identified host (max P(O|x) > threshold), simulating a +traditional host-identification approach. + +Scripts +======= + +Ready-to-run scripts using these modules are in ``zdm/scripts/Path/``: + +.. list-table:: + :header-rows: 1 + :widths: 40 60 + + * - Script + - Purpose + * - ``estimate_path_priors.py`` + - Demonstrate zdm-informed PATH priors on all CRAFT ICS FRBs; + compare flat vs. model priors; save posterior magnitude histogram + * - ``optimise_host_priors.py`` + - Fit host model parameters to the CRAFT ICS sample using + ``scipy.optimize.minimize`` + * - ``plot_host_models.py`` + - Visualise all three host models and compare their PATH posteriors + across the CRAFT ICS sample + +API Reference +============= + +optical_params +-------------- + +Parameter dataclasses for configuring host galaxy models. + +.. automodapi:: zdm.optical_params + :no-inheritance-diagram: + +optical +------- + +Host magnitude model classes and the ``model_wrapper`` PATH interface. + +.. automodapi:: zdm.optical + :no-inheritance-diagram: + +optical_numerics +---------------- + +Numerical evaluation, optimisation, and statistics for fitting host models. + +.. automodapi:: zdm.optical_numerics + :no-inheritance-diagram: + +References +========== + +- Marnoch et al. 2023, MNRAS 525, 994 — + FRB host galaxy r-band magnitude model + (https://doi.org/10.1093/mnras/stad2353) +- Macquart et al. 2020, Nature 581, 391 — + Macquart relation / p(DM | z) +- Aggarwal et al. 2021, ApJ 911, 95 — + PATH algorithm for probabilistic host association diff --git a/zdm/grid.py b/zdm/grid.py index be1ef34e..dea68aac 100644 --- a/zdm/grid.py +++ b/zdm/grid.py @@ -1260,7 +1260,45 @@ def chk_upd_param(self, param: str, vparams: dict, update=False): def smear_z(self,array,zsigma): """ - smears the z-grid according to a specified photometric error + Smear a 2-D z-DM grid along the redshift axis to account for + photometric redshift uncertainty. + + When a survey uses photometric rather than spectroscopic redshifts, + the true redshift of each FRB host is uncertain. This method convolves + each column of the grid (i.e. each fixed-DM slice along the z axis) + with a Gaussian kernel whose standard deviation equals ``zsigma``, + redistributing probability across neighbouring redshift bins. + + The kernel is truncated at ``state.photo.sigma_width`` standard + deviations on each side (default 6σ), rounded up to an odd number of + bins so that it is centred exactly on zero. + + Parameters + ---------- + array : np.ndarray, shape (Nz, NDM) + Input 2-D grid with redshift along axis 0 and DM along axis 1. + zsigma : float + Photometric redshift uncertainty (1σ), in the same units as + ``self.zvals`` (i.e. dimensionless redshift). + + Returns + ------- + smear_zgrid : np.ndarray, shape (Nz, NDM) + Copy of ``array`` with each DM column convolved along the z axis + by the Gaussian smearing kernel. Boundary effects are handled with + ``np.convolve`` mode ``"same"``, so the output has the same shape + as the input. + + Notes + ----- + The kernel width in grid bins is ``zsigma / self.dz``. Values near the + grid edges will be underestimated because the convolution truncates to + zero outside the grid; for well-chosen grid extents this edge effect is + negligible. + + In ``calc_rates``, this method is called with ``zsigma`` taken from + ``self.survey.survey_data.observing.Z_PHOTO`` and applied to + ``self.rates`` after the FRB rate grid has been computed. """ r,c=array.shape # get sigma in grid units diff --git a/zdm/optical.py b/zdm/optical.py index 26ffa022..eaa398c1 100644 --- a/zdm/optical.py +++ b/zdm/optical.py @@ -1,47 +1,58 @@ """ -This library contains routines that interact with -the FRB/astropath module and (optical) FRB host galaxy -information. +Optical FRB host galaxy models and PATH interface for zdm. -The philosophy of the module is this. The base class -is "host_model". This class is the top-level class -that contains base functions to e.g. calculate -p(m_r|DM). +This module connects zdm redshift-DM grids with the PATH +(Probabilistic Association of Transients to their Hosts) algorithm by +providing physically motivated priors on FRB host galaxy apparent +magnitudes, p(m_r | DM_EG). -However, no host_model class contains any astroiphysics. +Architecture +------------ +The module is built around a two-layer design: -Instead, it wraps an underlying set of possible class -objects that must have a specific set of callable functions -which each contain the relevant calculations. +**Host magnitude models** — each describes the intrinsic absolute +magnitude distribution of FRB host galaxies, p(M_r), and can convert +it to an apparent magnitude distribution p(m_r | z) at a given +redshift. Three models are provided: -The current set are the following: +- ``simple_host_model``: parametric histogram of p(M_r) with N free + amplitudes (default 10), interpolated linearly or via spline. An + optional power-law k-correction can be included. N (or N+1) + free parameters. -Simple_host_model: - Describes intrinsic host properties as a spline - interpolation between p(M_r) described by N - points. N parameters (e.g. 10). +- ``loudas_model``: precomputed p(m_r | z) tables from Nick Loudas, + constructed by weighting galaxies by stellar mass or star-formation + rate. Interpolated between tabulated redshift bins with a + luminosity-distance shift. Single free parameter ``fSFR`` sets the + SFR/mass mixing fraction. -Marnoch_model: - Fixed calculation of p(M_r) based on extrapolation - of known FRB host galaxies. No parameters. See - https://doi.org/10.1093/mnras/stad2353 - +- ``marnoch_model``: zero-parameter model. Fits a Gaussian to the + r-band magnitude distribution of known CRAFT ICS host galaxies from + Marnoch et al. 2023 (MNRAS 525, 994), using cubic splines for the + redshift-dependent mean and standard deviation. -Loudas_model: - Calculates p(M_r) via assigning a fraction of FRB - hosts to follow star-formation in galaxies, and - a fraction to stellar mass, then includes the modelled - evolution of these galaxies. 1 parameter. - -Each "host_model" class object above must provide functions to: - __init__ - calculate p(m_r|z,parameters) - -The wrapper class provides the following fubctions: - -- init_path_raw_prior_Oi(self,DM,grid): - (takes as input an FRB DM, and grid object) +**``model_wrapper``** — a survey-independent wrapper around any host +model. Given a zdm ``Grid`` and an observed DM_EG, it convolves +p(m_r | z) with the zdm p(z | DM_EG) posterior to produce a +DM-informed apparent magnitude prior for PATH. It also estimates +P_U, the prior probability that the true host is below the survey +detection limit. + +Typical usage +------------- +:: + + model = opt.marnoch_model() + wrapper = opt.model_wrapper(model, grid.zvals) + wrapper.init_path_raw_prior_Oi(DMEG, grid) # DM-specific initialisation + PU = wrapper.estimate_unseen_prior() + # pathpriors.USR_raw_prior_Oi is now set automatically by init_path_raw_prior_Oi +Module-level data +----------------- +``frblist`` : list of str + TNS names of CRAFT ICS FRBs for which PATH optical data are + available (used by the scripts in ``zdm/scripts/Path/``). """ @@ -109,7 +120,19 @@ def load_data(self): def process_rbands(self): """ - Returns parameters of the host magnitude distribution as a function of redshift + Build cubic spline fits to the mean and rms of p(m_r) as a function of z. + + Reads the per-FRB r-band magnitude columns from ``self.table``, + computes their mean (``Rbar``) and sample standard deviation (``Rrms``) + across all FRBs at each tabulated redshift, then fits two cubic splines: + + - ``self.sbar``: CubicSpline interpolating the mean apparent magnitude + as a function of redshift. + - ``self.srms``: CubicSpline interpolating the rms scatter as a + function of redshift. + + These splines are subsequently used by ``get_pmr_gz`` to evaluate + the Gaussian p(m_r | z) at arbitrary redshifts. """ table = self.table @@ -138,15 +161,32 @@ def process_rbands(self): #return Rbar,Rrms,zlist,sbar,srms - def get_pmr_gz(self,mrbins,z): # fsfr must be a self value z: float,fsfr: float): + def get_pmr_gz(self,mrbins,z): """ - Returns the p_mr distribution for a given redshift z and sfr fraction f_sfr - - Args: - mrbins (array of floats): list of r-band magnitude bins - z (float): redshift + Return the apparent magnitude probability distribution p(m_r | z). + + Evaluates a Gaussian distribution whose mean and standard deviation + are obtained from the cubic splines fit in ``process_rbands``, + and integrates it over the provided magnitude bins. + + This model has no free parameters; the Gaussian moments are fully + determined by the Marnoch et al. 2023 host galaxy data. + + Parameters + ---------- + mrbins : array-like of float, length N+1 + Edges of the apparent magnitude bins over which to compute the + probability. The output has length N (one value per bin). + z : float + Redshift at which to evaluate the magnitude distribution. + + Returns + ------- + pmr : np.ndarray, length N + Probability in each magnitude bin (sums to ≤ 1; may be less + than unity if the Gaussian extends beyond the bin range). """ - + mean = self.sbar(z) rms = self.srms(z) @@ -169,12 +209,16 @@ class loudas_model: def __init__(self,OpticalState=None,fname='p_mr_distributions_dz0.01_z_in_0_1.2.h5',data_dir=None,verbose=False): """ - initialises the model. Loads data provided by Nick Loudas - on mass- and sfr-weighted magnitudes. - + Initialise the Loudas model, loading precomputed p(m_r | z) tables. + Args: - fname [string]: h55 filename containing the data - datadir [string]: directory that the data is contained in. Defaults to None. + OpticalState (OpticalState, optional): optical parameter state. A + default ``OpticalState`` is created if not provided. + fname (str): HDF5 filename containing the Loudas p(m_r | z) tables. + Defaults to ``'p_mr_distributions_dz0.01_z_in_0_1.2.h5'``. + data_dir (str, optional): directory containing ``fname``. Defaults + to the package data directory ``zdm/data/optical/``. + verbose (bool): if True, print progress messages. Defaults to False. """ # uses the "simple hosts" descriptor @@ -246,15 +290,31 @@ def init_cubics(self): self.mass_splines = mass_splines self.sfr_splines = sfr_splines - def get_pmr_gz(self,mrbins,z): # fsfr must be a self value z: float,fsfr: float): + def get_pmr_gz(self,mrbins,z): """ - Returns the p_mr distribution for a given redshift z and sfr fraction f_sfr - Should be defined such that the sum over all mrbins is unity (or less, - if there is a limitation due to range) - - Args: - z (float): redshift - fsfr (float): fraction of population associated with star-formation + Return the apparent magnitude probability distribution p(m_r | z). + + Interpolates between the two nearest tabulated redshift bins (in + log-z space), applying a luminosity-distance shift to each tabulated + p(m_r) before combining them. The mass- and SFR-weighted distributions + are mixed according to ``self.fsfr`` (set via ``init_args``). + + The result is normalised so that the bin probabilities sum to unity + over the full magnitude range, provided the distribution does not + extend significantly beyond ``mrbins``. + + Parameters + ---------- + mrbins : array-like of float, length N+1 + Edges of the apparent magnitude bins. Output has length N. + z : float + Redshift at which to evaluate the distribution. Values outside + the tabulated range are extrapolated from the nearest edge bin. + + Returns + ------- + pmr : np.ndarray, length N + Probability in each apparent magnitude bin (sums to ≤ 1). """ fsfr = self.fsfr @@ -383,15 +443,18 @@ def load_p_mr_distributions(self,data_dir,fname: str = 'p_mr_distributions_dz0.0 def give_p_mr_mass(self,z: float): """ - Function to return p(mr|z) for mass-weighted population. + Return p(m_r | z) for the stellar-mass-weighted host population. + + Uses a nearest-bin lookup (no interpolation) in the tabulated redshift + grid. For interpolated results with luminosity-distance shifting, use + ``get_pmr_gz`` with ``fSFR=0`` instead. + Args: z (float): Redshift value. + Returns: - np.array: p(mr|z) values. - Note: - This function assumes that the redshift bins are defined in the `massweighted_population` data. - Given the fine discretization of redshift bins, it uses the nearest bin for the provided redshift value. - rmag_centers and p_mr_mass are defined in the outer scope of this function. + np.ndarray: p(m_r | z) values at the tabulated r-band magnitude + centres (``self.rmags``), normalised to sum to unity. """ # Find the appropriate redshift bin index idx = np.clip(np.searchsorted(self.zbins, z) - 1, 0, n_redshift_bins - 1) @@ -399,15 +462,18 @@ def give_p_mr_mass(self,z: float): def give_p_mr_sfr(self,z: float): """ - Function to return p(mr|z) for SFR-weighted population. + Return p(m_r | z) for the star-formation-rate-weighted host population. + + Uses a nearest-bin lookup (no interpolation) in the tabulated redshift + grid. For interpolated results with luminosity-distance shifting, use + ``get_pmr_gz`` with ``fSFR=1`` instead. + Args: z (float): Redshift value. + Returns: - np.array: p(mr|z) values. - Note: - This function assumes that the redshift bins are defined in the `sfrweighted_population` data. - Given the fine discretization of redshift bins, it uses the nearest bin for the provided redshift value. - rmag_centers and p_mr_sfr are defined in the outer scope of this function. + np.ndarray: p(m_r | z) values at the tabulated r-band magnitude + centres (``self.rmags``), normalised to sum to unity. """ # Find the appropriate redshift bin index idx = np.clip(np.searchsorted(self.zbins, z) - 1, 0, n_redshift_bins - 1) @@ -415,11 +481,14 @@ def give_p_mr_sfr(self,z: float): def init_args(self,fSFR): """ - Initialises prior based on sfr fraction - + Set the SFR/mass mixing fraction for the Loudas model. + Args: - opstate: optical model state. Grabs the Loudas parameters from there. - + fSFR (float or array-like of length 1): fraction of FRB hosts + that trace star-formation rate. ``fSFR=0`` gives a purely + mass-weighted population; ``fSFR=1`` gives a purely + SFR-weighted population. Intermediate values linearly mix + the two. If an array is passed, only the first element is used. """ # for numerical purposes, fSFR may have to be a vector if hasattr(fSFR,'__len__'): @@ -475,13 +544,15 @@ class simple_host_model: """ def __init__(self,OpticalState=None,verbose=False): """ - Class constructor. - + Initialise the simple host magnitude model. + Args: - opstate (class: Hosts, optional): class defining parameters - of optical state model - verbose (bool, optional): to be verbose y/n - + OpticalState (OpticalState, optional): optical parameter state + providing model configuration (magnitude ranges, number of + bins, interpolation scheme, k-correction flag). A default + ``OpticalState`` is created if not provided. + verbose (bool, optional): if True, print which sub-models are + being initialised. Defaults to False. """ # uses the "simple hosts" descriptor if OpticalState is None: @@ -647,25 +718,35 @@ def init_model_bins(self): def get_pmr_gz(self,mrbins,z): """ - For a set of redshifts, initialise mapping - between intrinsic magnitudes and apparent magnitudes - - This routine only needs to be called once, since the model - to convert absolute to apparent magnitudes is fixed - - It is not set automatically however, and needs to be called - with a set of z values. This is all for speedup purposes. - - Args: - mrbins (np.array, float, length N+1): array of apparent magnitudes (mr) - over which to calculate p(mr). These act as bins - in apparent magnitude mr for histogram purposes, - i.e. they are not probabilities *at* mr - zvals (float): redshifts at which - to map absolute to apparent magnitudes. - - Returns: - pmr: probability for each of the bins (length: N) + Return the apparent magnitude probability distribution p(m_r | z). + + Converts each apparent magnitude bin centre back to an absolute + magnitude using ``CalcAbsoluteMags``, then linearly interpolates + the absolute magnitude weight array (``self.AbsMagWeights``) to + obtain a probability density at each bin. + + Parameters + ---------- + mrbins : np.ndarray of float, length N+1 + Edges of the apparent magnitude bins. The output has length N, + with one probability value per bin centre. + z : float + Redshift at which to evaluate the distribution. + + Returns + ------- + pmr : np.ndarray, length N + Probability density at each apparent magnitude bin centre. + Values at the edges of the absolute magnitude range are + clamped to the nearest valid bin rather than extrapolated. + + Notes + ----- + The returned values are NOT renormalised to sum to unity; the sum + may be less than one if some absolute magnitudes lie outside the + model range ``[Absmin, Absmax]``. This is intentional: the + shortfall represents probability mass for hosts too faint or too + bright to appear in the apparent magnitude range. """ old = False @@ -921,23 +1002,28 @@ def init_zmapping(self,zvals): def init_path_raw_prior_Oi(self,DM,grid): """ - Initialises the priors for a particlar DM. - This performs a function very similar to - "get_posterior" except that it expicitly - only operates on a single DM, and saves the - information internally so that - path_raw_prior_Oi can be called for numerous - host galaxy candidates. - - It returns the priors distribution. - + Initialise the apparent magnitude prior for a single FRB DM. + + Computes p(m_r | DM_EG) by convolving the precomputed p(m_r | z) + grid (``self.p_mr_z``) with the zdm posterior p(z | DM_EG) extracted + from ``grid``. The result is stored internally so that + ``path_raw_prior_Oi`` can be called repeatedly for different host + galaxy candidates belonging to the same FRB without recomputing the + DM integral. + + Also computes the probability that the host is undetected: + - ``self.priors``: p(m_r | DM) weighted by the detection probability + p(detected | m_r). + - ``self.PUdist``: the magnitude-resolved contribution to P_U. + - ``self.PU``: scalar total prior probability that the host is unseen, + returned by ``estimate_unseen_prior()``. + + After this call, ``pathpriors.USR_raw_prior_Oi`` is automatically + pointed at ``self.path_raw_prior_Oi``. + Args: - DM [float]: dispersion measure of an FRB (pc cm-3) - grid (class grid): initialised grid object from which - to calculate priors - - Returns: - priors (float): vector of priors on host galaxy apparent magnitude + DM (float): extragalactic dispersion measure of the FRB (pc cm⁻³). + grid (Grid): initialised zdm grid object providing p(z, DM). """ # we start by getting the posterior distribution p(z) @@ -966,20 +1052,26 @@ def init_path_raw_prior_Oi(self,DM,grid): def get_posterior(self, grid, DM): """ - Similar functionality to init_path_raw_prior_Oi. May be legacy code. - - Returns posterior redshift distributiuon for a given grid, and DM - magnitude distribution, for FRBs of DM given a grid object. - Note: this calculates a prior for PATH, but is a posterior - from zDM's point of view. - + Return apparent magnitude and redshift posteriors for a given DM. + + Computes p(z | DM) from the grid and convolves it with the + precomputed ``self.maghist`` to obtain p(m_r | DM). + + Note: from PATH's perspective this is a prior on host magnitude, + but from zdm's perspective it is a posterior on redshift given DM. + + This method predates ``init_path_raw_prior_Oi`` and may not be + actively used in current scripts. + Args: - grid (class grid object): grid object defining p(z,DM) - DM (float, np.ndarray OR scalar): FRB DM(s) - + grid (Grid): initialised zdm grid object providing p(z, DM). + DM (float or np.ndarray): FRB extragalactic DM(s) in pc cm⁻³. + Returns: - papps (np.ndarray, floats): probability distribution of apparent magnitudes given DM - pz (np.ndarray, floats): probability distribution of redshift given DM + papps (np.ndarray): probability distribution of apparent magnitude + given DM, p(m_r | DM). + pz (np.ndarray): probability distribution of redshift given DM, + p(z | DM). """ # Step 1: get prior on z pz = get_pz_prior(grid,DM) @@ -1033,8 +1125,24 @@ def estimate_unseen_prior(self): def path_base_prior(self,mags): """ - Calculates base magnitude prior. Does NOT include - galaxy density factor + Evaluate the apparent magnitude prior p(m_r | DM) at a list of magnitudes. + + Linearly interpolates ``self.priors`` (which already incorporates the + detection probability p(detected | m_r)) at each requested magnitude, + converting from the internally normalised sum-to-unity convention to a + probability density by dividing by the bin width ``self.dAppmag``. + + Unlike ``path_raw_prior_Oi``, this method does NOT divide by the + galaxy surface density Sigma_m, so it returns the raw magnitude prior + without the PATH normalisation factor. + + Args: + mags (list or tuple of float): apparent r-band magnitudes of + candidate host galaxies at which to evaluate the prior. + + Returns: + Ois (np.ndarray): prior probability density p(m_r | DM) evaluated + at each magnitude in ``mags``. """ ngals = len(mags) Ois = [] @@ -1189,22 +1297,22 @@ def get_pz_prior(grid, DM): def SimplekApparentMags(Abs,k,zs): """ - Function to convert galaxy absolue to apparent magnitudes. - Same as simple apparent mags, but allows for a k-correction. - - Nominally, magnitudes are r-band magnitudes, but this function - is so simple it doesn't matter. - - Just applies a distance correction - no k-correction. - + Convert absolute to apparent magnitudes with a power-law k-correction. + + Applies the distance modulus plus a k-correction of the form + ``2.5 * k * log10(1 + z)``. + Args: - Abs (float or array of floats): intrinsic galaxy luminosities - k (float): k-correction - zs (float or array of floats): redshifts of galaxies - + Abs (float or np.ndarray): absolute magnitude(s) M_r. + k (float): k-correction power-law index. ``k=0`` reduces to a + pure distance modulus (identical to ``SimpleApparentMags``). + zs (float or np.ndarray): redshift(s) of the galaxies. + Returns: - ApparentMags: NAbs x NZ array of magnitudes, where these - are the dimensions of the inputs + ApparentMags: apparent magnitude(s). Scalar if both inputs are + scalar; 1-D array if one is scalar and one is an array; 2-D + array of shape (NAbs, Nz) if both are arrays (computed via + ``np.outer``). """ # calculates luminosity distances (Mpc) @@ -1239,22 +1347,22 @@ def SimplekApparentMags(Abs,k,zs): def SimplekAbsoluteMags(App,k,zs): """ - Function to convert galaxy apparent to absolute magnitudes. - Same as simple absolute mags mags, but allows for a k-correction. - - Nominally, magnitudes are r-band magnitudes, but this function - is so simple it doesn't matter. - - Just applies a distance correction - no k-correction. - + Convert apparent to absolute magnitudes with a power-law k-correction. + + Inverse of ``SimplekApparentMags``: subtracts the distance modulus and + k-correction ``2.5 * k * log10(1 + z)`` from the apparent magnitude. + Args: - App (float or array of floats): apparent galaxy luminosities - k (float): k-correction - zs (float or array of floats): redshifts of galaxies - + App (float or np.ndarray): apparent magnitude(s) m_r. + k (float): k-correction power-law index. ``k=0`` reduces to a + pure distance modulus (identical to ``SimpleAbsoluteMags``). + zs (float or np.ndarray): redshift(s) of the galaxies. + Returns: - AbsoluteMags: NAbs x NZ array of magnitudes, where these - are the dimensions of the inputs + AbsoluteMags: absolute magnitude(s). Scalar if both inputs are + scalar; 1-D array if one is scalar and one is an array; 2-D + array of shape (NApp, Nz) if both are arrays (computed via + ``np.outer``). """ # calculates luminosity distances (Mpc) @@ -1289,20 +1397,20 @@ def SimplekAbsoluteMags(App,k,zs): def SimpleAbsoluteMags(App,zs): """ - Function to convert galaxy apparent to absolute magnitudes. - - Nominally, magnitudes are r-band magnitudes, but this function - is so simple it doesn't matter. - - Just applies a distance correction - no k-correction. - + Convert apparent to absolute magnitudes using the distance modulus only. + + Subtracts ``5 * log10(D_L / 10 pc)`` from the apparent magnitude, where + D_L is the luminosity distance in Mpc. No k-correction is applied. + Args: - App (float or array of floats): apparent galaxy luminosities - zs (float or array of floats): redshifts of galaxies - + App (float or np.ndarray): apparent magnitude(s) m_r. + zs (float or np.ndarray): redshift(s) of the galaxies. + Returns: - AbsoluteMags: NAbs x NZ array of magnitudes, where these - are the dimensions of the inputs + AbsoluteMags: absolute magnitude(s). Scalar if both inputs are + scalar; 1-D array if one input is scalar and one is an array; + 2-D array of shape (NApp, Nz) if both are arrays (computed via + ``np.outer``). """ # calculates luminosity distances (Mpc) @@ -1332,20 +1440,20 @@ def SimpleAbsoluteMags(App,zs): def SimpleApparentMags(Abs,zs): """ - Function to convert galaxy absolue to apparent magnitudes. - - Nominally, magnitudes are r-band magnitudes, but this function - is so simple it doesn't matter. - - Just applies a distance correction - no k-correction. - + Convert absolute to apparent magnitudes using the distance modulus only. + + Adds ``5 * log10(D_L / 10 pc)`` to the absolute magnitude, where + D_L is the luminosity distance in Mpc. No k-correction is applied. + Args: - Abs (float or array of floats): intrinsic galaxy luminosities - zs (float or array of floats): redshifts of galaxies - + Abs (float or np.ndarray): absolute magnitude(s) M_r. + zs (float or np.ndarray): redshift(s) of the galaxies. + Returns: - ApparentMags: NAbs x NZ array of magnitudes, where these - are the dimensions of the inputs + ApparentMags: apparent magnitude(s). Scalar if both inputs are + scalar; 1-D array if one input is scalar and one is an array; + 2-D array of shape (NAbs, Nz) if both are arrays (computed via + ``np.outer``). """ # calculates luminosity distances (Mpc) @@ -1373,15 +1481,23 @@ def SimpleApparentMags(Abs,zs): def p_unseen_Marnoch(zvals,plot=False): """ - Returns probability of a hist being unseen in typical VLT - observations. - - Inputs: - zvals [float, array]: array of redshifts - + Return the probability that an FRB host galaxy is unseen in typical VLT observations. + + Digitises Figure 3 of Marnoch et al. 2023 (MNRAS 525, 994), which shows + p(U | z) — the cumulative probability that a host galaxy at redshift z + falls below the VLT/FORS2 R-band detection limit. A cubic polynomial is + fit to the digitised curve and evaluated at the requested redshifts. + Values are clamped to [0, 1]. + + Args: + zvals (float or np.ndarray): redshift(s) at which to evaluate p(U | z). + plot (bool): if True, save a diagnostic comparison plot of the raw + digitised data, linear interpolation, and polynomial fit to + ``p_unseen.pdf``. Defaults to False. + Returns: - fitv [float, array]: p(Unseen) for redshift zvals - + fitv (np.ndarray): p(U | z) evaluated at each element of ``zvals``, + clamped to the range [0, 1]. """ # approx digitisation of Figure 3 p(U|z) # from Marnoch et al. @@ -1426,7 +1542,20 @@ def p_unseen_Marnoch(zvals,plot=False): def simplify_name(TNSname): """ - Simplifies an FRB name to basics + Reduce a TNS FRB name to a six-character YYMMDD[L] identifier. + + Strips the leading ``FRB`` prefix (if present) and the year's + century digits, retaining only the six-digit date plus any + trailing letter suffix, to allow case-insensitive matching + between survey entries and external FRB catalogues. + + Args: + TNSname (str): FRB name in TNS format, e.g. ``'FRB20180924B'`` + or ``'20180924B'``. + + Returns: + name (str): simplified six-character identifier, e.g. ``'180924B'`` + (six digits plus optional letter). """ # reduces all FRBs to six integers @@ -1447,10 +1576,21 @@ def simplify_name(TNSname): def matchFRB(TNSname,survey): """ - Gets the FRB id from the survey list - Returns None if not in the survey - Used to match properties between a survey - and other FRB libraries + Find the index of an FRB in a survey by TNS name. + + Uses ``simplify_name`` to normalise both the query name and the survey + entries, so minor formatting differences (century digits, trailing + letters) do not prevent a match. + + Args: + TNSname (str): TNS name of the FRB to look up, e.g. + ``'FRB20180924B'``. + survey (Survey): loaded survey object whose ``frbs["TNS"]`` column + contains TNS names of detected FRBs. + + Returns: + int or None: index into ``survey.frbs`` of the first matching FRB, + or ``None`` if the FRB is not found in the survey. """ name = simplify_name(TNSname) @@ -1480,13 +1620,24 @@ def matchFRB(TNSname,survey): def plot_frb(name,ralist,declist,plist,opfile): """ - does an frb - - absolute [bool]: if True, treats rel_error as an absolute value - in arcseconds - - clist: list of astropy coordinates - plist: list of p(O|x) for candidates hosts + Plot an FRB localisation and its PATH host galaxy candidates. + + Produces a scatter plot showing the FRB position and a set of + deviated/sampled positions colour-coded by their PATH posterior + P(O|x), overlaid with circles representing candidate host galaxies + scaled by their angular size. All coordinates are shown in arcseconds + relative to the FRB position. + + Args: + name (str): TNS FRB name (e.g. ``'FRB20180924B'``), used to load + the FRB object and the corresponding PATH candidate table. + ralist (np.ndarray): right ascension values (degrees) of deviated + FRB positions to plot, colour-coded by ``plist``. + declist (np.ndarray): declination values (degrees) of deviated + FRB positions. + plist (np.ndarray): PATH posterior values P(O|x) for the deviated + positions, used to set the colour scale. + opfile (str): output file path for the saved figure. """ from frb.frb import FRB @@ -1543,15 +1694,25 @@ def plot_frb(name,ralist,declist,plist,opfile): def pUgm(mag,mean,width): """ - Function to describe probability of a galaxy being unidentified - in an optical image as a function of its magnitude - + Return the probability that a galaxy is undetected as a function of magnitude. + + Models the survey detection completeness as a logistic function that + transitions from ~0 (bright, always detected) to ~1 (faint, never + detected) with a smooth rolloff centred on ``mean``: + + p(U | m) = 1 / (1 + exp((mean - m) / width)) + Args: - mag (float or array of floats): magnitude(s) at which - to evaluate the function - mean: magnitude at which the probability is 50% - width: characteristic width of transition from 0-50 and - 50-100 % + mag (float or np.ndarray): r-band apparent magnitude(s) at which to + evaluate the detection-failure probability. + mean (float): magnitude at which p(U | m) = 0.5 (the 50% completeness + limit of the survey). + width (float): characteristic width of the completeness rolloff in + magnitudes. Smaller values give a sharper transition. + + Returns: + pU (float or np.ndarray): probability of non-detection at each + magnitude in ``mag``, in the range [0, 1]. """ # converts to a number relative to the mean. Will be weird for mags < 0. diff --git a/zdm/optical_numerics.py b/zdm/optical_numerics.py index 902ce5c5..d355a567 100644 --- a/zdm/optical_numerics.py +++ b/zdm/optical_numerics.py @@ -1,7 +1,30 @@ """ -Contains files related to numerical optimisation -of FRB host galaxy parameters. Similar to iteration.py -for the grid. +Numerical routines for evaluating and optimising FRB host galaxy magnitude models. + +This module is the numerical workhorse for the PATH integration in zdm, +analogous to ``iteration.py`` for the zdm grid. It provides: + +- **``function``** — objective function passed to ``scipy.optimize.minimize`` + that evaluates a goodness-of-fit statistic for a given set of host model + parameters against the CRAFT ICS optical data. + +- **``calc_path_priors``** — inner loop that runs PATH on a list of FRBs + across one or more surveys/grids, collecting priors, posteriors, and + undetected-host probabilities for each FRB. + +- **``run_path``** — runs the PATH algorithm for a single named FRB, + loading its candidate host galaxies from the ``frb`` package data + and applying colour corrections to convert to r-band. + +- **``calculate_likelihood_statistic``** and **``calculate_ks_statistic``** + — goodness-of-fit statistics comparing the model apparent magnitude prior + to the observed PATH posteriors across all FRBs. + +- **``make_cumulative_plots``** — plotting routine for visualising + cumulative magnitude distributions for one or more models simultaneously. + +- **``make_wrappers``**, **``make_cdf``**, **``flatten``**, + **``get_cand_properties``** — supporting utilities. """ import os @@ -19,20 +42,37 @@ def function(x,args): """ - This is a function for input into the scipi.optimize.minimise routine. - - It calculates a set of PATH priors for that model, and then calculates - a test statistic for that set. - - Args: - frblist: list of TNS FRB names - ss: list of surveys in which the FRB may exist - gs: list of grids corresponding to those surveys - model: optical model class which takes arguments x to be minimised. i.e. - the function call model.AbsPrior = x must fully specify the model. - istat [int]: which stat to use? 0 = ks stat. 1 = mak likelihood - - + Objective function for ``scipy.optimize.minimize`` over host model parameters. + + Updates the host magnitude model with parameter vector ``x``, runs PATH + on all FRBs, computes the chosen goodness-of-fit statistic, and returns + a scalar value suitable for minimisation (i.e. smaller is better). + + Parameters + ---------- + x : np.ndarray + Parameter vector passed to ``model.init_args(x)``. Its meaning + depends on the model (e.g. absolute magnitude bin weights for + ``simple_host_model``, or ``fSFR`` for ``loudas_model``). + args : list + Packed argument tuple with the following elements, in order: + + - ``frblist`` (list of str): TNS names of FRBs to evaluate. + - ``ss`` (list of Survey): surveys in which the FRBs may appear. + - ``gs`` (list of Grid): zdm grids corresponding to those surveys. + - ``model``: host magnitude model instance (must implement + ``init_args``). + - ``POxcut`` (float or None): if not None, restrict the statistic + to FRBs whose best host candidate has P(O|x) > POxcut. + - ``istat`` (int): statistic to use — 0 for KS-like statistic, + 1 for maximum-likelihood (returned as negative log-likelihood + so that minimisation maximises the likelihood). + + Returns + ------- + stat : float + Goodness-of-fit statistic (smaller is better). For ``istat=1`` + this is the negative log-likelihood. """ frblist = args[0] @@ -80,11 +120,32 @@ def make_wrappers(model,grids): return wrappers -def make_cdf(xs,ys,ws,norm = True): +def make_cdf(xs,ys,ws,norm=True): """ - makes a cumulative distribution in terms of - the x-values x, observed values y, and weights w - + Build a weighted empirical CDF evaluated on a fixed grid. + + For each grid point ``x`` in ``xs``, accumulates the weights ``ws[i]`` + of all observations ``ys[i]`` that fall below ``x``. The result is a + non-decreasing array that can be compared to a model prior CDF. + + Parameters + ---------- + xs : np.ndarray + Grid of x values at which to evaluate the CDF (e.g. apparent + magnitude bin centres). Must be sorted in ascending order. + ys : array-like + Observed data values (e.g. host galaxy apparent magnitudes). + ws : array-like + Weight for each observation in ``ys`` (e.g. PATH posteriors P_Ox). + Must have the same length as ``ys``. + norm : bool, optional + If True (default), normalise the CDF so that its maximum value + is 1. Set to False to preserve the raw cumulative weight sum. + + Returns + ------- + cdf : np.ndarray, shape (len(xs),) + Weighted empirical CDF evaluated at each point in ``xs``. """ cdf = np.zeros([xs.size]) for i,y in enumerate(ys): @@ -97,24 +158,68 @@ def make_cdf(xs,ys,ws,norm = True): def calc_path_priors(frblist,ss,gs,wrappers,verbose=True,usemodel=True,P_U=0.1): """ - Inner loop. Gets passed model parameters, but assumes everything is - initialsied from there. - - Inputs: - FRBLIST: list of FRBs to retrieve data for - ss: list of surveys modelling those FRBs (searches for FRB in data) - gs: list of zDM grids modelling those surveys - wrappers: list of optical wrapper class objects used to calculate priors on magnitude - verbose (bool): Set to true to generate further output - - Returns: - Number of FRBs fitted - AppMags: list of apparent magnitudes used internally in the model - allMagPriors: summed array of magnitude priors calculated by the model - allObsMags: list of observed magnitudes of candidate hosts - allPOx: list of posterior probabilities calculated by the model - allPU: summed values of unobserved prior - allPUx: summed values of posterior of being unobserved + Run PATH on a list of FRBs and return priors, posteriors, and P_U values. + + For each FRB in ``frblist``, searches all surveys in ``ss`` for a match, + computes the zdm-derived apparent magnitude prior (if ``usemodel=True``), + and runs PATH to produce host association posteriors. Results for all FRBs + are collected into parallel lists (one entry per FRB). + + Also writes a CSV file ``allgalaxies.csv`` (if it does not already exist) + containing the magnitude and VLT/FORS2 R-band columns for all candidate + host galaxies across all FRBs. + + Parameters + ---------- + frblist : list of str + TNS names of FRBs to process (e.g. ``['FRB20180924B', ...]``). + ss : list of Survey + Survey objects to search for each FRB. The first survey containing + a given FRB is used. + gs : list of Grid + zdm grids corresponding to each survey in ``ss``. + wrappers : list of model_wrapper + One ``model_wrapper`` per grid (from ``make_wrappers``), used to + compute DM-dependent apparent magnitude priors. + verbose : bool, optional + If True, print a warning for each FRB not found in any survey. + Defaults to True. + usemodel : bool, optional + If True, use the zdm-derived magnitude prior from ``wrappers`` and + estimate P_U from the model. If False, use PATH's built-in inverse + prior and the supplied fixed ``P_U``. Defaults to True. + P_U : float, optional + Fixed prior probability that the host galaxy is undetected. Only + used when ``usemodel=False``. Defaults to 0.1. + + Returns + ------- + nfitted : int + Number of FRBs successfully matched to a survey and processed. + AppMags : np.ndarray + Internal apparent magnitude grid (from the last processed wrapper). + allMagPriors : list of np.ndarray + One array per FRB giving p(m_r | DM_EG) on the ``AppMags`` grid. + Entries are ``None`` when ``usemodel=False``. + allObsMags : list of np.ndarray + One array per FRB listing the r-band magnitudes of PATH candidate + host galaxies. + allPO : list of np.ndarray + One array per FRB giving the PATH prior P_O for each candidate. + allPOx : list of np.ndarray + One array per FRB giving the PATH posterior P(O|x) for each candidate. + allPU : list of float + Prior P_U (probability of unseen host) for each FRB. + allPUx : list of float + Posterior P(U|x) (probability host is unseen, given data) for each FRB. + sumPU : float + Sum of ``allPU`` across all FRBs. + sumPUx : float + Sum of ``allPUx`` across all FRBs. + frbs : list of str + TNS names of the FRBs that were successfully matched and processed. + dms : list of float + Extragalactic DM (pc cm⁻³) for each FRB in ``frbs``. """ NFRB = len(frblist) @@ -243,26 +348,48 @@ def calc_path_priors(frblist,ss,gs,wrappers,verbose=True,usemodel=True,P_U=0.1): def calculate_likelihood_statistic(NFRB,AppMags,AppMagPriors,ObsMags,ObsPosteriors,PUobs, PUprior,plotfile=None,POxcut=None): """ - Calculates a likelihood for each of the FRBs, and returns the log-likelihood. - - The inputs are in two categories. One is a form of lists of lists, where there is one list for - each FRB, and one entry in that list for each host galaxy candidate. Size is NFRB x NCAND - - The other input is where the length of the list matches the internal array size used to - calculate priors on host magnitudes. Size is either NMAG or NFRBxNMAG - - Inputs: - AppMags [array of floats: NMAG]: array listing apparent magnitudes used to calculate priors - AppMagPrior [array of floats NFRB xNMAG]: array giving prior on AppMags - ObsMags: list of lists of floats giving observed magnitudes m_r of host candidates - ObsPosteriors: list of lists float of posterior values P(O|x) corresponding to ObsMags - PUobs [float]: posterior on unseen probability - PUprior [float]: prior on PU - plotfile: set to name of output file for comparison plot - POxcut: if not None, cut data to fixed POx. Used to simulate current techniques - - Returns: - log likelihood of the observation + Compute the total log-likelihood of the observed PATH posteriors given the model prior. + + For each FRB, evaluates log10(Σ P(O_i|x) / rescale + P_U_prior), where the + rescale factor accounts for PATH's internal renormalisation of posteriors + relative to the model prior. Summing over all FRBs gives the total + log-likelihood returned to the caller. + + Parameters + ---------- + NFRB : int + Number of FRBs to sum over. + AppMags : np.ndarray, shape (NMAG,) + Apparent magnitude grid used to compute the model prior (not used + directly in this function, but kept for API consistency with + ``calculate_ks_statistic``). + AppMagPriors : list of np.ndarray, length NFRB + Model prior p(m_r | DM_EG) on the ``AppMags`` grid, one array per FRB. + ObsMags : list of np.ndarray, length NFRB + Observed r-band magnitudes of PATH candidate host galaxies, one array + per FRB (length NCAND varies by FRB). + ObsPosteriors : list of np.ndarray, length NFRB + PATH posterior P(O_i|x) for each candidate, one array per FRB. + PUobs : list of float, length NFRB + PATH posterior P(U|x) — probability that the true host is undetected — + for each FRB, as returned by PATH after renormalisation. + PUprior : list of float, length NFRB + Model prior P_U for each FRB, as estimated by + ``wrapper.estimate_unseen_prior()``. + plotfile : str or None, optional + If provided, save a diagnostic plot comparing prior and posterior + magnitude distributions to this file path. Defaults to None. + POxcut : float or None, optional + If not None, restrict the statistic to FRBs whose maximum P(O|x) + exceeds this threshold (simulates requiring a confident host ID). + Defaults to None. + + Returns + ------- + stat : float + Total log10-likelihood summed over all NFRB FRBs. Larger values + indicate a better fit. Multiply by -1 for use as a minimisation + objective. """ # calculates log-likelihood of observation stat=0 @@ -296,35 +423,58 @@ def calculate_likelihood_statistic(NFRB,AppMags,AppMagPriors,ObsMags,ObsPosterio def flatten(xss): """ - Turns a list of lists into a single list + Flatten a list of lists into a single flat list. """ return [x for xs in xss for x in xs] def calculate_ks_statistic(NFRB,AppMags,AppMagPriors,ObsMags,ObsPosteriors,PUobs, PUprior,plotfile=None,POxcut=None,plotlabel=None,abc=None,tag=""): """ - Calculates a ks-like statistic to be proxy for goodness-of-fit - We must set each AppMagPriors to 1.-PUprior at the limiting magnitude for each observation, - and sum the ObsPosteriors to be equal to 1.-PUobs at that magnitude. - Then these are what gets summed. - - This can be readily done by combining all ObsMags and ObsPosteriors into a single long list, - since this should already be correctly normalised. Priors require their own weight. - - Inputs: - AppMags: array listing apparent magnitudes - AppMagPriors: list of lists giving priors on AppMags for each FRB - ObsMags: list of observed magnitudes - ObsPosteriors: list of posterior values corresponding to ObsMags - PUobs: posterior on unseen probability - PUprior: prior on PU - plotfile: set to name of output file for comparison plot - POxcut: if not None, cut data to fixed POx. Used to simulate current techniques - abc [None]: add label, e.g. (a), to upper left - tag [string]: string to prefix labels - - Returns: - k-like statistic of biggest obs/prior difference + Compute a KS-like goodness-of-fit statistic between model prior and observed posteriors. + + Builds cumulative magnitude distributions for both the model prior and the + PATH posteriors, normalised by the number of FRBs, and returns the maximum + absolute difference between them — analogous to the KS test statistic. + + Optionally produces a plot comparing the two cumulative distributions. + + Parameters + ---------- + NFRB : int + Number of FRBs used for normalisation. + AppMags : np.ndarray, shape (NMAG,) + Apparent magnitude grid on which priors are defined. + AppMagPriors : list of np.ndarray, length NFRB + Model prior p(m_r | DM_EG) on ``AppMags``, one array per FRB. + ObsMags : list of np.ndarray, length NFRB + Observed r-band magnitudes of PATH candidate galaxies, one per FRB. + ObsPosteriors : list of np.ndarray, length NFRB + PATH posteriors P(O_i|x) for each candidate, one array per FRB. + PUobs : list of float + Posterior P(U|x) for each FRB (not used directly in the statistic, + kept for API consistency). + PUprior : list of float + Prior P_U for each FRB (not used directly, kept for API consistency). + plotfile : str or None, optional + If provided, save a CDF comparison plot to this path. Defaults to None. + POxcut : float or None, optional + If not None, restrict to candidates with P(O|x) > POxcut and + normalise both CDFs to unity (simulates the approach of selecting + only confidently identified hosts). Defaults to None. + plotlabel : str or None, optional + Text label placed in the centre-bottom of the plot. Defaults to None. + abc : str or None, optional + Panel label (e.g. ``'(a)'``) placed in the upper-left corner of the + figure in figure-coordinate space. Defaults to None. + tag : str, optional + String prefix added to the legend labels ``"Observed"`` and + ``"Prior"``. Defaults to ``""``. + + Returns + ------- + stat : float + Maximum absolute difference between the observed and prior cumulative + distributions. Smaller values indicate a better fit. """ # sums the apparent mag priors over all FRBs to create a cumulative distribution fAppMagPriors = np.zeros([len(AppMags)]) @@ -402,28 +552,58 @@ def make_cumulative_plots(NMODELS,NFRB,AppMags,AppMagPriors,ObsMags,ObsPosterior PUprior,plotfile,plotlabel,POxcut=None,abc=None,onlyobs=None, greyscale=[],addpriorlabel=True): """ - Creates cumulative plots of KS-like behaviour for multiple fit outcomes - - Inputs: see "calculate_ks_statistic" except: - - NMODELS (int): number of models to plot - - abc remains unchanged - - NFRB remains unchanged - - plotfile remains unchanged - - onlyobs (int): if not None, only plot observed distribution for this case - - all other parameters have a leading dimension of NMODELS - - Inputs from "calculate_ks_statistic" with extra NMODELS dimension: - AppMags: array listing apparent magnitudes - AppMagPriors: list of lists giving priors on AppMags for each FRB - ObsMags: list of observed magnitudes - ObsPosteriors: list of posterior values corresponding to ObsMags - PUobs: posterior on unseen probability - PUprior: prior on PU - POxcut: if not None, cut data to fixed POx. Used to simulate current techniques - greyscale (list of ints): if present, defines which models to plot as background greyscales - addpriorlabel (bool): if True, add "prior" to label of priors - Returns: - None + Plot cumulative apparent magnitude distributions for multiple host models on one figure. + + Computes the same normalised prior and observed CDFs as + ``calculate_ks_statistic``, but for ``NMODELS`` models simultaneously, + overlaying them on a single figure with distinct line styles. + + All list-valued parameters that appear in ``calculate_ks_statistic`` + gain an additional leading dimension of size ``NMODELS`` here. + + Parameters + ---------- + NMODELS : int + Number of models to plot. + NFRB : list of int, length NMODELS + Number of FRBs for each model, used for normalisation. + AppMags : list of np.ndarray, length NMODELS + Apparent magnitude grid for each model. + AppMagPriors : list of lists of np.ndarray, shape (NMODELS, NFRB, NMAG) + Model prior p(m_r | DM_EG) for each model and FRB. + ObsMags : list of lists of np.ndarray, shape (NMODELS, NFRB, NCAND) + Observed candidate magnitudes for each model and FRB. + ObsPosteriors : list of lists of np.ndarray, shape (NMODELS, NFRB, NCAND) + PATH posteriors P(O_i|x) for each model and FRB. + PUobs : list, length NMODELS + Posterior P(U|x) per model (not used directly in the plot). + PUprior : list, length NMODELS + Prior P_U per model (not used directly in the plot). + plotfile : str + Output file path for the saved figure. + plotlabel : list of str, length NMODELS + Legend label prefix for each model. + POxcut : float or None, optional + If not None, restrict to candidates with P(O|x) > POxcut and + normalise CDFs to unity. Defaults to None. + abc : str or None, optional + Panel label (e.g. ``'(a)'``) placed in the upper-left corner in + figure-coordinate space. Defaults to None. + onlyobs : int or None, optional + If not None, only draw the observed CDF for the model with this + index (useful when all models share the same observations). The + observed line is then labelled ``"Observed"`` without a model prefix. + Defaults to None. + greyscale : list of int, optional + Indices of models whose observed CDF should additionally be drawn + in grey (for background reference). Defaults to ``[]``. + addpriorlabel : bool, optional + If True (default), append ``": Prior"`` to each model's legend entry. + Set to False to use only ``plotlabel[imodel]`` as the label. + + Returns + ------- + None """ # arrays to hold created observed and prior distributions @@ -532,13 +712,20 @@ def make_cumulative_plots(NMODELS,NFRB,AppMags,AppMagPriors,ObsMags,ObsPosterior def get_cand_properties(frblist): """ - Returns properties of galaxy candidates for FRBs - + Load PATH candidate host galaxy properties for a list of FRBs. + + Reads the pre-generated PATH CSV files from the ``frb`` package data + directory (``frb/data/Galaxies/PATH/_PATH.csv``) and extracts + the columns ``['ang_size', 'mag', 'ra', 'dec', 'separation']`` for + each FRB. + Args: - frblist: list of strings giving FRB names - + frblist (list of str): TNS FRB names (e.g. ``['FRB20180924B', ...]``). + Returns: - all_candidates: list of pandas dataframes containing candidate info + all_candidates (list of pd.DataFrame): one DataFrame per FRB, + each with columns ``ang_size``, ``mag``, ``ra``, ``dec``, + and ``separation``. """ all_candidates=[] @@ -556,16 +743,54 @@ def get_cand_properties(frblist): all_candidates.append(candidates) return all_candidates -def run_path(name,P_U=0.1,usemodel = False, sort=False): +def run_path(name,P_U=0.1,usemodel=False,sort=False): """ - evaluates PATH on an FRB - - Args: - name [string]: TNS name of FRB - P_U [float]: unseen prior - usemodel [bool]: if True, use user-defined P_O|x model - sort [bool]: if True, sort candidates by posterior - + Run the PATH algorithm on a single FRB and return host association results. + + Loads the FRB object and its pre-generated PATH candidate table from the + ``frb`` package, applies colour corrections to convert candidate magnitudes + to r-band (using fixed offsets: I → R: +0.65, g → R: −0.65), sets up the + FRB localisation ellipse and offset prior, and evaluates PATH posteriors. + + The magnitude prior used for the candidates is: + + - ``usemodel=False``: PATH's built-in ``'inverse'`` prior (uniform in log + surface density). + - ``usemodel=True``: the ``'user'`` prior, which must be set externally by + pointing ``pathpriors.USR_raw_prior_Oi`` at a ``model_wrapper`` method + before calling this function (typically done by + ``wrapper.init_path_raw_prior_Oi``). + + The offset prior is always the ``'exp'`` model from PATH's ``'adopted'`` + standard priors, with scale 0.5 arcsec. + + Parameters + ---------- + name : str + TNS name of the FRB (e.g. ``'FRB20180924B'``). + P_U : float, optional + Prior probability that the true host galaxy is undetected. Defaults + to 0.1. + usemodel : bool, optional + If True, use the externally set user prior for candidate magnitudes. + Defaults to False. + sort : bool, optional + If True, sort the returned arrays by P(O|x) in ascending order. + Defaults to False. + + Returns + ------- + P_O : np.ndarray + Prior probability P(O_i) for each candidate host galaxy. + P_Ox : np.ndarray + Posterior probability P(O_i|x) for each candidate. + P_Ux : float + Posterior probability P(U|x) that the true host is undetected. + mags : np.ndarray + R-band apparent magnitudes of the candidates (after colour correction). + ptbl : pd.DataFrame + Full PATH candidate table loaded from the CSV file, with an + additional ``'frb'`` column set to ``name``. """ ######### Loads FRB, and modifes properties ######### diff --git a/zdm/scripts/Path/estimate_path_priors.py b/zdm/scripts/Path/estimate_path_priors.py index 36cf0d84..b3b1155e 100644 --- a/zdm/scripts/Path/estimate_path_priors.py +++ b/zdm/scripts/Path/estimate_path_priors.py @@ -1,10 +1,39 @@ """ -Script showing how to use zDM as priors for CRAFT -host galaxy magnitudes. +Estimate zdm-informed PATH priors for CRAFT/ICS FRB host galaxies. -It requirses the FRB and astropath modules to be installed. +This script demonstrates how to incorporate zdm-derived p(z|DM) predictions +as priors for the PATH (Probabilistic Association of Transients to their Hosts) +algorithm applied to CRAFT ICS FRBs. -This does NOT include optimisation of any parameters +For each FRB in the CRAFT ICS sample (`opt.frblist`), the script runs PATH +twice and compares results: + +1. **Baseline run**: PATH with a flat (uninformative) prior on host galaxy + apparent magnitude, and a fixed prior P_U=0.1 on the host being below + the detection threshold. + +2. **zdm-informed run**: PATH using a physically motivated prior on host + apparent magnitude derived from the Marnoch+2023 host galaxy luminosity + model combined with the zdm p(z|DM_EG) probability distribution. The + probability P_U that the true host is undetected is also estimated from + the model rather than set by hand. + +The output is a weighted histogram of posterior host galaxy apparent +magnitudes (P_Ox) across all FRBs, saved to ``posterior_pOx.png``. + +Note: This script does NOT optimise any zdm or host galaxy model parameters. +It uses the CRAFT_ICS_1300 survey grid with default zdm parameter values. + +Requirements +------------ +- ``astropath`` package (PATH implementation) +- ``frb`` package (FRB utilities and optical data) +- PATH-compatible optical data for each FRB in ``opt.frblist`` + +References +---------- +- Marnoch et al. 2023, MNRAS 525, 994 (host galaxy luminosity model) +- Macquart et al. 2020 (Macquart relation / p(DM|z)) """ #standard Python imports @@ -24,7 +53,35 @@ def calc_path_priors(): """ - Loops over all ICS FRBs + Run PATH on all CRAFT ICS FRBs with and without zdm-derived priors. + + Initialises a zdm grid for the CRAFT_ICS_1300 survey and the Marnoch+2023 + host galaxy luminosity model. For each FRB in ``opt.frblist``: + + - Matches the FRB to the CRAFT_ICS_1300 survey to retrieve its + extragalactic dispersion measure (DM_EG). + - Runs PATH with a flat apparent-magnitude prior and fixed P_U=0.1 + (``usemodel=False``), giving baseline posteriors P_Ox1. + - Uses the zdm model to compute a physically motivated prior on apparent + host magnitude, p(m_r | DM_EG), via ``wrapper.init_path_raw_prior_Oi``, + and estimates P_U from the fraction of the magnitude prior that falls + below the survey detection limit via ``wrapper.estimate_unseen_prior``. + - Runs PATH again with the zdm-derived prior (``usemodel=True``) to give + updated posteriors P_Ox2. + + After processing all FRBs, produces a weighted histogram of the posterior + host apparent magnitudes (P_Ox2) across the whole sample and saves it to + ``posterior_pOx.png``. + + Notes + ----- + FRBs not found in the CRAFT_ICS_1300 survey (e.g. because they were + detected by a different instrument configuration) are skipped with a + warning. + + The zdm model parameters are held fixed at their default values; no + parameter optimisation is performed here. See + ``optimise_host_priors.py`` for the equivalent script with optimisation. """ frblist = opt.frblist diff --git a/zdm/scripts/Path/optimise_host_priors.py b/zdm/scripts/Path/optimise_host_priors.py index a1da24fb..99d9569b 100644 --- a/zdm/scripts/Path/optimise_host_priors.py +++ b/zdm/scripts/Path/optimise_host_priors.py @@ -1,19 +1,50 @@ """ -This file illustrates how to optimise the host prior -distribution by fitting to CRAFT ICS optical observations. -It fits a model of absolute galaxy magnitude distributions, -uses zDM to predict redshifts and hence apparent magntidues, -runs PATH using that prior, and tries to get priors to match posteriors. +Optimise FRB host galaxy magnitude priors using zdm predictions and PATH. -WARNING: this is NOT the optimal method! That would require using -a catalogue of galaxies to sample from to generate fake optical fields. -But nonetheless, this tests the power of estimating FRB host galaxy -contributions using zDM to set priors for apparent magnitudes. +This script fits a parametric model of FRB host galaxy absolute magnitude +distributions to the CRAFT ICS optical observations. It works by: -WARNING2: To do this properly also requires inputting the posterior POx -for host galaxies into zDM! This simulation does not do that either. +1. Initialising zdm grids for the three CRAFT ICS survey bands (892, 1300, + and 1632 MHz) using the HoffmannHalo25 parameter state. +2. Constructing a host galaxy model (``simple`` or ``loudas``) that predicts + apparent r-band magnitudes by convolving the absolute magnitude distribution + with the zdm p(z|DM_EG) redshift prior, optionally including a k-correction. +3. Running PATH with those zdm-derived apparent magnitude priors to obtain + posterior host association probabilities P_Ox for each CRAFT ICS FRB. +4. Optimising the model parameters with ``scipy.optimize.minimize`` by + minimising either a maximum-likelihood statistic or a KS-like goodness-of-fit + statistic against the observed PATH posteriors. -WARNING3: this program can take a while to run, if optimising the simple model. +After optimisation the script: + +- Saves the best-fit parameters to ``_output/best_fit_params.npy``. +- Plots the predicted vs observed apparent magnitude distributions for the + best-fit model (``best_fit_apparent_magnitudes.png``). +- Re-runs PATH with the original (flat) priors for comparison and produces a + scatter plot of best-fit vs original posteriors + (``Scatter_plot_comparison.png``). + +Limitations +----------- +- The optimal approach would sample galaxy candidates from a real photometric + catalogue to construct proper optical fields; this script uses a parametric + model instead. +- Host identification posteriors (P_Ox) are not fed back into the zdm + likelihood; a self-consistent joint fit is not performed. +- Runtime can be significant when optimising the ``simple`` model (10 free + parameters by default). + +Usage +----- +Set ``minimise = True`` (default) to run the optimiser, or ``False`` to load +previously saved parameters from ``_output/best_fit_params.npy``. +Switch between host models by changing ``modelname`` to ``"simple"`` or +``"loudas"``. + +Requirements +------------ +- ``astropath`` package (PATH implementation) +- ``frb`` package (FRB utilities and optical data) """ @@ -48,13 +79,40 @@ def main(): """ - Main function - Contains outer loop to iterate over parameters - + Optimise host galaxy magnitude model parameters and compare with baseline PATH. + + Workflow: + + 1. Load the CRAFT ICS FRB list and initialise zdm grids for the 892, 1300, + and 1632 MHz survey bands using the HoffmannHalo25 cosmological/FRB state. + 2. Select a host magnitude model (``"simple"`` or ``"loudas"``) and configure + its parameter bounds and initial values. + 3. If ``minimise=True``, call ``scipy.optimize.minimize`` with + ``on.function`` as the objective, minimising either the maximum-likelihood + statistic (``istat=1``) or the KS-like statistic (``istat=0``) over all + CRAFT ICS FRBs. Best-fit parameters are saved to + ``_output/best_fit_params.npy``. + 4. Re-evaluate PATH at the best-fit parameters and compute both the + likelihood and KS statistics; save the apparent magnitude comparison + plot to ``_output/best_fit_apparent_magnitudes.png``. + 5. Re-run PATH with the original flat priors (``usemodel=False``) and save + a scatter plot comparing original vs best-fit P_Ox posteriors to + ``_output/Scatter_plot_comparison.png``. + + Configuration knobs (edit at the top of the function body): + + - ``istat``: 0 = KS statistic, 1 = maximum-likelihood statistic. + - ``dok``: whether to include a k-correction in the apparent magnitude model. + - ``modelname``: ``"simple"`` for the parametric histogram model or + ``"loudas"`` for the Loudas single-parameter model. + - ``POxcut``: optional float (e.g. 0.9) to exclude low-confidence FRBs + from the model comparison. + - ``minimise``: set to ``False`` to skip optimisation and load saved + parameters instead. """ ######### List of all ICS FRBs for which we can run PATH ####### - # hard-coded list of FRBs with PATH data in ice paper + # hard-coded list of FRBs with PATH data in ICE paper frblist=opt.frblist # Initlisation of zDM grid @@ -176,11 +234,31 @@ def main(): plt.close() -def make_cdf_for_plotting(xvals,weights=None): +def make_cdf_for_plotting(xvals, weights=None): """ - Creates a cumulative distribution function - - xvals,yvals: values of data points + Build a step-function CDF suitable for plotting. + + Converts an array of data values (and optional weights) into paired + (x, y) arrays that trace the cumulative distribution as a staircase, + with two points per input value so that horizontal steps are rendered + correctly by matplotlib. + + Parameters + ---------- + xvals : np.ndarray + 1-D array of data values. Will be sorted in ascending order. + weights : np.ndarray, optional + 1-D array of weights with the same length as ``xvals``. If provided, + the CDF is computed as the normalised cumulative sum of the sorted + weights. If ``None``, a uniform CDF over ``N`` points is used, + with steps at ``0, 1/N, 2/N, ..., 1``. + + Returns + ------- + cx : np.ndarray + x-coordinates of the staircase CDF (length ``2 * N``). + cy : np.ndarray + y-coordinates of the staircase CDF (length ``2 * N``). """ N = xvals.size cx = np.zeros([2*N]) diff --git a/zdm/scripts/Path/plot_host_models.py b/zdm/scripts/Path/plot_host_models.py index b40d7a89..7ae403c3 100644 --- a/zdm/scripts/Path/plot_host_models.py +++ b/zdm/scripts/Path/plot_host_models.py @@ -1,16 +1,52 @@ """ -Script showing how to load different FRB host models. +Plot and compare FRB host galaxy magnitude models against CRAFT ICS PATH posteriors. -It does this for both my simple naive model, and Nick Loudas's -model based on mass- or sfr-weightings, and Lachlan's -evolution of galaxy spectra. +This script demonstrates how to load, configure, and visualise the three +available FRB host galaxy magnitude models, then evaluate PATH host association +posteriors for all CRAFT ICS FRBs using each model in turn. -We then evaluate P(O|x) for CRAFT FRBs in the CRAFT 1300 MHz survey +Model comparison +---------------- +Three host magnitude models are loaded and plotted: -It's not the simplest script, but it should show how to do a whole bunch of stuff +- **Simple model** (``opt.simple_host_model``): a parametric histogram of + absolute magnitudes M_r, linearly interpolated, with an optional + k-correction. Parameters are set by hand here (not from a fitted result). +- **Loudas model** (``opt.loudas_model``): predicts apparent magnitudes from + a galaxy luminosity function weighted by stellar mass (``fSFR=0``) or + star-formation rate (``fSFR=1``), based on Nick Loudas's galaxy model. + The mixing parameter ``fSFR`` is varied to show sensitivity. +- **Marnoch model** (``opt.marnoch_model``): predicts apparent magnitudes + from the galaxy spectral evolution model of Marnoch et al. 2023 + (MNRAS 525, 994). -NOTE: this does NOT use the best-fit distributions form the recent paper. +Diagnostic plots produced in ``Plots/`` +----------------------------------------- +- ``simple_model_mags.png``: absolute magnitude prior p(M_r) for the simple + model, showing the interpolated curve and the raw histogram bin values. +- ``loudas_model_mags.png``: apparent magnitude distributions p(m_r) for the + Loudas model at several redshifts, comparing mass- vs SFR-weighted variants. +- ``loudas_fsfr_interpolation.png``: sensitivity of the Loudas model to the + ``fSFR`` mixing parameter at z=0.5, illustrating the full allowed range. +- ``all_model_apparent_mags.png``: side-by-side comparison of p(m_r | z) for + all three models at z = 0.1, 0.5, and 2.0. +- ``all_model_mag_priors_dm.png``: PATH apparent magnitude priors p(m_r | DM) + for all three models at DM = 200, 600, and 1000 pc/cm³, using the + CRAFT_ICS_1300 zdm grid to convert DM to a redshift prior. +- ``posterior_comparison.png``: scatter plot of PATH host posteriors P(O|x) + from the original flat-prior run vs each of the four zdm-informed model + runs, across all CRAFT ICS FRBs. +Note +---- +Parameter values used here are illustrative initial estimates, not best-fit +results from the published analysis. See ``optimise_host_priors.py`` for the +fitting procedure. + +Requirements +------------ +- ``astropath`` package (PATH implementation) +- ``frb`` package (FRB utilities and optical data) """ #standard Python imports @@ -42,7 +78,46 @@ def calc_path_priors(): """ - Loops over all ICS FRBs + Generate diagnostic plots for all host models and compare PATH posteriors. + + Workflow: + + 1. **Model initialisation**: Loads the simple, Loudas (mass-weighted, + ``fSFR=0``), and Marnoch host magnitude models with illustrative + parameter values. + + 2. **Intrinsic magnitude plots**: Plots the absolute magnitude prior + p(M_r) for the simple model and apparent magnitude distributions + p(m_r) for the Loudas model at several redshifts. + + 3. **fSFR sensitivity**: Plots p(m_r | z=0.5) for the Loudas model + across a wide range of ``fSFR`` values to illustrate model behaviour + beyond the physically motivated [0, 1] range. + + 4. **Model comparison at fixed z**: Compares p(m_r | z) across all three + models at z = 0.1, 0.5, and 2.0. + + 5. **DM-dependent priors**: Loads the CRAFT_ICS_1300 zdm grid and wraps + each model in a ``model_wrapper`` to compute the PATH apparent magnitude + prior p(m_r | DM) at DM = 200, 600, and 1000 pc/cm³. + + 6. **PATH evaluation over CRAFT ICS FRBs**: For each FRB in + ``opt.frblist`` that is found in the CRAFT_ICS_1300 survey: + + - Runs PATH with a flat prior (``usemodel=False``, P_U=0.1) as the + baseline. + - Runs PATH four more times, one per zdm-informed model variant + (simple; Loudas mass-weighted; Loudas SFR-weighted; Marnoch), each + with P_U estimated from ``wrapper.estimate_unseen_prior()``. + - Prints diagnostic output for candidate host galaxies that flip above + P_Ox=0.5 relative to the baseline (simple model only). + + 7. **Posterior scatter plot**: Produces ``Plots/posterior_comparison.png`` + showing P(O|x) from each zdm-informed model against the flat-prior + baseline across all FRBs. + + Output files are written to the ``Plots/`` subdirectory, which is created + if it does not already exist. """ opdir = "Plots/"