Coverage for bilby/gw/prior.py: 74%
683 statements
« prev ^ index » next coverage.py v7.6.1, created at 2025-05-06 04:57 +0000
1import os
2import copy
4import numpy as np
5from scipy.interpolate import InterpolatedUnivariateSpline, interp1d
6from scipy.special import hyp2f1
7from scipy.stats import norm
9from ..core.prior import (
10 PriorDict, Uniform, Prior, DeltaFunction, Gaussian, Interped, Constraint,
11 conditional_prior_factory, PowerLaw, ConditionalLogUniform,
12 ConditionalPriorDict, ConditionalBasePrior, BaseJointPriorDist, JointPrior,
13 JointPriorDistError,
14)
15from ..core.utils import infer_args_from_method, logger, random
16from .conversion import (
17 convert_to_lal_binary_black_hole_parameters,
18 convert_to_lal_binary_neutron_star_parameters, generate_mass_parameters,
19 generate_tidal_parameters, fill_from_fixed_priors,
20 generate_all_bbh_parameters,
21 chirp_mass_and_mass_ratio_to_total_mass,
22 total_mass_and_mass_ratio_to_component_masses)
23from .cosmology import get_cosmology, z_at_value
24from .source import PARAMETER_SETS
25from .utils import calculate_time_to_merger
28DEFAULT_PRIOR_DIR = os.path.join(os.path.dirname(__file__), 'prior_files')
31class BilbyPriorConversionError(Exception):
32 pass
35def convert_to_flat_in_component_mass_prior(result, fraction=0.25):
36 """ Converts samples with a defined prior in chirp-mass and mass-ratio to flat in component mass by resampling with
37 the posterior with weights defined as ratio in new:old prior values times the jacobian which for
38 F(mc, q) -> G(m1, m2) is defined as J := m1^2 / mc
40 Parameters
41 ==========
42 result: bilby.core.result.Result
43 The output result complete with priors and posteriors
44 fraction: float [0, 1]
45 The fraction of samples to draw (default=0.25). Note, if too high a
46 fraction of samples are draw, the reweighting will not be applied in
47 effect.
49 """
50 if getattr(result, "priors") is not None:
51 for key in ['chirp_mass', 'mass_ratio']:
52 if key not in result.priors.keys():
53 BilbyPriorConversionError("{} Prior not found in result object".format(key))
54 if isinstance(result.priors[key], Constraint):
55 BilbyPriorConversionError("{} Prior should not be a Constraint".format(key))
56 for key in ['mass_1', 'mass_2']:
57 if not isinstance(result.priors[key], Constraint):
58 BilbyPriorConversionError("{} Prior should be a Constraint Prior".format(key))
59 else:
60 BilbyPriorConversionError("No prior in the result: unable to convert")
62 result = copy.copy(result)
63 priors = result.priors
64 old_priors = copy.copy(result.priors)
65 posterior = result.posterior
67 for key in ['chirp_mass', 'mass_ratio']:
68 priors[key] = Constraint(priors[key].minimum, priors[key].maximum, key, latex_label=priors[key].latex_label)
69 for key in ['mass_1', 'mass_2']:
70 priors[key] = Uniform(priors[key].minimum, priors[key].maximum, key, latex_label=priors[key].latex_label,
71 unit=r"$M_{\odot}$")
73 weights = np.array(result.get_weights_by_new_prior(old_priors, priors,
74 prior_names=['chirp_mass', 'mass_ratio', 'mass_1', 'mass_2']))
75 jacobian = posterior["mass_1"] ** 2 / posterior["chirp_mass"]
76 weights = jacobian * weights
77 result.posterior = posterior.sample(frac=fraction, weights=weights)
79 logger.info("Resampling posterior to flat-in-component mass")
80 effective_sample_size = sum(weights)**2 / sum(weights**2)
81 n_posterior = len(posterior)
82 if fraction > effective_sample_size / n_posterior:
83 logger.warning(
84 "Sampling posterior of length {} with fraction {}, but "
85 "effective_sample_size / len(posterior) = {}. This may produce "
86 "biased results"
87 .format(n_posterior, fraction, effective_sample_size / n_posterior)
88 )
89 result.posterior = posterior.sample(frac=fraction, weights=weights, replace=True)
90 result.meta_data["reweighted_to_flat_in_component_mass"] = True
91 return result
94class Cosmological(Interped):
95 """
96 Base prior class for cosmologically aware distance parameters.
97 """
99 @property
100 def _default_args_dict(self):
101 from astropy import units
102 return dict(
103 redshift=dict(name='redshift', latex_label='$z$', unit=None),
104 luminosity_distance=dict(
105 name='luminosity_distance', latex_label='$d_L$', unit=units.Mpc),
106 comoving_distance=dict(
107 name='comoving_distance', latex_label='$d_C$', unit=units.Mpc))
109 def __init__(self, minimum, maximum, cosmology=None, name=None,
110 latex_label=None, unit=None, boundary=None):
111 """
113 Parameters
114 ----------
115 minimum: float
116 The minimum value for the distribution.
117 maximum: float
118 The maximum value for the distribution.
119 cosmology: Union[None, str, dict, astropy.cosmology.FLRW]
120 The cosmology to use, see `bilby.gw.cosmology.set_cosmology`
121 for more details.
122 If :code:`None`, the default cosmology will be used.
123 name: str
124 The name of the parameter, should be one of
125 :code:`[luminosity_distance, comoving_distance, redshift]`.
126 latex_label: str
127 Latex string to use in plotting, if :code:`None`, this will
128 be set automatically using the name.
129 unit: Union[str, astropy.units.Unit]
130 The unit of the maximum and minimum values, :code:`default=Mpc`.
131 boundary: str
132 The boundary condition to apply to the prior when sampling.
133 """
134 from astropy import units
135 self.cosmology = get_cosmology(cosmology)
136 if name not in self._default_args_dict:
137 raise ValueError(
138 "Name {} not recognised. Must be one of luminosity_distance, "
139 "comoving_distance, redshift".format(name))
140 self.name = name
141 label_args = self._default_args_dict[self.name]
142 if latex_label is not None:
143 label_args['latex_label'] = latex_label
144 if unit is not None:
145 if not isinstance(unit, units.Unit):
146 unit = units.Unit(unit)
147 label_args['unit'] = unit
148 self.unit = label_args['unit']
149 self._minimum = dict()
150 self._maximum = dict()
151 self.minimum = minimum
152 self.maximum = maximum
153 if name == 'redshift':
154 xx, yy = self._get_redshift_arrays()
155 elif name == 'comoving_distance':
156 xx, yy = self._get_comoving_distance_arrays()
157 elif name == 'luminosity_distance':
158 xx, yy = self._get_luminosity_distance_arrays()
159 else:
160 raise ValueError('Name {} not recognized.'.format(name))
161 super(Cosmological, self).__init__(xx=xx, yy=yy, minimum=minimum, maximum=maximum,
162 boundary=boundary, **label_args)
164 @property
165 def minimum(self):
166 return self._minimum[self.name]
168 @minimum.setter
169 def minimum(self, minimum):
170 if (self.name in self._minimum) and (minimum < self.minimum):
171 self._set_limit(value=minimum, limit_dict=self._minimum, recalculate_array=True)
172 else:
173 self._set_limit(value=minimum, limit_dict=self._minimum)
175 @property
176 def maximum(self):
177 return self._maximum[self.name]
179 @maximum.setter
180 def maximum(self, maximum):
181 if (self.name in self._maximum) and (maximum > self.maximum):
182 self._set_limit(value=maximum, limit_dict=self._maximum, recalculate_array=True)
183 else:
184 self._set_limit(value=maximum, limit_dict=self._maximum)
186 def _set_limit(self, value, limit_dict, recalculate_array=False):
187 """
188 Set either of the limits for redshift, luminosity, and comoving distances
190 Parameters
191 ==========
192 value: float
193 Limit value in current class' parameter
194 limit_dict: dict
195 The limit dictionary to modify in place
196 recalculate_array: boolean
197 Determines if the distance arrays are recalculated
198 """
199 cosmology = get_cosmology(self.cosmology)
200 limit_dict[self.name] = value
201 if self.name == 'redshift':
202 limit_dict['luminosity_distance'] = \
203 cosmology.luminosity_distance(value).value
204 limit_dict['comoving_distance'] = \
205 cosmology.comoving_distance(value).value
206 elif self.name == 'luminosity_distance':
207 if value == 0:
208 limit_dict['redshift'] = 0
209 else:
210 limit_dict['redshift'] = z_at_value(
211 cosmology.luminosity_distance, value * self.unit
212 )
213 limit_dict['comoving_distance'] = (
214 cosmology.comoving_distance(limit_dict['redshift']).value
215 )
216 elif self.name == 'comoving_distance':
217 if value == 0:
218 limit_dict['redshift'] = 0
219 else:
220 limit_dict['redshift'] = z_at_value(
221 cosmology.comoving_distance, value * self.unit
222 )
223 limit_dict['luminosity_distance'] = (
224 cosmology.luminosity_distance(limit_dict['redshift']).value
225 )
226 if recalculate_array:
227 if self.name == 'redshift':
228 self.xx, self.yy = self._get_redshift_arrays()
229 elif self.name == 'comoving_distance':
230 self.xx, self.yy = self._get_comoving_distance_arrays()
231 elif self.name == 'luminosity_distance':
232 self.xx, self.yy = self._get_luminosity_distance_arrays()
233 try:
234 self._update_instance()
235 except (AttributeError, KeyError):
236 pass
238 def get_corresponding_prior(self, name=None, unit=None):
239 """
240 Utility method to convert between equivalent redshift, comoving
241 distance, and luminosity distance priors.
243 Parameters
244 ----------
245 name: str
246 The name of the new prior, this should be one of
247 :code:`[luminosity_distance, comoving_distance, redshift]`.
248 unit: Union[str, astropy.units.Unit]
249 The unit of the output prior.
251 Returns
252 -------
253 The corresponding prior.
254 """
255 subclass_args = infer_args_from_method(self.__init__)
256 args_dict = {key: getattr(self, key) for key in subclass_args}
257 self._convert_to(new=name, args_dict=args_dict)
258 if unit is not None:
259 args_dict['unit'] = unit
260 return self.__class__(**args_dict)
262 def _convert_to(self, new, args_dict):
263 args_dict.update(self._default_args_dict[new])
264 args_dict['minimum'] = self._minimum[args_dict['name']]
265 args_dict['maximum'] = self._maximum[args_dict['name']]
267 def _get_comoving_distance_arrays(self):
268 zs, p_dz = self._get_redshift_arrays()
269 dc_of_z = self.cosmology.comoving_distance(zs).value
270 ddc_dz = np.gradient(dc_of_z, zs)
271 p_dc = p_dz / ddc_dz
272 return dc_of_z, p_dc
274 def _get_luminosity_distance_arrays(self):
275 zs, p_dz = self._get_redshift_arrays()
276 dl_of_z = self.cosmology.luminosity_distance(zs).value
277 ddl_dz = np.gradient(dl_of_z, zs)
278 p_dl = p_dz / ddl_dz
279 return dl_of_z, p_dl
281 def _get_redshift_arrays(self):
282 raise NotImplementedError
284 @classmethod
285 def from_repr(cls, string):
286 if "FlatLambdaCDM" in string:
287 logger.warning(
288 "Cosmological priors cannot be loaded from a string. "
289 "If the prior has a name, use that instead."
290 )
291 return string
292 else:
293 return cls._from_repr(string)
295 def get_instantiation_dict(self):
296 from astropy import units
297 from astropy.cosmology.realizations import available
298 instantiation_dict = super().get_instantiation_dict()
299 if self.cosmology.name in available:
300 instantiation_dict['cosmology'] = self.cosmology.name
301 if isinstance(self.unit, units.Unit):
302 instantiation_dict['unit'] = self.unit.to_string()
303 return instantiation_dict
306class UniformComovingVolume(Cosmological):
307 r"""
308 Prior that is uniform in the comoving volume.
310 Note that this does not account for time dilation in the source frame
311 use `bilby.gw.prior.UniformSourceFrame` to include that effect.
313 .. math::
315 p(z) \propto \frac{d V_{c}}{dz}
316 """
318 def _get_redshift_arrays(self):
319 zs = np.linspace(self._minimum['redshift'] * 0.99,
320 self._maximum['redshift'] * 1.01, 1000)
321 p_dz = self.cosmology.differential_comoving_volume(zs).value
322 return zs, p_dz
325class UniformSourceFrame(Cosmological):
326 r"""
327 Prior for redshift which is uniform in comoving volume and source frame
328 time.
330 .. math::
332 p(z) \propto \frac{1}{1 + z}\frac{dV_{c}}{dz}
334 The extra :math:`1+z` is due to doppler shifting of the source frame time.
335 """
337 def _get_redshift_arrays(self):
338 zs = np.linspace(self._minimum['redshift'] * 0.99,
339 self._maximum['redshift'] * 1.01, 1000)
340 p_dz = self.cosmology.differential_comoving_volume(zs).value / (1 + zs)
341 return zs, p_dz
344class UniformInComponentsChirpMass(PowerLaw):
345 r"""
346 Prior distribution for chirp mass which is uniform in component masses.
348 This is useful when chirp mass and mass ratio are sampled while the
349 prior is uniform in component masses.
351 .. math::
353 p({\cal M}) \propto {\cal M}
355 Notes
356 -----
357 This prior is intended to be used in conjunction with the corresponding
358 :code:`bilby.gw.prior.UniformInComponentsMassRatio`.
359 """
361 def __init__(self, minimum, maximum, name='chirp_mass',
362 latex_label=r'$\mathcal{M}$', unit=None, boundary=None):
363 """
364 Parameters
365 ==========
366 minimum : float
367 The minimum of chirp mass
368 maximum : float
369 The maximum of chirp mass
370 name: see superclass
371 latex_label: see superclass
372 unit: see superclass
373 boundary: see superclass
374 """
375 super(UniformInComponentsChirpMass, self).__init__(
376 alpha=1., minimum=minimum, maximum=maximum,
377 name=name, latex_label=latex_label, unit=unit, boundary=boundary)
380class WrappedInterp1d(interp1d):
381 """ A wrapper around scipy interp1d which sets equality-by-instantiation """
382 def __eq__(self, other):
384 for key in self.__dict__:
385 if type(self.__dict__[key]) is np.ndarray:
386 if not np.array_equal(self.__dict__[key], other.__dict__[key]):
387 return False
388 elif key == "_spline":
389 pass
390 elif getattr(self, key) != getattr(other, key):
391 return False
392 return True
395class UniformInComponentsMassRatio(Prior):
396 r"""
397 Prior distribution for chirp mass which is uniform in component masses.
399 This is useful when chirp mass and mass ratio are sampled while the
400 prior is uniform in component masses.
402 .. math::
404 p(q) \propto \frac{(1 + q)^{2/5}}{q^{6/5}}
406 Notes
407 -----
408 This prior is intended to be used in conjunction with the corresponding
409 :code:`bilby.gw.prior.UniformInComponentsChirpMass`.
410 """
412 def __init__(self, minimum, maximum, name='mass_ratio', latex_label='$q$',
413 unit=None, boundary=None, equal_mass=False):
414 """
415 Parameters
416 ==========
417 minimum : float
418 The minimum of mass ratio
419 maximum : float
420 The maximum of mass ratio
421 name: see superclass
422 latex_label: see superclass
423 unit: see superclass
424 boundary: see superclass
425 equal_mass: bool
426 Whether the likelihood being considered is expected to peak at
427 equal masses. If True, the mapping described in Appendix A of
428 arXiv:2111.13619 is used in the :code:`rescale` method.
429 default=False
431 """
432 self.equal_mass = equal_mass
433 super(UniformInComponentsMassRatio, self).__init__(
434 minimum=minimum, maximum=maximum, name=name,
435 latex_label=latex_label, unit=unit, boundary=boundary)
436 self.norm = self._integral(maximum) - self._integral(minimum)
437 qs = np.linspace(minimum, maximum, 1000)
438 self.icdf = WrappedInterp1d(
439 self.cdf(qs), qs, kind='cubic',
440 bounds_error=False, fill_value=(minimum, maximum))
442 @staticmethod
443 def _integral(q):
444 return -5. * q**(-1. / 5.) * hyp2f1(-2. / 5., -1. / 5., 4. / 5., -q)
446 def cdf(self, val):
447 return (self._integral(val) - self._integral(self.minimum)) / self.norm
449 def rescale(self, val):
450 if self.equal_mass:
451 val = 2 * np.minimum(val, 1 - val)
452 resc = self.icdf(val)
453 if resc.ndim == 0:
454 return resc.item()
455 else:
456 return resc
458 def prob(self, val):
459 in_prior = (val >= self.minimum) & (val <= self.maximum)
460 with np.errstate(invalid="ignore"):
461 prob = (1. + val)**(2. / 5.) / (val**(6. / 5.)) / self.norm * in_prior
462 return prob
464 def ln_prob(self, val):
465 with np.errstate(divide="ignore"):
466 return np.log(self.prob(val))
469class AlignedSpin(Interped):
470 r"""
471 Prior distribution for the aligned (z) component of the spin.
473 This takes prior distributions for the magnitude and cosine of the tilt
474 and forms a compound prior using a numerical convolution integral.
476 .. math::
478 \pi(\chi) = \int_{0}^{1} da \int_{-1}^{1} d\cos\theta
479 \pi(a) \pi(\cos\theta) \delta(\chi - a \cos\theta)
481 This is useful when using aligned-spin only waveform approximants.
483 This is an extension of e.g., (A7) of https://arxiv.org/abs/1805.10457.
484 """
486 def __init__(self, a_prior=Uniform(0, 1), z_prior=Uniform(-1, 1),
487 name=None, latex_label=None, unit=None, boundary=None,
488 minimum=np.nan, maximum=np.nan):
489 """
490 Parameters
491 ==========
492 a_prior: Prior
493 Prior distribution for spin magnitude
494 z_prior: Prior
495 Prior distribution for cosine spin tilt
496 name: see superclass
497 latex_label: see superclass
498 unit: see superclass
499 """
500 self.a_prior = a_prior
501 self.z_prior = z_prior
502 chi_min = min(a_prior.maximum * z_prior.minimum,
503 a_prior.minimum * z_prior.maximum)
504 chi_max = a_prior.maximum * z_prior.maximum
505 xx = np.linspace(chi_min, chi_max, 800)
506 a_prior_minimum = a_prior.minimum
507 if a_prior_minimum == 0:
508 a_prior_minimum += 1e-32
509 aas = np.linspace(a_prior_minimum, a_prior.maximum, 1000)
510 yy = [np.trapz(np.nan_to_num(a_prior.prob(aas) / aas *
511 z_prior.prob(x / aas)), aas) for x in xx]
512 super(AlignedSpin, self).__init__(xx=xx, yy=yy, name=name,
513 latex_label=latex_label, unit=unit,
514 boundary=boundary, minimum=minimum,
515 maximum=maximum)
518class ConditionalChiUniformSpinMagnitude(ConditionalLogUniform):
519 r"""
520 This prior characterizes the conditional prior on the spin magnitude given
521 the aligned component of the spin such that the marginal prior is uniform
522 if the distribution of spin orientations is isotropic.
524 .. math::
526 p(a) &= \frac{1}{a_{\max}} \\
527 p(\chi) &= - \frac{\ln(|\chi|)}{2 a_{\max}} \\
528 p(a | \chi) &\propto \frac{1}{a}
529 """
531 def __init__(self, minimum, maximum, name, latex_label=None, unit=None, boundary=None):
532 super(ConditionalChiUniformSpinMagnitude, self).__init__(
533 minimum=minimum, maximum=maximum, name=name, latex_label=latex_label, unit=unit, boundary=boundary,
534 condition_func=self._condition_function)
535 self._required_variables = [name.replace("a", "chi")]
536 self.__class__.__name__ = "ConditionalChiUniformSpinMagnitude"
537 self.__class__.__qualname__ = "ConditionalChiUniformSpinMagnitude"
539 def _condition_function(self, reference_params, **kwargs):
540 return dict(minimum=np.abs(kwargs[self._required_variables[0]]), maximum=reference_params["maximum"])
542 def __repr__(self):
543 return Prior.__repr__(self)
545 def get_instantiation_dict(self):
546 instantiation_dict = Prior.get_instantiation_dict(self)
547 for key, value in self.reference_params.items():
548 if key in instantiation_dict:
549 instantiation_dict[key] = value
550 return instantiation_dict
553class ConditionalChiInPlane(ConditionalBasePrior):
554 r"""
555 This prior characterizes the conditional prior on the in-plane spin magnitude
556 given the aligned component of the spin such that the marginal prior is uniform
557 if the distribution of spin orientations is isotropic.
559 .. math::
561 p(a) &= \frac{1}{a_{\max}} \\
562 p(\chi_\perp) &= \frac{2 N \chi_\perp}{\chi^2 + \chi_\perp^2} \\
563 N^{-1} &= 2 \ln\left(\frac{a_\max}{|\chi|}\right)
564 """
566 def __init__(self, minimum, maximum, name, latex_label=None, unit=None, boundary=None):
567 super(ConditionalChiInPlane, self).__init__(
568 minimum=minimum, maximum=maximum,
569 name=name, latex_label=latex_label,
570 unit=unit, boundary=boundary,
571 condition_func=self._condition_function
572 )
573 self._required_variables = [name[:5]]
574 self._reference_maximum = maximum
575 self.__class__.__name__ = "ConditionalChiInPlane"
576 self.__class__.__qualname__ = "ConditionalChiInPlane"
578 def prob(self, val, **required_variables):
579 self.update_conditions(**required_variables)
580 chi_aligned = abs(required_variables[self._required_variables[0]])
581 return (
582 (val >= self.minimum) * (val <= self.maximum)
583 * val
584 / (chi_aligned ** 2 + val ** 2)
585 / np.log(self._reference_maximum / chi_aligned)
586 )
588 def ln_prob(self, val, **required_variables):
589 with np.errstate(divide="ignore"):
590 return np.log(self.prob(val, **required_variables))
592 def cdf(self, val, **required_variables):
593 r"""
594 .. math::
595 \text{CDF}(\chi_\per) = N ln(1 + (\chi_\perp / \chi) ** 2)
597 Parameters
598 ----------
599 val: (float, array-like)
600 The value at which to evaluate the CDF
601 required_variables: dict
602 A dictionary containing the aligned component of the spin
604 Returns
605 -------
606 (float, array-like)
607 The value of the CDF
609 """
610 self.update_conditions(**required_variables)
611 chi_aligned = abs(required_variables[self._required_variables[0]])
612 return np.maximum(np.minimum(
613 (val >= self.minimum) * (val <= self.maximum)
614 * np.log(1 + (val / chi_aligned) ** 2)
615 / 2 / np.log(self._reference_maximum / chi_aligned)
616 , 1
617 ), 0)
619 def rescale(self, val, **required_variables):
620 r"""
621 .. math::
622 \text{PPF}(\chi_\perp) = ((a_\max / \chi) ** (2x) - 1) ** 0.5 * \chi
624 Parameters
625 ----------
626 val: (float, array-like)
627 The value to rescale
628 required_variables: dict
629 Dictionary containing the aligned spin component
631 Returns
632 -------
633 (float, array-like)
634 The in-plane component of the spin
635 """
636 self.update_conditions(**required_variables)
637 chi_aligned = abs(required_variables[self._required_variables[0]])
638 return chi_aligned * ((self._reference_maximum / chi_aligned) ** (2 * val) - 1) ** 0.5
640 def _condition_function(self, reference_params, **kwargs):
641 with np.errstate(invalid="ignore"):
642 maximum = np.sqrt(
643 self._reference_maximum ** 2 - kwargs[self._required_variables[0]] ** 2
644 )
645 return dict(minimum=0, maximum=maximum)
647 def __repr__(self):
648 return Prior.__repr__(self)
650 def get_instantiation_dict(self):
651 instantiation_dict = Prior.get_instantiation_dict(self)
652 for key, value in self.reference_params.items():
653 if key in instantiation_dict:
654 instantiation_dict[key] = value
655 return instantiation_dict
658class EOSCheck(Constraint):
659 def __init__(self, minimum=-np.inf, maximum=np.inf):
660 """
661 Constraint used for EoS sampling. Converts the result of various
662 checks on the EoS and its parameters into a prior that can reject
663 unphysical samples. Necessary for EoS sampling.
664 """
666 super().__init__(minimum=minimum, maximum=maximum, name=None, latex_label=None, unit=None)
668 def prob(self, val):
669 """
670 Returns the result of the equation of state check in the conversion function.
671 """
672 return val
674 def ln_prob(self, val):
676 if val:
677 result = 0.0
678 elif not val:
679 result = -np.inf
681 return result
684class CBCPriorDict(ConditionalPriorDict):
685 @property
686 def minimum_chirp_mass(self):
687 mass_1 = None
688 mass_2 = None
689 if "chirp_mass" in self:
690 return self["chirp_mass"].minimum
691 elif "mass_1" in self:
692 mass_1 = self['mass_1'].minimum
693 if "mass_2" in self:
694 mass_2 = self['mass_2'].minimum
695 elif "mass_ratio" in self:
696 mass_2 = mass_1 * self["mass_ratio"].minimum
697 if mass_1 is not None and mass_2 is not None:
698 s = generate_mass_parameters(dict(mass_1=mass_1, mass_2=mass_2))
699 return s["chirp_mass"]
700 else:
701 logger.warning("Unable to determine minimum chirp mass")
702 return None
704 @property
705 def maximum_chirp_mass(self):
706 mass_1 = None
707 mass_2 = None
708 if "chirp_mass" in self:
709 return self["chirp_mass"].maximum
710 elif "mass_1" in self:
711 mass_1 = self['mass_1'].maximum
712 if "mass_2" in self:
713 mass_2 = self['mass_2'].maximum
714 elif "mass_ratio" in self:
715 mass_2 = mass_1 * self["mass_ratio"].maximum
716 if mass_1 is not None and mass_2 is not None:
717 s = generate_mass_parameters(dict(mass_1=mass_1, mass_2=mass_2))
718 return s["chirp_mass"]
719 else:
720 logger.warning("Unable to determine maximum chirp mass")
721 return None
723 @property
724 def minimum_component_mass(self):
725 """
726 The minimum component mass allowed for the prior dictionary.
728 This property requires either:
729 * a prior for :code:`mass_2`
730 * priors for :code:`chirp_mass` and :code:`mass_ratio`
732 Returns
733 -------
734 mass_2: float
735 The minimum allowed component mass.
736 """
737 if "mass_2" in self:
738 return self["mass_2"].minimum
739 if "chirp_mass" in self and "mass_ratio" in self:
740 total_mass = chirp_mass_and_mass_ratio_to_total_mass(
741 self["chirp_mass"].minimum, self["mass_ratio"].minimum)
742 _, mass_2 = total_mass_and_mass_ratio_to_component_masses(
743 self["mass_ratio"].minimum, total_mass)
744 return mass_2
745 else:
746 logger.warning("Unable to determine minimum component mass")
747 return None
749 def is_nonempty_intersection(self, pset):
750 """ Check if keys in self exist in the parameter set
752 Parameters
753 ----------
754 pset: str, set
755 Either a string referencing a parameter set in PARAMETER_SETS or
756 a set of keys
757 """
758 if isinstance(pset, str) and pset in PARAMETER_SETS:
759 check_set = PARAMETER_SETS[pset]
760 elif isinstance(pset, set):
761 check_set = pset
762 else:
763 raise ValueError(f"pset {pset} not understood")
764 if len(check_set.intersection(self.non_fixed_keys)) > 0:
765 return True
766 else:
767 return False
769 @property
770 def spin(self):
771 """ Return true if priors include any spin parameters """
772 return self.is_nonempty_intersection("spin")
774 @property
775 def precession(self):
776 """ Return true if priors include any precession parameters """
777 return self.is_nonempty_intersection("precession_only")
779 @property
780 def measured_spin(self):
781 """ Return true if priors include any measured_spin parameters """
782 return self.is_nonempty_intersection("measured_spin")
784 @property
785 def intrinsic(self):
786 """ Return true if priors include any intrinsic parameters """
787 return self.is_nonempty_intersection("intrinsic")
789 @property
790 def extrinsic(self):
791 """ Return true if priors include any extrinsic parameters """
792 return self.is_nonempty_intersection("extrinsic")
794 @property
795 def sky(self):
796 """ Return true if priors include any extrinsic parameters """
797 return self.is_nonempty_intersection("sky")
799 @property
800 def distance_inclination(self):
801 """ Return true if priors include any extrinsic parameters """
802 return self.is_nonempty_intersection("distance_inclination")
804 @property
805 def mass(self):
806 """ Return true if priors include any mass parameters """
807 return self.is_nonempty_intersection("mass")
809 @property
810 def phase(self):
811 """ Return true if priors include phase parameters """
812 return self.is_nonempty_intersection("phase")
814 def validate_prior(self, duration, minimum_frequency, N=1000, error=True, warning=False):
815 """ Validate the prior is suitable for use
817 Parameters
818 ==========
819 duration: float
820 The data duration in seconds
821 minimum_frequency: float
822 The minimum frequency in Hz of the analysis
823 N: int
824 The number of samples to draw when checking
825 error: bool
826 Whether to raise a ValueError on failure.
827 warning: bool
828 Whether to log a warning on failure.
830 Returns
831 =======
832 bool: whether the template will fit within the segment duration
833 """
834 samples = self.sample(N)
835 samples = generate_all_bbh_parameters(samples)
836 durations = np.array([
837 calculate_time_to_merger(
838 frequency=minimum_frequency,
839 mass_1=mass_1,
840 mass_2=mass_2,
841 )
842 for (mass_1, mass_2) in zip(samples["mass_1"], samples["mass_2"])
843 ])
844 longest_duration = max(durations)
845 if longest_duration < duration:
846 return True
847 message = (
848 "Prior contains regions of parameter space in which the signal"
849 f" is longer than the data duration {duration}s."
850 f" Maximum estimated duration = {longest_duration:.1f}s."
851 )
852 if warning:
853 logger.warning(message)
854 return False
855 if error:
856 raise ValueError(message)
859class BBHPriorDict(CBCPriorDict):
860 def __init__(self, dictionary=None, filename=None, aligned_spin=False,
861 conversion_function=None):
862 """ Initialises a Prior set for Binary Black holes
864 Parameters
865 ==========
866 dictionary: dict, optional
867 See superclass
868 filename: str, optional
869 See superclass
870 conversion_function: func
871 Function to convert between sampled parameters and constraints.
872 By default this generates many additional parameters, see
873 BBHPriorDict.default_conversion_function
874 """
875 if dictionary is None and filename is None:
876 if aligned_spin:
877 fname = 'aligned_spins_bbh.prior'
878 logger.info('Using aligned spin prior')
879 else:
880 fname = 'precessing_spins_bbh.prior'
881 filename = os.path.join(DEFAULT_PRIOR_DIR, fname)
882 logger.info('No prior given, using default BBH priors in {}.'.format(filename))
883 elif filename is not None:
884 if not os.path.isfile(filename):
885 filename = os.path.join(DEFAULT_PRIOR_DIR, filename)
886 super(BBHPriorDict, self).__init__(dictionary=dictionary, filename=filename,
887 conversion_function=conversion_function)
889 def default_conversion_function(self, sample):
890 """
891 Default parameter conversion function for BBH signals.
893 This generates:
894 - the parameters passed to source.lal_binary_black_hole
895 - all mass parameters
897 It does not generate:
898 - component spins
899 - source-frame parameters
901 Parameters
902 ==========
903 sample: dict
904 Dictionary to convert
906 Returns
907 =======
908 sample: dict
909 Same as input
910 """
911 out_sample = fill_from_fixed_priors(sample, self)
912 out_sample, _ = convert_to_lal_binary_black_hole_parameters(out_sample)
913 out_sample = generate_mass_parameters(out_sample)
915 return out_sample
917 def test_redundancy(self, key, disable_logging=False):
918 """
919 Test whether adding the key would add be redundant.
920 Already existing keys return True.
922 Parameters
923 ==========
924 key: str
925 The key to test.
926 disable_logging: bool, optional
927 Disable logging in this function call. Default is False.
929 Returns
930 =======
931 redundant: bool
932 Whether the key is redundant or not
933 """
934 if key in self:
935 logger.debug('{} already in prior'.format(key))
936 return True
938 sampling_parameters = {key for key in self if not isinstance(
939 self[key], (DeltaFunction, Constraint))}
941 mass_parameters = {'mass_1', 'mass_2', 'chirp_mass', 'total_mass', 'mass_ratio', 'symmetric_mass_ratio'}
942 spin_tilt_1_parameters = {'tilt_1', 'cos_tilt_1'}
943 spin_tilt_2_parameters = {'tilt_2', 'cos_tilt_2'}
944 spin_azimuth_parameters = {'phi_1', 'phi_2', 'phi_12', 'phi_jl'}
945 inclination_parameters = {'theta_jn', 'cos_theta_jn'}
946 distance_parameters = {'luminosity_distance', 'comoving_distance', 'redshift'}
948 for independent_parameters, parameter_set in \
949 zip([2, 2, 1, 1, 1, 1],
950 [mass_parameters, spin_azimuth_parameters,
951 spin_tilt_1_parameters, spin_tilt_2_parameters,
952 inclination_parameters, distance_parameters]):
953 if key in parameter_set:
954 if len(parameter_set.intersection(
955 sampling_parameters)) >= independent_parameters:
956 logger.disabled = disable_logging
957 logger.warning('{} already in prior. '
958 'This may lead to unexpected behaviour.'
959 .format(parameter_set.intersection(self)))
960 logger.disabled = False
961 return True
962 return False
965class BNSPriorDict(CBCPriorDict):
967 def __init__(self, dictionary=None, filename=None, aligned_spin=True,
968 conversion_function=None):
969 """ Initialises a Prior set for Binary Neutron Stars
971 Parameters
972 ==========
973 dictionary: dict, optional
974 See superclass
975 filename: str, optional
976 See superclass
977 conversion_function: func
978 Function to convert between sampled parameters and constraints.
979 By default this generates many additional parameters, see
980 BNSPriorDict.default_conversion_function
981 """
982 if aligned_spin:
983 default_file = 'aligned_spins_bns_tides_on.prior'
984 else:
985 default_file = 'precessing_spins_bns_tides_on.prior'
986 if dictionary is None and filename is None:
987 filename = os.path.join(DEFAULT_PRIOR_DIR, default_file)
988 logger.info('No prior given, using default BNS priors in {}.'.format(filename))
989 elif filename is not None:
990 if not os.path.isfile(filename):
991 filename = os.path.join(DEFAULT_PRIOR_DIR, filename)
992 super(BNSPriorDict, self).__init__(dictionary=dictionary, filename=filename,
993 conversion_function=conversion_function)
995 def default_conversion_function(self, sample):
996 """
997 Default parameter conversion function for BNS signals.
999 This generates:
1001 * the parameters passed to source.lal_binary_neutron_star
1002 * all mass parameters
1003 * all tidal parameters
1005 It does not generate:
1007 * component spins
1008 * source-frame parameters
1010 Parameters
1011 ==========
1012 sample: dict
1013 Dictionary to convert
1015 Returns
1016 =======
1017 sample: dict
1018 Same as input
1019 """
1020 out_sample = fill_from_fixed_priors(sample, self)
1021 out_sample, _ = convert_to_lal_binary_neutron_star_parameters(out_sample)
1022 out_sample = generate_mass_parameters(out_sample)
1023 out_sample = generate_tidal_parameters(out_sample)
1024 return out_sample
1026 def test_redundancy(self, key, disable_logging=False):
1027 logger.disabled = disable_logging
1028 logger.info("Performing redundancy check using BBHPriorDict(self).test_redundancy")
1029 bbh_redundancy = BBHPriorDict(self).test_redundancy(key)
1030 logger.disabled = False
1032 if bbh_redundancy:
1033 return True
1034 redundant = False
1036 sampling_parameters = {key for key in self if not isinstance(
1037 self[key], (DeltaFunction, Constraint))}
1039 tidal_parameters = \
1040 {'lambda_1', 'lambda_2', 'lambda_tilde', 'delta_lambda_tilde'}
1042 if key in tidal_parameters:
1043 if len(tidal_parameters.intersection(sampling_parameters)) > 2:
1044 redundant = True
1045 logger.disabled = disable_logging
1046 logger.warning('{} already in prior. '
1047 'This may lead to unexpected behaviour.'
1048 .format(tidal_parameters.intersection(self)))
1049 logger.disabled = False
1050 elif len(tidal_parameters.intersection(sampling_parameters)) == 2:
1051 redundant = True
1052 return redundant
1054 @property
1055 def tidal(self):
1056 """ Return true if priors include phase parameters """
1057 return self.is_nonempty_intersection("tidal")
1060Prior._default_latex_labels = {
1061 'mass_1': '$m_1$',
1062 'mass_2': '$m_2$',
1063 'total_mass': '$M$',
1064 'chirp_mass': r'$\mathcal{M}$',
1065 'mass_ratio': '$q$',
1066 'symmetric_mass_ratio': r'$\eta$',
1067 'a_1': '$a_1$',
1068 'a_2': '$a_2$',
1069 'tilt_1': r'$\theta_1$',
1070 'tilt_2': r'$\theta_2$',
1071 'cos_tilt_1': r'$\cos\theta_1$',
1072 'cos_tilt_2': r'$\cos\theta_2$',
1073 'phi_12': r'$\Delta\phi$',
1074 'phi_jl': r'$\phi_{JL}$',
1075 'luminosity_distance': '$d_L$',
1076 'dec': r'$\mathrm{DEC}$',
1077 'ra': r'$\mathrm{RA}$',
1078 'iota': r'$\iota$',
1079 'cos_iota': r'$\cos\iota$',
1080 'theta_jn': r'$\theta_{JN}$',
1081 'cos_theta_jn': r'$\cos\theta_{JN}$',
1082 'psi': r'$\psi$',
1083 'phase': r'$\phi$',
1084 'geocent_time': '$t_c$',
1085 'time_jitter': '$t_j$',
1086 'lambda_1': r'$\Lambda_1$',
1087 'lambda_2': r'$\Lambda_2$',
1088 'lambda_tilde': r'$\tilde{\Lambda}$',
1089 'delta_lambda_tilde': r'$\delta\tilde{\Lambda}$',
1090 'chi_1': r'$\chi_1$',
1091 'chi_2': r'$\chi_2$',
1092 'chi_1_in_plane': r'$\chi_{1, \perp}$',
1093 'chi_2_in_plane': r'$\chi_{2, \perp}$',
1094}
1097class CalibrationPriorDict(PriorDict):
1098 """
1099 Prior dictionary class for calibration parameters.
1100 This has methods for simplifying the creation of priors for the large
1101 numbers of parameters used with the spline model.
1102 """
1104 def __init__(self, dictionary=None, filename=None):
1105 """
1106 Initialises a Prior dictionary for calibration parameters
1108 Parameters
1109 ==========
1110 dictionary: dict, optional
1111 See superclass
1112 filename: str, optional
1113 See superclass
1114 """
1115 if dictionary is None and filename is not None:
1116 filename = os.path.join(DEFAULT_PRIOR_DIR, filename)
1117 super(CalibrationPriorDict, self).__init__(dictionary=dictionary, filename=filename)
1118 self.source = None
1120 def to_file(self, outdir, label):
1121 """
1122 Write the prior to file. This includes information about the source if
1123 possible.
1125 Parameters
1126 ==========
1127 outdir: str
1128 Output directory.
1129 label: str
1130 Label for prior.
1131 """
1132 PriorDict.to_file(self, outdir=outdir, label=label)
1133 if self.source is not None:
1134 prior_file = os.path.join(outdir, "{}.prior".format(label))
1135 with open(prior_file, "a") as outfile:
1136 outfile.write("# prior source file is {}".format(self.source))
1138 @staticmethod
1139 def from_envelope_file(envelope_file, minimum_frequency,
1140 maximum_frequency, n_nodes, label,
1141 boundary="reflective"):
1142 """
1143 Load in the calibration envelope.
1145 This is a text file with columns
1147 ::
1149 frequency median-amplitude median-phase -1-sigma-amplitude
1150 -1-sigma-phase +1-sigma-amplitude +1-sigma-phase
1152 Parameters
1153 ==========
1154 envelope_file: str
1155 Name of file to read in.
1156 minimum_frequency: float
1157 Minimum frequency for the spline.
1158 maximum_frequency: float
1159 Minimum frequency for the spline.
1160 n_nodes: int
1161 Number of nodes for the spline.
1162 label: str
1163 Label for the names of the parameters, e.g., `recalib_H1_`
1164 boundary: None, 'reflective', 'periodic'
1165 The type of prior boundary to assign
1167 Returns
1168 =======
1169 prior: PriorDict
1170 Priors for the relevant parameters.
1171 This includes the frequencies of the nodes which are _not_ sampled.
1172 """
1173 calibration_data = np.genfromtxt(envelope_file).T
1174 log_frequency_array = np.log(calibration_data[0])
1175 amplitude_median = calibration_data[1] - 1
1176 phase_median = calibration_data[2]
1177 amplitude_sigma = (calibration_data[5] - calibration_data[3]) / 2
1178 phase_sigma = (calibration_data[6] - calibration_data[4]) / 2
1180 log_nodes = np.linspace(np.log(minimum_frequency),
1181 np.log(maximum_frequency), n_nodes)
1183 amplitude_mean_nodes = \
1184 InterpolatedUnivariateSpline(log_frequency_array, amplitude_median)(log_nodes)
1185 amplitude_sigma_nodes = \
1186 InterpolatedUnivariateSpline(log_frequency_array, amplitude_sigma)(log_nodes)
1187 phase_mean_nodes = \
1188 InterpolatedUnivariateSpline(log_frequency_array, phase_median)(log_nodes)
1189 phase_sigma_nodes = \
1190 InterpolatedUnivariateSpline(log_frequency_array, phase_sigma)(log_nodes)
1192 prior = CalibrationPriorDict()
1193 for ii in range(n_nodes):
1194 name = "recalib_{}_amplitude_{}".format(label, ii)
1195 latex_label = "$A^{}_{}$".format(label, ii)
1196 prior[name] = Gaussian(mu=amplitude_mean_nodes[ii],
1197 sigma=amplitude_sigma_nodes[ii],
1198 name=name, latex_label=latex_label,
1199 boundary=boundary)
1200 for ii in range(n_nodes):
1201 name = "recalib_{}_phase_{}".format(label, ii)
1202 latex_label = r"$\phi^{}_{}$".format(label, ii)
1203 prior[name] = Gaussian(mu=phase_mean_nodes[ii],
1204 sigma=phase_sigma_nodes[ii],
1205 name=name, latex_label=latex_label,
1206 boundary=boundary)
1207 for ii in range(n_nodes):
1208 name = "recalib_{}_frequency_{}".format(label, ii)
1209 latex_label = "$f^{}_{}$".format(label, ii)
1210 prior[name] = DeltaFunction(peak=np.exp(log_nodes[ii]), name=name,
1211 latex_label=latex_label)
1212 prior.source = os.path.abspath(envelope_file)
1213 return prior
1215 @staticmethod
1216 def constant_uncertainty_spline(
1217 amplitude_sigma, phase_sigma, minimum_frequency, maximum_frequency,
1218 n_nodes, label, boundary="reflective"):
1219 """
1220 Make prior assuming constant in frequency calibration uncertainty.
1222 This assumes Gaussian fluctuations about 0.
1224 Parameters
1225 ==========
1226 amplitude_sigma: float
1227 Uncertainty in the amplitude.
1228 phase_sigma: float
1229 Uncertainty in the phase.
1230 minimum_frequency: float
1231 Minimum frequency for the spline.
1232 maximum_frequency: float
1233 Minimum frequency for the spline.
1234 n_nodes: int
1235 Number of nodes for the spline.
1236 label: str
1237 Label for the names of the parameters, e.g., `recalib_H1_`
1238 boundary: None, 'reflective', 'periodic'
1239 The type of prior boundary to assign
1241 Returns
1242 =======
1243 prior: PriorDict
1244 Priors for the relevant parameters.
1245 This includes the frequencies of the nodes which are _not_ sampled.
1246 """
1247 nodes = np.logspace(np.log10(minimum_frequency),
1248 np.log10(maximum_frequency), n_nodes)
1250 amplitude_mean_nodes = [0] * n_nodes
1251 amplitude_sigma_nodes = [amplitude_sigma] * n_nodes
1252 phase_mean_nodes = [0] * n_nodes
1253 phase_sigma_nodes = [phase_sigma] * n_nodes
1255 prior = CalibrationPriorDict()
1256 for ii in range(n_nodes):
1257 name = "recalib_{}_amplitude_{}".format(label, ii)
1258 latex_label = "$A^{}_{}$".format(label, ii)
1259 prior[name] = Gaussian(mu=amplitude_mean_nodes[ii],
1260 sigma=amplitude_sigma_nodes[ii],
1261 name=name, latex_label=latex_label,
1262 boundary=boundary)
1263 for ii in range(n_nodes):
1264 name = "recalib_{}_phase_{}".format(label, ii)
1265 latex_label = r"$\phi^{}_{}$".format(label, ii)
1266 prior[name] = Gaussian(mu=phase_mean_nodes[ii],
1267 sigma=phase_sigma_nodes[ii],
1268 name=name, latex_label=latex_label,
1269 boundary=boundary)
1270 for ii in range(n_nodes):
1271 name = "recalib_{}_frequency_{}".format(label, ii)
1272 latex_label = "$f^{}_{}$".format(label, ii)
1273 prior[name] = DeltaFunction(peak=nodes[ii], name=name,
1274 latex_label=latex_label)
1276 return prior
1279def secondary_mass_condition_function(reference_params, mass_1):
1280 """
1281 Condition function to use for a prior that is conditional on the value
1282 of the primary mass, for example, a prior on the secondary mass that is
1283 bounded by the primary mass.
1285 .. code-block:: python
1287 import bilby
1288 priors = bilby.gw.prior.CBCPriorDict()
1289 priors["mass_1"] = bilby.core.prior.Uniform(5, 50)
1290 priors["mass_2"] = bilby.core.prior.ConditionalUniform(
1291 condition_func=bilby.gw.prior.secondary_mass_condition_function,
1292 minimum=5,
1293 maximum=50,
1294 )
1296 Parameters
1297 ----------
1298 reference_params: dict
1299 Reference parameter dictionary, this should contain a "minimum"
1300 attribute.
1301 mass_1: float
1302 The primary mass to use as the upper limit for the prior.
1304 Returns
1305 -------
1306 dict:
1307 Updated prior limits given the provided primary mass.
1308 """
1309 return dict(minimum=reference_params['minimum'], maximum=mass_1)
1312ConditionalCosmological = conditional_prior_factory(Cosmological)
1313ConditionalUniformComovingVolume = conditional_prior_factory(UniformComovingVolume)
1314ConditionalUniformSourceFrame = conditional_prior_factory(UniformSourceFrame)
1317class HealPixMapPriorDist(BaseJointPriorDist):
1318 """
1319 Class defining prior according to given HealPix Map, defaults to 2D in ra and dec but can be set to include
1320 Distance as well. This only works with skymaps that include the 2D joint probability in ra/dec and that use the
1321 normal LALInference type skymaps where each pixel has a DISTMU, DISTSIGMA, and DISTNORM defining the conditional
1322 distance distribution along a given line of sight.
1324 Parameters
1325 ==========
1327 hp_file : file path to .fits file
1328 .fits file that contains the 2D or 3D Healpix Map
1329 names : list (optional)
1330 list of names of parameters included in the JointPriorDist, defaults to ['ra', 'dec']
1331 bounds : dict or list (optional)
1332 dictionary or list with given prior bounds. defaults to normal bounds on ra, dev and 0, inf for distance
1333 if this is for a 3D map
1335 Returns
1336 =======
1338 PriorDist : `bilby.gw.prior.HealPixMapPriorDist`
1339 A JointPriorDist object to store the joint prior distribution according to passed healpix map
1340 """
1341 def __init__(self, hp_file, names=None, bounds=None, distance=False):
1342 self.hp = self._check_imports()
1343 self.hp_file = hp_file
1344 if names is None:
1345 names = ["ra", "dec"]
1346 if bounds is None:
1347 bounds = [[0, 2 * np.pi], [-np.pi / 2.0, np.pi / 2.0]]
1348 elif isinstance(bounds, dict):
1349 bs = [[] for _ in bounds.keys()]
1350 for i, key in enumerate(bounds.keys()):
1351 bs[i] = (bounds[key][0], bounds[key][1])
1352 bounds = bs
1353 if distance:
1354 if len(names) == 2:
1355 names.append("distance")
1356 if len(bounds) == 2:
1357 bounds.append([0, np.inf])
1358 self.distance = True
1359 self.prob, self.distmu, self.distsigma, self.distnorm = self.hp.read_map(
1360 hp_file, field=range(4)
1361 )
1362 else:
1363 self.distance = False
1364 self.prob = self.hp.read_map(hp_file)
1365 self.prob = self._check_norm(self.prob)
1367 super(HealPixMapPriorDist, self).__init__(names=names, bounds=bounds)
1368 self.distname = "hpmap"
1369 self.npix = len(self.prob)
1370 self.nside = self.hp.npix2nside(self.npix)
1371 self.pixel_area = self.hp.nside2pixarea(self.nside)
1372 self.pixel_length = self.pixel_area ** (1 / 2.0)
1373 self.pix_xx = np.arange(self.npix)
1374 self._all_interped = interp1d(x=self.pix_xx, y=self.prob, bounds_error=False, fill_value=0)
1375 self.inverse_cdf = None
1376 self.distance_pdf = None
1377 self.distance_icdf = None
1378 self._build_attributes()
1379 name = self.names[-1]
1380 if self.bounds[name][1] != np.inf and self.bounds[name][0] != -np.inf:
1381 self.rs = np.linspace(self.bounds[name][0], self.bounds[name][1], 1000)
1382 else:
1383 self.rs = np.linspace(0, 5000, 1000)
1385 def _build_attributes(self):
1386 """
1387 Method that builds the inverse cdf of the P(pixel) distribution for rescaling
1388 """
1389 from scipy.integrate import cumulative_trapezoid
1390 yy = self._all_interped(self.pix_xx)
1391 yy /= np.trapz(yy, self.pix_xx)
1392 YY = cumulative_trapezoid(yy, self.pix_xx, initial=0)
1393 YY[-1] = 1
1394 self.inverse_cdf = interp1d(x=YY, y=self.pix_xx, bounds_error=True)
1396 @staticmethod
1397 def _check_imports():
1398 """
1399 Static method to check that healpy is installed on the machine running bibly
1400 """
1401 try:
1402 import healpy
1403 except Exception:
1404 raise ImportError("Must have healpy installed on this machine to use HealPixMapPrior")
1405 return healpy
1407 def _rescale(self, samp, **kwargs):
1408 """
1409 Overwrites the _rescale method of BaseJoint Prior to rescale a single value from the unitcube onto
1410 two values (ra, dec) or 3 (ra, dec, dist) if distance is included
1412 Parameters
1413 ==========
1414 samp : float, int
1415 must take in single value for pixel on unitcube to recale onto ra, dec (distance), for the map Prior
1416 kwargs : dict
1417 kwargs are all passed to _rescale() method
1419 Returns
1420 =======
1421 rescaled_sample : array_like
1422 sample to rescale onto the prior
1423 """
1424 if self.distance:
1425 dist_samp = samp[:, -1]
1426 samp = samp[:, 0]
1427 else:
1428 samp = samp[:, 0]
1429 pix_rescale = self.inverse_cdf(samp)
1430 sample = np.empty((len(pix_rescale), 2))
1431 dist_samples = np.empty((len(pix_rescale)))
1432 for i, val in enumerate(pix_rescale):
1433 theta, ra = self.hp.pix2ang(self.nside, int(round(val)))
1434 dec = 0.5 * np.pi - theta
1435 sample[i, :] = self.draw_from_pixel(ra, dec, int(round(val)))
1436 if self.distance:
1437 self.update_distance(int(round(val)))
1438 dist_samples[i] = self.distance_icdf(dist_samp[i])
1439 if self.distance:
1440 sample = np.vstack([sample[:, 0], sample[:, 1], dist_samples])
1441 return sample.reshape((-1, self.num_vars))
1443 def update_distance(self, pix_idx):
1444 """
1445 Method to update the conditional distance distributions at given pixel used for distance handling in the
1446 JointPrior Parameters. This function updates the current distance pdf, inverse_cdf, and sampler according to
1447 given pixel or line of sight.
1449 Parameters
1450 ==========
1451 pix_idx : int
1452 pixel index value to create the distribution for
1454 Returns
1455 =======
1456 None : None
1457 just updates these functions at new pixel values
1458 """
1459 self.distance_pdf = lambda r: self.distnorm[pix_idx] * norm(
1460 loc=self.distmu[pix_idx], scale=self.distsigma[pix_idx]
1461 ).pdf(r)
1462 pdfs = self.rs ** 2 * self.distance_pdf(self.rs)
1463 cdfs = np.cumsum(pdfs) / np.sum(pdfs)
1464 self.distance_icdf = interp1d(cdfs, self.rs)
1466 @staticmethod
1467 def _check_norm(array):
1468 """
1469 static method to check if array is properlly normalized and if not to normalize it.
1471 Parameters
1472 ==========
1473 array : array_like
1474 input array we want to renormalize if not already normalized
1476 Returns
1477 =======
1478 normed_array : array_like
1479 returns input array normalized
1480 """
1481 norm = np.linalg.norm(array, ord=1)
1482 if norm == 0:
1483 norm = np.finfo(array.dtype).eps
1484 return array / norm
1486 def _sample(self, size, **kwargs):
1487 """
1488 Overwrites the _sample method of BaseJoint Prior. Picks a pixel value according to their probabilities, then
1489 uniformly samples ra, and decs that are contained in chosen pixel. If the PriorDist includes distance it then
1490 updates the distance distributions and will sample according to the conditional distance distribution along a
1491 given line of sight
1493 Parameters
1494 ==========
1495 size : int
1496 number of samples we want to draw
1497 kwargs : dict
1498 kwargs are all passed to be used
1500 Returns
1501 =======
1502 sample : array_like
1503 sample of ra, and dec (and distance if 3D=True)
1504 """
1505 sample_pix = random.rng.choice(self.npix, size=size, p=self.prob, replace=True)
1506 sample = np.empty((size, self.num_vars))
1507 for samp in range(size):
1508 theta, ra = self.hp.pix2ang(self.nside, sample_pix[samp])
1509 dec = 0.5 * np.pi - theta
1510 if self.distance:
1511 self.update_distance(sample_pix[samp])
1512 dist = self.draw_distance(sample_pix[samp])
1513 ra_dec = self.draw_from_pixel(ra, dec, sample_pix[samp])
1514 sample[samp, :] = [ra_dec[0], ra_dec[1], dist]
1515 else:
1516 sample[samp, :] = self.draw_from_pixel(ra, dec, sample_pix[samp])
1517 return sample.reshape((-1, self.num_vars))
1519 def draw_distance(self, pix):
1520 """
1521 Method to recursively draw a distance value from the given set distance distribution and check that it is in
1522 the bounds
1524 Parameters
1525 ==========
1527 pix : int
1528 integer for pixel to draw a distance from
1530 Returns
1531 =======
1532 dist : float
1533 sample drawn from the distance distribution at set pixel index
1534 """
1535 if self.distmu[pix] == np.inf or self.distmu[pix] <= 0:
1536 return 0
1537 dist = self.distance_icdf(random.rng.uniform(0, 1))
1538 name = self.names[-1]
1539 if (dist > self.bounds[name][1]) | (dist < self.bounds[name][0]):
1540 self.draw_distance(pix)
1541 else:
1542 return dist
1544 def draw_from_pixel(self, ra, dec, pix):
1545 """
1546 Recursive function to uniformly draw ra, and dec values that are located in the given pixel
1548 Parameters
1549 ==========
1550 ra : float, int
1551 value drawn for rightascension
1552 dec : float, int
1553 value drawn for declination
1554 pix : int
1555 pixel index for given pixel we want to get ra, and dec from
1557 Returns
1558 =======
1559 ra_dec : tuple
1560 this returns a tuple of ra, and dec sampled uniformly that are in the pixel given
1561 """
1562 if not self.check_in_pixel(ra, dec, pix):
1563 self.draw_from_pixel(ra, dec, pix)
1564 return np.array(
1565 [
1566 random.rng.uniform(ra - self.pixel_length, ra + self.pixel_length),
1567 random.rng.uniform(dec - self.pixel_length, dec + self.pixel_length),
1568 ]
1569 )
1571 def check_in_pixel(self, ra, dec, pix):
1572 """
1573 Method that checks if given rightacension and declination values are within the given pixel index and the bounds
1575 Parameters
1576 ==========
1577 ra : float, int
1578 rightascension value to check
1579 dec : float, int
1580 declination value to check
1581 pix : int
1582 index for pixel we want to check in
1584 Returns
1585 =======
1586 bool :
1587 returns True if values inside pixel, False if not
1588 """
1589 for val, name in zip([ra, dec], self.names):
1590 if (val < self.bounds[name][0]) or (val > self.bounds[name][1]):
1591 return False
1592 phi, theta = ra, 0.5 * np.pi - dec
1593 pixel = self.hp.ang2pix(self.nside, theta, phi)
1594 return pix == pixel
1596 def _ln_prob(self, samp, lnprob, outbounds):
1597 """
1598 Overwrites the _lnprob method of BaseJoint Prior
1600 Parameters
1601 ==========
1602 samp : array_like
1603 samples of ra, dec to evaluate the lnprob at
1604 lnprob : array_like
1605 array of correct length we want to populate with lnprob values
1606 outbounds : boolean array
1607 boolean array that flags samples that are out of the given bounds
1609 Returns
1610 =======
1611 lnprob : array_like
1612 lnprob values at each sample
1613 """
1614 for i in range(samp.shape[0]):
1615 if not outbounds[i]:
1616 if self.distance:
1617 phi, dec, dist = samp[0]
1618 else:
1619 phi, dec = samp[0]
1620 theta = 0.5 * np.pi - dec
1621 pixel = self.hp.ang2pix(self.nside, theta, phi)
1622 lnprob[i] = np.log(self.prob[pixel] / self.pixel_area)
1623 if self.distance:
1624 self.update_distance(pixel)
1625 lnprob[i] += np.log(self.distance_pdf(dist) * dist ** 2)
1626 lnprob[outbounds] = -np.inf
1627 return lnprob
1629 def __eq__(self, other):
1630 skip_keys = ["_all_interped", "inverse_cdf", "distance_pdf", "distance_icdf"]
1631 if self.__class__ != other.__class__:
1632 return False
1633 if sorted(self.__dict__.keys()) != sorted(other.__dict__.keys()):
1634 return False
1635 for key in self.__dict__:
1636 if key in skip_keys:
1637 continue
1638 if key == "hp_file":
1639 if self.__dict__[key] != other.__dict__[key]:
1640 return False
1641 elif isinstance(self.__dict__[key], (np.ndarray, list)):
1642 thisarr = np.asarray(self.__dict__[key])
1643 otherarr = np.asarray(other.__dict__[key])
1644 if thisarr.dtype == float and otherarr.dtype == float:
1645 fin1 = np.isfinite(np.asarray(self.__dict__[key]))
1646 fin2 = np.isfinite(np.asarray(other.__dict__[key]))
1647 if not np.array_equal(fin1, fin2):
1648 return False
1649 if not np.allclose(thisarr[fin1], otherarr[fin2], atol=1e-15):
1650 return False
1651 else:
1652 if not np.array_equal(thisarr, otherarr):
1653 return False
1654 else:
1655 if not self.__dict__[key] == other.__dict__[key]:
1656 return False
1657 return True
1660class HealPixPrior(JointPrior):
1661 """
1662 A prior distribution that follows a user-provided HealPix map for one
1663 parameter.
1665 See :code:`bilby.gw.prior.HealPixMapPriorDist` for more details of how to
1666 instantiate the prior.
1667 """
1668 def __init__(self, dist, name=None, latex_label=None, unit=None):
1669 """
1671 Parameters
1672 ----------
1673 dist: bilby.gw.prior.HealPixMapPriorDist
1674 The base joint probability.
1675 name: str
1676 The name of the parameter, it should be contained in the map.
1677 One of ["ra", "dec", "luminosity_distance"].
1678 latex_label: str
1679 Latex label used for plotting, will be read from default values if
1680 not provided.
1681 unit: str
1682 The unit of the parameter.
1683 """
1684 if not isinstance(dist, HealPixMapPriorDist):
1685 raise JointPriorDistError("dist object must be instance of HealPixMapPriorDist")
1686 super(HealPixPrior, self).__init__(dist=dist, name=name, latex_label=latex_label, unit=unit)