Coverage for pesummary/gw/plots/plot.py: 81.1%
681 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-05-02 08:42 +0000
1# Licensed under an MIT style license -- see LICENSE.md
3from pesummary.utils.utils import (
4 logger, number_of_columns_for_legend, _check_latex_install,
5)
6from pesummary.utils.decorators import no_latex_plot
7from pesummary.gw.plots.bounds import default_bounds
8from pesummary.core.plots.figure import figure, subplots, ExistingFigure
9from pesummary.core.plots.plot import _default_legend_kwargs
10from pesummary import conf
12import os
13import matplotlib.style
14import numpy as np
15import math
16from scipy.ndimage import gaussian_filter
17from astropy.time import Time
19_check_latex_install()
21from lal import MSUN_SI, PC_SI
23__author__ = ["Charlie Hoy <charlie.hoy@ligo.org>"]
24try:
25 import lalsimulation as lalsim
26 LALSIMULATION = True
27except ImportError:
28 LALSIMULATION = None
31def _return_bounds(param, samples, comparison=False):
32 """Return the bounds for a given param
34 Parameters
35 ----------
36 param: str
37 name of the parameter you wish to get bounds for
38 samples: list/np.ndarray
39 array or list of array of posterior samples for param
40 comparison: Bool, optional
41 True if samples is a list of array's of posterior samples
42 """
43 xlow, xhigh = None, None
44 if param in default_bounds.keys():
45 bounds = default_bounds[param]
46 if "low" in bounds.keys():
47 xlow = bounds["low"]
48 if "high" in bounds.keys():
49 if isinstance(bounds["high"], str) and "mass_1" in bounds["high"]:
50 if comparison:
51 xhigh = np.max([np.max(i) for i in samples])
52 else:
53 xhigh = np.max(samples)
54 else:
55 xhigh = bounds["high"]
56 return xlow, xhigh
59def _add_default_bounds_to_kde_kwargs_dict(
60 kde_kwargs, param, samples, comparison=False
61):
62 """Add default kde bounds to the a dictionary of kwargs
64 Parameters
65 ----------
66 kde_kwargs: dict
67 dictionary of kwargs to pass to the kde class
68 param: str
69 name of the parameter you wish to plot
70 samples: list
71 list of samples for param
72 """
73 from pesummary.utils.bounded_1d_kde import bounded_1d_kde
75 xlow, xhigh = _return_bounds(param, samples, comparison=comparison)
76 kde_kwargs["xlow"] = xlow
77 kde_kwargs["xhigh"] = xhigh
78 kde_kwargs["kde_kernel"] = bounded_1d_kde
79 return kde_kwargs
82def _1d_histogram_plot(
83 param, samples, *args, kde_kwargs={}, bounded=True, **kwargs
84):
85 """Generate the 1d histogram plot for a given parameter for a given
86 approximant.
88 Parameters
89 ----------
90 *args: tuple
91 all args passed directly to pesummary.core.plots.plot._1d_histogram_plot
92 function
93 kde_kwargs: dict, optional
94 optional kwargs passed to the kde class
95 bounded: Bool, optional
96 if True, pass default 'xlow' and 'xhigh' arguments to the kde class
97 **kwargs: dict, optional
98 all additional kwargs passed to the
99 pesummary.core.plots.plot._1d_histogram_plot function
100 """
101 from pesummary.core.plots.plot import _1d_histogram_plot
103 if bounded:
104 kde_kwargs = _add_default_bounds_to_kde_kwargs_dict(
105 kde_kwargs, param, samples
106 )
107 return _1d_histogram_plot(
108 param, samples, *args, kde_kwargs=kde_kwargs, **kwargs
109 )
112def _1d_histogram_plot_mcmc(
113 param, samples, *args, kde_kwargs={}, bounded=True, **kwargs
114):
115 """Generate the 1d histogram plot for a given parameter for set of
116 mcmc chains
118 Parameters
119 ----------
120 *args: tuple
121 all args passed directly to
122 pesummary.core.plots.plot._1d_histogram_plot_mcmc function
123 kde_kwargs: dict, optional
124 optional kwargs passed to the kde class
125 bounded: Bool, optional
126 if True, pass default 'xlow' and 'xhigh' arguments to the kde class
127 **kwargs: dict, optional
128 all additional kwargs passed to the
129 pesummary.core.plots.plot._1d_histogram_plot_mcmc function
130 """
131 from pesummary.core.plots.plot import _1d_histogram_plot_mcmc
133 if bounded:
134 kde_kwargs = _add_default_bounds_to_kde_kwargs_dict(
135 kde_kwargs, param, samples, comparison=True
136 )
137 return _1d_histogram_plot_mcmc(
138 param, samples, *args, kde_kwargs=kde_kwargs, **kwargs
139 )
142def _1d_histogram_plot_bootstrap(
143 param, samples, *args, kde_kwargs={}, bounded=True, **kwargs
144):
145 """Generate a bootstrapped 1d histogram plot for a given parameter
147 Parameters
148 ----------
149 param: str
150 name of the parameter that you wish to plot
151 samples: np.ndarray
152 array of samples for param
153 args: tuple
154 all args passed to
155 pesummary.core.plots.plot._1d_histogram_plot_bootstrap function
156 kde_kwargs: dict, optional
157 optional kwargs passed to the kde class
158 bounded: Bool, optional
159 if True, pass default 'xlow' and 'xhigh' arguments to the kde class
160 **kwargs: dict, optional
161 all additional kwargs passed to the
162 pesummary.core.plots.plot._1d_histogram_plot_bootstrap function
163 """
164 from pesummary.core.plots.plot import _1d_histogram_plot_bootstrap
166 if bounded:
167 kde_kwargs = _add_default_bounds_to_kde_kwargs_dict(
168 kde_kwargs, param, samples
169 )
170 return _1d_histogram_plot_bootstrap(
171 param, samples, *args, kde_kwargs=kde_kwargs, **kwargs
172 )
175def _1d_comparison_histogram_plot(
176 param, samples, *args, kde_kwargs={}, bounded=True, max_vline=2,
177 legend_kwargs=_default_legend_kwargs, **kwargs
178):
179 """Generate the a plot to compare the 1d_histogram plots for a given
180 parameter for different approximants.
182 Parameters
183 ----------
184 *args: tuple
185 all args passed directly to
186 pesummary.core.plots.plot._1d_comparisonhistogram_plot function
187 kde_kwargs: dict, optional
188 optional kwargs passed to the kde class
189 bounded: Bool, optional
190 if True, pass default 'xlow' and 'xhigh' arguments to the kde class
191 max_vline: int, optional
192 if number of peaks < max_vline draw peaks as vertical lines rather
193 than histogramming the data
194 **kwargs: dict, optional
195 all additional kwargs passed to the
196 pesummary.core.plots.plot._1d_comparison_histogram_plot function
197 """
198 from pesummary.core.plots.plot import _1d_comparison_histogram_plot
200 if bounded:
201 kde_kwargs = _add_default_bounds_to_kde_kwargs_dict(
202 kde_kwargs, param, samples, comparison=True
203 )
204 return _1d_comparison_histogram_plot(
205 param, samples, *args, kde_kwargs=kde_kwargs, max_vline=max_vline,
206 legend_kwargs=legend_kwargs, **kwargs
207 )
210def _make_corner_plot(samples, latex_labels, corner_parameters=None, **kwargs):
211 """Generate the corner plots for a given approximant
213 Parameters
214 ----------
215 opts: argparse
216 argument parser object to hold all information from the command line
217 samples: nd list
218 nd list of samples for each parameter for a given approximant
219 params: list
220 list of parameters associated with each element in samples
221 approximant: str
222 name of approximant that was used to generate the samples
223 latex_labels: dict
224 dictionary of latex labels for each parameter
225 """
226 from pesummary.core.plots.plot import _make_corner_plot
228 if corner_parameters is None:
229 corner_parameters = conf.gw_corner_parameters
231 return _make_corner_plot(
232 samples, latex_labels, corner_parameters=corner_parameters, **kwargs
233 )
236def _make_source_corner_plot(samples, latex_labels, **kwargs):
237 """Generate the corner plots for a given approximant
239 Parameters
240 ----------
241 opts: argparse
242 argument parser object to hold all information from the command line
243 samples: nd list
244 nd list of samples for each parameter for a given approximant
245 params: list
246 list of parameters associated with each element in samples
247 approximant: str
248 name of approximant that was used to generate the samples
249 latex_labels: dict
250 dictionary of latex labels for each parameter
251 """
252 from pesummary.core.plots.plot import _make_corner_plot
254 return _make_corner_plot(
255 samples, latex_labels,
256 corner_parameters=conf.gw_source_frame_corner_parameters, **kwargs
257 )[0]
260def _make_extrinsic_corner_plot(samples, latex_labels, **kwargs):
261 """Generate the corner plots for a given approximant
263 Parameters
264 ----------
265 opts: argparse
266 argument parser object to hold all information from the command line
267 samples: nd list
268 nd list of samples for each parameter for a given approximant
269 params: list
270 list of parameters associated with each element in samples
271 approximant: str
272 name of approximant that was used to generate the samples
273 latex_labels: dict
274 dictionary of latex labels for each parameter
275 """
276 from pesummary.core.plots.plot import _make_corner_plot
278 return _make_corner_plot(
279 samples, latex_labels,
280 corner_parameters=conf.gw_extrinsic_corner_parameters, **kwargs
281 )[0]
284def _make_comparison_corner_plot(
285 samples, latex_labels, corner_parameters=None, colors=conf.corner_colors,
286 **kwargs
287):
288 """Generate a corner plot which contains multiple datasets
290 Parameters
291 ----------
292 samples: dict
293 nested dictionary containing the label as key and SamplesDict as item
294 for each dataset you wish to plot
295 latex_labels: dict
296 dictionary of latex labels for each parameter
297 corner_parameters: list, optional
298 corner parameters you wish to include in the plot
299 colors: list, optional
300 unique colors for each dataset
301 **kwargs: dict
302 all kwargs are passed to `corner.corner`
303 """
304 from pesummary.core.plots.plot import _make_comparison_corner_plot
306 if corner_parameters is None:
307 corner_parameters = conf.gw_corner_parameters
309 return _make_comparison_corner_plot(
310 samples, latex_labels, corner_parameters=corner_parameters,
311 colors=colors, **kwargs
312 )
315def __antenna_response(name, ra, dec, psi, time_gps):
316 """Calculate the antenna response function
318 Parameters
319 ----------
320 name: str
321 name of the detector you wish to calculate the antenna response
322 function for
323 ra: float
324 right ascension of the source
325 dec: float
326 declination of the source
327 psi: float
328 polarisation of the source
329 time_gps: float
330 gps time of merger
331 """
332 # Following 8 lines taken from pycbc.detector.Detector
333 from astropy.units.si import sday
334 reference_time = 1126259462.0
335 gmst_reference = Time(
336 reference_time, format='gps', scale='utc', location=(0, 0)
337 ).sidereal_time('mean').rad
338 dphase = (time_gps - reference_time) / float(sday.si.scale) * (2.0 * np.pi)
339 gmst = (gmst_reference + dphase) % (2.0 * np.pi)
340 corrected_ra = gmst - ra
341 if not LALSIMULATION:
342 raise Exception("lalsimulation could not be imported. please install "
343 "lalsuite to be able to use all features")
344 detector = lalsim.DetectorPrefixToLALDetector(str(name))
346 x0 = -np.cos(psi) * np.sin(corrected_ra) - \
347 np.sin(psi) * np.cos(corrected_ra) * np.sin(dec)
348 x1 = -np.cos(psi) * np.cos(corrected_ra) + \
349 np.sin(psi) * np.sin(corrected_ra) * np.sin(dec)
350 x2 = np.sin(psi) * np.cos(dec)
351 x = np.array([x0, x1, x2])
352 dx = detector.response.dot(x)
354 y0 = np.sin(psi) * np.sin(corrected_ra) - \
355 np.cos(psi) * np.cos(corrected_ra) * np.sin(dec)
356 y1 = np.sin(psi) * np.cos(corrected_ra) + \
357 np.cos(psi) * np.sin(corrected_ra) * np.sin(dec)
358 y2 = np.cos(psi) * np.cos(dec)
359 y = np.array([y0, y1, y2])
360 dy = detector.response.dot(y)
362 if hasattr(dx, "shape"):
363 fplus = (x * dx - y * dy).sum(axis=0)
364 fcross = (x * dy + y * dx).sum(axis=0)
365 else:
366 fplus = (x * dx - y * dy).sum()
367 fcross = (x * dy + y * dx).sum()
369 return fplus, fcross
372@no_latex_plot
373def _waveform_plot(detectors, maxL_params, **kwargs):
374 """Plot the maximum likelihood waveform for a given approximant.
376 Parameters
377 ----------
378 detectors: list
379 list of detectors that you want to generate waveforms for
380 maxL_params: dict
381 dictionary of maximum likelihood parameter values
382 kwargs: dict
383 dictionary of optional keyword arguments
384 """
385 from gwpy.plot.colors import GW_OBSERVATORY_COLORS
386 if math.isnan(maxL_params["mass_1"]):
387 return
388 logger.debug("Generating the maximum likelihood waveform plot")
389 if not LALSIMULATION:
390 raise Exception("lalsimulation could not be imported. please install "
391 "lalsuite to be able to use all features")
392 delta_frequency = kwargs.get("delta_f", 1. / 256)
393 minimum_frequency = kwargs.get("f_min", 5.)
394 maximum_frequency = kwargs.get("f_max", 1000.)
395 frequency_array = np.arange(minimum_frequency, maximum_frequency,
396 delta_frequency)
398 approx = lalsim.GetApproximantFromString(maxL_params["approximant"])
399 mass_1 = maxL_params["mass_1"] * MSUN_SI
400 mass_2 = maxL_params["mass_2"] * MSUN_SI
401 luminosity_distance = maxL_params["luminosity_distance"] * PC_SI * 10**6
402 if "phi_jl" in maxL_params.keys():
403 iota, S1x, S1y, S1z, S2x, S2y, S2z = \
404 lalsim.SimInspiralTransformPrecessingNewInitialConditions(
405 maxL_params["theta_jn"], maxL_params["phi_jl"], maxL_params["tilt_1"],
406 maxL_params["tilt_2"], maxL_params["phi_12"], maxL_params["a_1"],
407 maxL_params["a_2"], mass_1, mass_2, kwargs.get("f_ref", 10.),
408 maxL_params["phase"])
409 else:
410 iota, S1x, S1y, S1z, S2x, S2y, S2z = maxL_params["iota"], 0., 0., 0., \
411 0., 0., 0.
412 phase = maxL_params["phase"] if "phase" in maxL_params.keys() else 0.0
413 for func in [lalsim.SimInspiralChooseFDWaveform, lalsim.SimInspiralFD]:
414 try:
415 h_plus, h_cross = func(
416 mass_1, mass_2, S1x, S1y, S1z, S2x, S2y, S2z, luminosity_distance, iota,
417 phase, 0.0, 0.0, 0.0, delta_frequency, minimum_frequency,
418 maximum_frequency, kwargs.get("f_ref", 10.), None, approx)
419 except Exception:
420 continue
421 break
422 h_plus = h_plus.data.data
423 h_cross = h_cross.data.data
424 h_plus = h_plus[:len(frequency_array)]
425 h_cross = h_cross[:len(frequency_array)]
426 fig, ax = figure(gca=True)
427 colors = [GW_OBSERVATORY_COLORS[i] for i in detectors]
428 for num, i in enumerate(detectors):
429 ar = __antenna_response(i, maxL_params["ra"], maxL_params["dec"],
430 maxL_params["psi"], maxL_params["geocent_time"])
431 ax.plot(frequency_array, abs(h_plus * ar[0] + h_cross * ar[1]),
432 color=colors[num], linewidth=1.0, label=i)
433 ax.set_xscale("log")
434 ax.set_yscale("log")
435 ax.set_xlabel(r"Frequency $[Hz]$")
436 ax.set_ylabel(r"Strain")
437 ax.grid(visible=True)
438 ax.legend(loc="best")
439 fig.tight_layout()
440 return fig
443@no_latex_plot
444def _waveform_comparison_plot(maxL_params_list, colors, labels,
445 **kwargs):
446 """Generate a plot which compares the maximum likelihood waveforms for
447 each approximant.
449 Parameters
450 ----------
451 maxL_params_list: list
452 list of dictionaries containing the maximum likelihood parameter
453 values for each approximant
454 colors: list
455 list of colors to be used to differentiate the different approximants
456 approximant_labels: list, optional
457 label to prepend the approximant in the legend
458 kwargs: dict
459 dictionary of optional keyword arguments
460 """
461 logger.debug("Generating the maximum likelihood waveform comparison plot "
462 "for H1")
463 if not LALSIMULATION:
464 raise Exception("LALSimulation could not be imported. Please install "
465 "LALSuite to be able to use all features")
466 delta_frequency = kwargs.get("delta_f", 1. / 256)
467 minimum_frequency = kwargs.get("f_min", 5.)
468 maximum_frequency = kwargs.get("f_max", 1000.)
469 frequency_array = np.arange(minimum_frequency, maximum_frequency,
470 delta_frequency)
472 fig, ax = figure(gca=True)
473 for num, i in enumerate(maxL_params_list):
474 if math.isnan(i["mass_1"]):
475 continue
476 approx = lalsim.GetApproximantFromString(i["approximant"])
477 mass_1 = i["mass_1"] * MSUN_SI
478 mass_2 = i["mass_2"] * MSUN_SI
479 luminosity_distance = i["luminosity_distance"] * PC_SI * 10**6
480 if "phi_jl" in i.keys():
481 iota, S1x, S1y, S1z, S2x, S2y, S2z = \
482 lalsim.SimInspiralTransformPrecessingNewInitialConditions(
483 i["theta_jn"], i["phi_jl"], i["tilt_1"],
484 i["tilt_2"], i["phi_12"], i["a_1"],
485 i["a_2"], mass_1, mass_2, kwargs.get("f_ref", 10.),
486 i["phase"])
487 else:
488 iota, S1x, S1y, S1z, S2x, S2y, S2z = i["iota"], 0., 0., 0., \
489 0., 0., 0.
490 phase = i["phase"] if "phase" in i.keys() else 0.0
491 for func in [lalsim.SimInspiralChooseFDWaveform, lalsim.SimInspiralFD]:
492 try:
493 h_plus, h_cross = func(
494 mass_1, mass_2, S1x, S1y, S1z, S2x, S2y, S2z, luminosity_distance,
495 iota, phase, 0.0, 0.0, 0.0, delta_frequency, minimum_frequency,
496 maximum_frequency, kwargs.get("f_ref", 10.), None, approx)
497 except Exception:
498 continue
499 break
500 h_plus = h_plus.data.data
501 h_cross = h_cross.data.data
502 h_plus = h_plus[:len(frequency_array)]
503 h_cross = h_cross[:len(frequency_array)]
504 ar = __antenna_response("H1", i["ra"], i["dec"], i["psi"],
505 i["geocent_time"])
506 ax.plot(frequency_array, abs(h_plus * ar[0] + h_cross * ar[1]),
507 color=colors[num], label=labels[num], linewidth=2.0)
508 ax.set_xscale("log")
509 ax.set_yscale("log")
510 ax.grid(visible=True)
511 ax.legend(loc="best")
512 ax.set_xlabel(r"Frequency $[Hz]$")
513 ax.set_ylabel(r"Strain")
514 fig.tight_layout()
515 return fig
518def _ligo_skymap_plot(ra, dec, dist=None, savedir="./", nprocess=1,
519 downsampled=False, label="pesummary", time=None,
520 distance_map=True, multi_resolution=True,
521 injection=None, **kwargs):
522 """Plot the sky location of the source for a given approximant using the
523 ligo.skymap package
525 Parameters
526 ----------
527 ra: list
528 list of samples for right ascension
529 dec: list
530 list of samples for declination
531 dist: list
532 list of samples for the luminosity distance
533 savedir: str
534 path to the directory where you would like to save the output files
535 nprocess: Bool
536 Boolean for whether to use multithreading or not
537 downsampled: Bool
538 Boolean for whether the samples have been downsampled or not
539 distance_map: Bool
540 Boolean for whether or not to produce a distance map
541 multi_resolution: Bool
542 Boolean for whether or not to generate a multiresolution HEALPix map
543 injection: list, optional
544 List containing RA and DEC of the injection. Both must be in radians
545 kwargs: dict
546 optional keyword arguments
547 """
548 from ligo.skymap.bayestar import rasterize
549 from ligo.skymap import io
550 from ligo.skymap.kde import Clustered2DSkyKDE, Clustered2Plus1DSkyKDE
552 if dist is not None and distance_map:
553 pts = np.column_stack((ra, dec, dist))
554 cls = Clustered2Plus1DSkyKDE
555 else:
556 pts = np.column_stack((ra, dec))
557 cls = Clustered2DSkyKDE
558 skypost = cls(pts, trials=5, jobs=nprocess)
559 hpmap = skypost.as_healpix()
560 if not multi_resolution:
561 hpmap = rasterize(hpmap)
562 hpmap.meta['creator'] = "pesummary"
563 hpmap.meta['origin'] = 'LIGO/Virgo'
564 hpmap.meta['gps_creation_time'] = Time.now().gps
565 if dist is not None:
566 hpmap.meta["distmean"] = float(np.mean(dist))
567 hpmap.meta["diststd"] = float(np.std(dist))
568 if time is not None:
569 if isinstance(time, (float, int)):
570 _time = time
571 else:
572 _time = np.mean(time)
573 hpmap.meta["gps_time"] = _time
575 io.write_sky_map(
576 os.path.join(savedir, "%s_skymap.fits" % (label)), hpmap, nest=True
577 )
578 skymap, metadata = io.fits.read_sky_map(
579 os.path.join(savedir, "%s_skymap.fits" % (label)), nest=None
580 )
581 return _ligo_skymap_plot_from_array(
582 skymap, nsamples=len(ra), downsampled=downsampled, injection=injection
583 )[0]
586def _ligo_skymap_plot_from_array(
587 skymap, nsamples=None, downsampled=False, contour=[50, 90],
588 annotate=True, ax=None, colors="k", injection=None
589):
590 """Generate a skymap with `ligo.skymap` based on an array of probabilities
592 Parameters
593 ----------
594 skymap: np.array
595 array of probabilities
596 nsamples: int, optional
597 number of samples used
598 downsampled: Bool, optional
599 If True, add a header to the skymap saying that this plot is downsampled
600 contour: list, optional
601 list of contours to be plotted on the skymap. Default 50, 90
602 annotate: Bool, optional
603 If True, annotate the figure by adding the 90% and 50% sky areas
604 by default
605 ax: matplotlib.axes._subplots.AxesSubplot, optional
606 Existing axis to add the plot to
607 colors: str/list
608 colors to use for the contours
609 injection: list, optional
610 List containing RA and DEC of the injection. Both must be in radians
611 """
612 import healpy as hp
613 from ligo.skymap import plot
615 if ax is None:
616 fig = figure(gca=False)
617 ax = fig.add_subplot(111, projection='astro hours mollweide')
618 ax.grid(visible=True)
620 nside = hp.npix2nside(len(skymap))
621 deg2perpix = hp.nside2pixarea(nside, degrees=True)
622 probperdeg2 = skymap / deg2perpix
624 if downsampled:
625 ax.set_title("Downsampled to %s" % (nsamples), fontdict={'fontsize': 11})
627 vmax = probperdeg2.max()
628 ax.imshow_hpx((probperdeg2, 'ICRS'), nested=True, vmin=0.,
629 vmax=vmax, cmap="cylon")
630 cls, cs = _ligo_skymap_contours(ax, skymap, contour=contour, colors=colors)
631 if annotate:
632 text = []
633 pp = np.round(contour).astype(int)
634 ii = np.round(
635 np.searchsorted(np.sort(cls), contour) * deg2perpix).astype(int)
636 for i, p in zip(ii, pp):
637 text.append(u'{:d}% area: {:d} deg²'.format(p, i, grouping=True))
638 ax.text(1, 1.05, '\n'.join(text), transform=ax.transAxes, ha='right',
639 fontsize=10)
640 plot.outline_text(ax)
641 if injection is not None and len(injection) == 2:
642 from astropy.coordinates import SkyCoord
643 from astropy import units as u
645 _inj = SkyCoord(*injection, unit=u.rad)
646 ax.scatter(
647 _inj.ra.value, _inj.dec.value, marker="*", color="orange",
648 edgecolors='k', linewidth=1.75, s=100, zorder=100,
649 transform=ax.get_transform('world')
650 )
651 return ExistingFigure(fig), ax
654def _ligo_skymap_comparion_plot_from_array(
655 skymaps, colors, labels, contour=[50, 90], show_probability_map=False,
656 injection=None
657):
658 """Generate a skymap with `ligo.skymap` based which compares arrays of
659 probabilities
661 Parameters
662 ----------
663 skymaps: list
664 list of skymap probabilities
665 colors: list
666 list of colors to use for each skymap
667 labels: list
668 list of labels associated with each skymap
669 contour: list, optional
670 contours you wish to display on the comparison plot
671 show_probability_map: int, optional
672 the index of the skymap you wish to show the probability
673 map for. Default False
674 injection: list, optional
675 List containing RA and DEC of the injection. Both must be in radians
676 """
677 ncols = number_of_columns_for_legend(labels)
678 fig = figure(gca=False)
679 ax = fig.add_subplot(111, projection='astro hours mollweide')
680 ax.grid(visible=True)
681 for num, skymap in enumerate(skymaps):
682 if isinstance(show_probability_map, int) and show_probability_map == num:
683 _, ax = _ligo_skymap_plot_from_array(
684 skymap, nsamples=None, downsampled=False, contour=contour,
685 annotate=False, ax=ax, colors=colors[num], injection=injection,
686 )
687 cls, cs = _ligo_skymap_contours(
688 ax, skymap, contour=contour, colors=colors[num]
689 )
690 cs.collections[0].set_label(labels[num])
691 ax.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, borderaxespad=0.,
692 mode="expand", ncol=ncols)
693 return fig
696def _ligo_skymap_contours(ax, skymap, contour=[50, 90], colors='k'):
697 """Plot contours on a ligo.skymap skymap
699 Parameters
700 ----------
701 ax: matplotlib.axes._subplots.AxesSubplot, optional
702 Existing axis to add the plot to
703 skymap: np.array
704 array of probabilities
705 contour: list
706 list contours you wish to plot
707 colors: str/list
708 colors to use for the contours
709 """
710 from ligo.skymap import postprocess
712 cls = 100 * postprocess.find_greedy_credible_levels(skymap)
713 cs = ax.contour_hpx((cls, 'ICRS'), nested=True, colors=colors,
714 linewidths=0.5, levels=contour)
715 ax.clabel(cs, fmt=r'%g\%%', fontsize=6, inline=True)
716 return cls, cs
719def _default_skymap_plot(ra, dec, weights=None, injection=None, **kwargs):
720 """Plot the default sky location of the source for a given approximant
722 Parameters
723 ----------
724 ra: list
725 list of samples for right ascension
726 dec: list
727 list of samples for declination
728 injection: list, optional
729 list containing the injected value of ra and dec
730 kwargs: dict
731 optional keyword arguments
732 """
733 from .cmap import register_cylon, unregister_cylon
734 # register the cylon cmap
735 register_cylon()
736 ra = [-i + np.pi for i in ra]
737 logger.debug("Generating the sky map plot")
738 fig, ax = figure(gca=True)
739 ax = fig.add_subplot(
740 111, projection="mollweide",
741 facecolor=(1.0, 0.939165516411, 0.880255669068)
742 )
743 ax.cla()
744 ax.set_title("Preliminary", fontdict={'fontsize': 11})
745 ax.grid(visible=True)
746 ax.set_xticklabels([
747 r"$2^{h}$", r"$4^{h}$", r"$6^{h}$", r"$8^{h}$", r"$10^{h}$",
748 r"$12^{h}$", r"$14^{h}$", r"$16^{h}$", r"$18^{h}$", r"$20^{h}$",
749 r"$22^{h}$"])
750 levels = [0.9, 0.5]
752 if weights is None:
753 H, X, Y = np.histogram2d(ra, dec, bins=50)
754 else:
755 H, X, Y = np.histogram2d(ra, dec, bins=50, weights=weights)
756 H = gaussian_filter(H, kwargs.get("smooth", 0.9))
757 Hflat = H.flatten()
758 indicies = np.argsort(Hflat)[::-1]
759 Hflat = Hflat[indicies]
761 CF = np.cumsum(Hflat)
762 CF /= CF[-1]
764 V = np.empty(len(levels))
765 for num, i in enumerate(levels):
766 try:
767 V[num] = Hflat[CF <= i][-1]
768 except Exception:
769 V[num] = Hflat[0]
770 V.sort()
771 m = np.diff(V) == 0
772 while np.any(m):
773 V[np.where(m)[0][0]] *= 1.0 - 1e-4
774 m = np.diff(V) == 0
775 V.sort()
776 X1, Y1 = 0.5 * (X[1:] + X[:-1]), 0.5 * (Y[1:] + Y[:-1])
778 H2 = H.min() + np.zeros((H.shape[0] + 4, H.shape[1] + 4))
779 H2[2:-2, 2:-2] = H
780 H2[2:-2, 1] = H[:, 0]
781 H2[2:-2, -2] = H[:, -1]
782 H2[1, 2:-2] = H[0]
783 H2[-2, 2:-2] = H[-1]
784 H2[1, 1] = H[0, 0]
785 H2[1, -2] = H[0, -1]
786 H2[-2, 1] = H[-1, 0]
787 H2[-2, -2] = H[-1, -1]
788 X2 = np.concatenate([X1[0] + np.array([-2, -1]) * np.diff(X1[:2]), X1,
789 X1[-1] + np.array([1, 2]) * np.diff(X1[-2:]), ])
790 Y2 = np.concatenate([Y1[0] + np.array([-2, -1]) * np.diff(Y1[:2]), Y1,
791 Y1[-1] + np.array([1, 2]) * np.diff(Y1[-2:]), ])
793 ax.pcolormesh(X2, Y2, H2.T, vmin=0., vmax=H2.T.max(), cmap="cylon")
794 cs = ax.contour(X2, Y2, H2.T, V, colors="k", linewidths=0.5)
795 if injection is not None:
796 ax.scatter(
797 -injection[0] + np.pi, injection[1], marker="*",
798 color=conf.injection_color, edgecolors='k', linewidth=1.75, s=100
799 )
800 fmt = {l: s for l, s in zip(cs.levels, [r"$90\%$", r"$50\%$"])}
801 ax.clabel(cs, fmt=fmt, fontsize=8, inline=True)
802 text = []
803 for i, j in zip(cs.collections, [90, 50]):
804 area = 0.
805 for k in i.get_paths():
806 x = k.vertices[:, 0]
807 y = k.vertices[:, 1]
808 area += 0.5 * np.sum(y[:-1] * np.diff(x) - x[:-1] * np.diff(y))
809 area = int(np.abs(area) * (180 / np.pi) * (180 / np.pi))
810 text.append(u'{:d}% area: {:d} deg²'.format(
811 int(j), area, grouping=True))
812 ax.text(1, 1.05, '\n'.join(text[::-1]), transform=ax.transAxes, ha='right',
813 fontsize=10)
814 xticks = np.arange(-np.pi, np.pi + np.pi / 6, np.pi / 4)
815 ax.set_xticks(xticks)
816 ax.set_yticks([-np.pi / 3, -np.pi / 6, 0, np.pi / 6, np.pi / 3])
817 labels = [r"$%s^{h}$" % (int(np.round((i + np.pi) * 3.82, 1))) for i in xticks]
818 ax.set_xticklabels(labels[::-1], fontsize=10)
819 ax.set_yticklabels([r"$-60^{\circ}$", r"$-30^{\circ}$", r"$0^{\circ}$",
820 r"$30^{\circ}$", r"$60^{\circ}$"], fontsize=10)
821 ax.grid(visible=True)
822 # unregister the cylon cmap
823 unregister_cylon()
824 return fig
827def _sky_map_comparison_plot(ra_list, dec_list, labels, colors, **kwargs):
828 """Generate a plot that compares the sky location for multiple approximants
830 Parameters
831 ----------
832 ra_list: 2d list
833 list of samples for right ascension for each approximant
834 dec_list: 2d list
835 list of samples for declination for each approximant
836 approximants: list
837 list of approximants used to generate the samples
838 colors: list
839 list of colors to be used to differentiate the different approximants
840 approximant_labels: list, optional
841 label to prepend the approximant in the legend
842 kwargs: dict
843 optional keyword arguments
844 """
845 ra_list = [[-i + np.pi for i in j] for j in ra_list]
846 logger.debug("Generating the sky map comparison plot")
847 fig = figure(gca=False)
848 ax = fig.add_subplot(
849 111, projection="mollweide",
850 facecolor=(1.0, 0.939165516411, 0.880255669068)
851 )
852 ax.cla()
853 ax.grid(visible=True)
854 ax.set_xticklabels([
855 r"$2^{h}$", r"$4^{h}$", r"$6^{h}$", r"$8^{h}$", r"$10^{h}$",
856 r"$12^{h}$", r"$14^{h}$", r"$16^{h}$", r"$18^{h}$", r"$20^{h}$",
857 r"$22^{h}$"])
858 levels = [0.9, 0.5]
859 for num, i in enumerate(ra_list):
860 H, X, Y = np.histogram2d(i, dec_list[num], bins=50)
861 H = gaussian_filter(H, kwargs.get("smooth", 0.9))
862 Hflat = H.flatten()
863 indicies = np.argsort(Hflat)[::-1]
864 Hflat = Hflat[indicies]
866 CF = np.cumsum(Hflat)
867 CF /= CF[-1]
869 V = np.empty(len(levels))
870 for num2, j in enumerate(levels):
871 try:
872 V[num2] = Hflat[CF <= j][-1]
873 except Exception:
874 V[num2] = Hflat[0]
875 V.sort()
876 m = np.diff(V) == 0
877 while np.any(m):
878 V[np.where(m)[0][0]] *= 1.0 - 1e-4
879 m = np.diff(V) == 0
880 V.sort()
881 X1, Y1 = 0.5 * (X[1:] + X[:-1]), 0.5 * (Y[1:] + Y[:-1])
883 H2 = H.min() + np.zeros((H.shape[0] + 4, H.shape[1] + 4))
884 H2[2:-2, 2:-2] = H
885 H2[2:-2, 1] = H[:, 0]
886 H2[2:-2, -2] = H[:, -1]
887 H2[1, 2:-2] = H[0]
888 H2[-2, 2:-2] = H[-1]
889 H2[1, 1] = H[0, 0]
890 H2[1, -2] = H[0, -1]
891 H2[-2, 1] = H[-1, 0]
892 H2[-2, -2] = H[-1, -1]
893 X2 = np.concatenate([X1[0] + np.array([-2, -1]) * np.diff(X1[:2]), X1,
894 X1[-1] + np.array([1, 2]) * np.diff(X1[-2:]), ])
895 Y2 = np.concatenate([Y1[0] + np.array([-2, -1]) * np.diff(Y1[:2]), Y1,
896 Y1[-1] + np.array([1, 2]) * np.diff(Y1[-2:]), ])
897 CS = ax.contour(X2, Y2, H2.T, V, colors=colors[num], linewidths=2.0)
898 CS.collections[0].set_label(labels[num])
899 ncols = number_of_columns_for_legend(labels)
900 ax.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, borderaxespad=0.,
901 mode="expand", ncol=ncols)
902 xticks = np.arange(-np.pi, np.pi + np.pi / 6, np.pi / 4)
903 ax.set_xticks(xticks)
904 ax.set_yticks([-np.pi / 3, -np.pi / 6, 0, np.pi / 6, np.pi / 3])
905 labels = [r"$%s^{h}$" % (int(np.round((i + np.pi) * 3.82, 1))) for i in xticks]
906 ax.set_xticklabels(labels[::-1], fontsize=10)
907 ax.set_yticklabels([r"$-60^\degree$", r"$-30^\degree$", r"$0^\degree$",
908 r"$30^\degree$", r"$60^\degree$"], fontsize=10)
909 ax.grid(visible=True)
910 return fig
913def __get_cutoff_indices(flow, fhigh, df, N):
914 """
915 Gets the indices of a frequency series at which to stop an overlap
916 calculation.
918 Parameters
919 ----------
920 flow: float
921 The frequency (in Hz) of the lower index.
922 fhigh: float
923 The frequency (in Hz) of the upper index.
924 df: float
925 The frequency step (in Hz) of the frequency series.
926 N: int
927 The number of points in the **time** series. Can be odd
928 or even.
930 Returns
931 -------
932 kmin: int
933 kmax: int
934 """
935 if flow:
936 kmin = int(flow / df)
937 else:
938 kmin = 1
939 if fhigh:
940 kmax = int(fhigh / df)
941 else:
942 kmax = int((N + 1) / 2.)
943 return kmin, kmax
946@no_latex_plot
947def _sky_sensitivity(network, resolution, maxL_params, **kwargs):
948 """Generate the sky sensitivity for a given network
950 Parameters
951 ----------
952 network: list
953 list of detectors you want included in your sky sensitivity plot
954 resolution: float
955 resolution of the skymap
956 maxL_params: dict
957 dictionary of waveform parameters for the maximum likelihood waveform
958 """
959 logger.debug("Generating the sky sensitivity for %s" % (network))
960 if not LALSIMULATION:
961 raise Exception("LALSimulation could not be imported. Please install "
962 "LALSuite to be able to use all features")
963 delta_frequency = kwargs.get("delta_f", 1. / 256)
964 minimum_frequency = kwargs.get("f_min", 20.)
965 maximum_frequency = kwargs.get("f_max", 1000.)
966 frequency_array = np.arange(minimum_frequency, maximum_frequency,
967 delta_frequency)
969 approx = lalsim.GetApproximantFromString(maxL_params["approximant"])
970 mass_1 = maxL_params["mass_1"] * MSUN_SI
971 mass_2 = maxL_params["mass_2"] * MSUN_SI
972 luminosity_distance = maxL_params["luminosity_distance"] * PC_SI * 10**6
973 iota, S1x, S1y, S1z, S2x, S2y, S2z = \
974 lalsim.SimInspiralTransformPrecessingNewInitialConditions(
975 maxL_params["iota"], maxL_params["phi_jl"], maxL_params["tilt_1"],
976 maxL_params["tilt_2"], maxL_params["phi_12"], maxL_params["a_1"],
977 maxL_params["a_2"], mass_1, mass_2, kwargs.get("f_ref", 10.),
978 maxL_params["phase"])
979 h_plus, h_cross = lalsim.SimInspiralChooseFDWaveform(
980 mass_1, mass_2, S1x, S1y, S1z, S2x, S2y, S2z, luminosity_distance, iota,
981 maxL_params["phase"], 0.0, 0.0, 0.0, delta_frequency, minimum_frequency,
982 maximum_frequency, kwargs.get("f_ref", 10.), None, approx)
983 h_plus = h_plus.data.data
984 h_cross = h_cross.data.data
985 h_plus = h_plus[:len(frequency_array)]
986 h_cross = h_cross[:len(frequency_array)]
987 psd = {}
988 psd["H1"] = psd["L1"] = np.array([
989 lalsim.SimNoisePSDaLIGOZeroDetHighPower(i) for i in frequency_array])
990 psd["V1"] = np.array([lalsim.SimNoisePSDVirgo(i) for i in frequency_array])
991 kmin, kmax = __get_cutoff_indices(minimum_frequency, maximum_frequency,
992 delta_frequency, (len(h_plus) - 1) * 2)
993 ra = np.arange(-np.pi, np.pi, resolution)
994 dec = np.arange(-np.pi, np.pi, resolution)
995 X, Y = np.meshgrid(ra, dec)
996 N = np.zeros([len(dec), len(ra)])
998 indices = np.ndindex(len(ra), len(dec))
999 for ind in indices:
1000 ar = {}
1001 SNR = {}
1002 for i in network:
1003 ard = __antenna_response(i, ra[ind[0]], dec[ind[1]],
1004 maxL_params["psi"], maxL_params["geocent_time"])
1005 ar[i] = [ard[0], ard[1]]
1006 strain = np.array(h_plus * ar[i][0] + h_cross * ar[i][1])
1007 integrand = np.conj(strain[kmin:kmax]) * strain[kmin:kmax] / psd[i][kmin:kmax]
1008 integrand = integrand[:-1]
1009 SNR[i] = np.sqrt(4 * delta_frequency * np.sum(integrand).real)
1010 ar[i][0] *= SNR[i]
1011 ar[i][1] *= SNR[i]
1012 numerator = 0.0
1013 denominator = 0.0
1014 for i in network:
1015 numerator += sum(i**2 for i in ar[i])
1016 denominator += SNR[i]**2
1017 N[ind[1]][ind[0]] = (((numerator / denominator)**0.5))
1018 fig = figure(gca=False)
1019 ax = fig.add_subplot(111, projection="hammer")
1020 ax.cla()
1021 ax.grid(visible=True)
1022 ax.pcolormesh(X, Y, N)
1023 ax.set_xticklabels([
1024 r"$22^{h}$", r"$20^{h}$", r"$18^{h}$", r"$16^{h}$", r"$14^{h}$",
1025 r"$12^{h}$", r"$10^{h}$", r"$8^{h}$", r"$6^{h}$", r"$4^{h}$",
1026 r"$2^{h}$"])
1027 return fig
1030@no_latex_plot
1031def _time_domain_waveform(detectors, maxL_params, **kwargs):
1032 """
1033 Plot the maximum likelihood waveform for a given approximant
1034 in the time domain.
1036 Parameters
1037 ----------
1038 detectors: list
1039 list of detectors that you want to generate waveforms for
1040 maxL_params: dict
1041 dictionary of maximum likelihood parameter values
1042 kwargs: dict
1043 dictionary of optional keyword arguments
1044 """
1045 from gwpy.timeseries import TimeSeries
1046 from gwpy.plot.colors import GW_OBSERVATORY_COLORS
1047 if math.isnan(maxL_params["mass_1"]):
1048 return
1049 logger.debug("Generating the maximum likelihood waveform time domain plot")
1050 if not LALSIMULATION:
1051 raise Exception("lalsimulation could not be imported. please install "
1052 "lalsuite to be able to use all features")
1053 delta_t = 1. / 4096.
1054 minimum_frequency = kwargs.get("f_min", 5.)
1055 t_start = maxL_params['geocent_time']
1056 t_finish = maxL_params['geocent_time'] + 4.
1057 time_array = np.arange(t_start, t_finish, delta_t)
1059 approx = lalsim.GetApproximantFromString(maxL_params["approximant"])
1060 mass_1 = maxL_params["mass_1"] * MSUN_SI
1061 mass_2 = maxL_params["mass_2"] * MSUN_SI
1062 luminosity_distance = maxL_params["luminosity_distance"] * PC_SI * 10**6
1063 if "phi_jl" in maxL_params.keys():
1064 iota, S1x, S1y, S1z, S2x, S2y, S2z = \
1065 lalsim.SimInspiralTransformPrecessingNewInitialConditions(
1066 maxL_params["theta_jn"], maxL_params["phi_jl"], maxL_params["tilt_1"],
1067 maxL_params["tilt_2"], maxL_params["phi_12"], maxL_params["a_1"],
1068 maxL_params["a_2"], mass_1, mass_2, kwargs.get("f_ref", 10.),
1069 maxL_params["phase"])
1070 else:
1071 iota, S1x, S1y, S1z, S2x, S2y, S2z = maxL_params["iota"], 0., 0., 0., \
1072 0., 0., 0.
1073 phase = maxL_params["phase"] if "phase" in maxL_params.keys() else 0.0
1074 chirptime = lalsim.SimIMRPhenomXASDuration(
1075 mass_1, mass_2, S1z, S2z, minimum_frequency
1076 )
1077 duration = np.max([2**np.ceil(np.log2(chirptime)), 1.0])
1078 h_plus, h_cross = lalsim.SimInspiralChooseTDWaveform(
1079 mass_1, mass_2, S1x, S1y, S1z, S2x, S2y, S2z, luminosity_distance, iota,
1080 phase, 0.0, 0.0, 0.0, delta_t, minimum_frequency,
1081 kwargs.get("f_ref", 10.), None, approx)
1083 fig, ax = figure(gca=True)
1084 colors = [GW_OBSERVATORY_COLORS[i] for i in detectors]
1085 for num, i in enumerate(detectors):
1086 ar = __antenna_response(i, maxL_params["ra"], maxL_params["dec"],
1087 maxL_params["psi"], maxL_params["geocent_time"])
1088 h_t = h_plus.data.data * ar[0] + h_cross.data.data * ar[1]
1089 h_t = TimeSeries(h_t[:], dt=h_plus.deltaT, t0=h_plus.epoch)
1090 h_t.times = [float(np.array(i)) + t_start for i in h_t.times]
1091 ax.plot(h_t.times, h_t,
1092 color=colors[num], linewidth=1.0, label=i)
1093 ax.set_xlim([t_start - 0.75 * duration, t_start + duration / 4])
1094 ax.set_xlabel(r"Time $[s]$")
1095 ax.set_ylabel(r"Strain")
1096 ax.grid(visible=True)
1097 ax.legend(loc="best")
1098 fig.tight_layout()
1099 return fig
1102@no_latex_plot
1103def _time_domain_waveform_comparison_plot(maxL_params_list, colors, labels,
1104 **kwargs):
1105 """Generate a plot which compares the maximum likelihood waveforms for
1106 each approximant.
1108 Parameters
1109 ----------
1110 maxL_params_list: list
1111 list of dictionaries containing the maximum likelihood parameter
1112 values for each approximant
1113 colors: list
1114 list of colors to be used to differentiate the different approximants
1115 approximant_labels: list, optional
1116 label to prepend the approximant in the legend
1117 kwargs: dict
1118 dictionary of optional keyword arguments
1119 """
1120 from gwpy.timeseries import TimeSeries
1121 logger.debug("Generating the maximum likelihood time domain waveform "
1122 "comparison plot for H1")
1123 if not LALSIMULATION:
1124 raise Exception("LALSimulation could not be imported. Please install "
1125 "LALSuite to be able to use all features")
1126 delta_t = 1. / 4096.
1127 minimum_frequency = kwargs.get("f_min", 5.)
1129 fig, ax = figure(gca=True)
1130 for num, i in enumerate(maxL_params_list):
1131 if math.isnan(i["mass_1"]):
1132 continue
1133 t_start = i['geocent_time']
1134 t_finish = i['geocent_time'] + 4.
1135 time_array = np.arange(t_start, t_finish, delta_t)
1137 approx = lalsim.GetApproximantFromString(i["approximant"])
1138 mass_1 = i["mass_1"] * MSUN_SI
1139 mass_2 = i["mass_2"] * MSUN_SI
1140 luminosity_distance = i["luminosity_distance"] * PC_SI * 10**6
1141 if "phi_jl" in i.keys():
1142 iota, S1x, S1y, S1z, S2x, S2y, S2z = \
1143 lalsim.SimInspiralTransformPrecessingNewInitialConditions(
1144 i["theta_jn"], i["phi_jl"], i["tilt_1"],
1145 i["tilt_2"], i["phi_12"], i["a_1"],
1146 i["a_2"], mass_1, mass_2, kwargs.get("f_ref", 10.),
1147 i["phase"])
1148 else:
1149 iota, S1x, S1y, S1z, S2x, S2y, S2z = i["iota"], 0., 0., 0., \
1150 0., 0., 0.
1151 phase = i["phase"] if "phase" in i.keys() else 0.0
1152 h_plus, h_cross = lalsim.SimInspiralChooseTDWaveform(
1153 mass_1, mass_2, S1x, S1y, S1z, S2x, S2y, S2z, luminosity_distance,
1154 iota, phase, 0.0, 0.0, 0.0, delta_t, minimum_frequency,
1155 kwargs.get("f_ref", 10.), None, approx)
1157 ar = __antenna_response("H1", i["ra"], i["dec"], i["psi"],
1158 i["geocent_time"])
1159 h_t = h_plus.data.data * ar[0] + h_cross.data.data * ar[1]
1160 h_t = TimeSeries(h_t[:], dt=h_plus.deltaT, t0=h_plus.epoch)
1161 h_t.times = [float(np.array(i)) + t_start for i in h_t.times]
1163 ax.plot(h_t.times, h_t,
1164 color=colors[num], label=labels[num], linewidth=2.0)
1165 ax.set_xlabel(r"Time $[s]$")
1166 ax.set_ylabel(r"Strain")
1167 ax.set_xlim([t_start - 3, t_start + 0.5])
1168 ax.grid(visible=True)
1169 ax.legend(loc="best")
1170 fig.tight_layout()
1171 return fig
1174def _psd_plot(frequencies, strains, colors=None, labels=None, fmin=None, fmax=None):
1175 """Superimpose all PSD plots onto a single figure.
1177 Parameters
1178 ----------
1179 frequencies: nd list
1180 list of all frequencies used for each psd file
1181 strains: nd list
1182 list of all strains used for each psd file
1183 colors: optional, list
1184 list of colors to be used to differentiate the different PSDs
1185 labels: optional, list
1186 list of lavels for each PSD
1187 fmin: optional, float
1188 starting frequency of the plot
1189 fmax: optional, float
1190 maximum frequency of the plot
1191 """
1192 from gwpy.plot.colors import GW_OBSERVATORY_COLORS
1193 fig, ax = figure(gca=True)
1194 if not colors and all(i in GW_OBSERVATORY_COLORS.keys() for i in labels):
1195 colors = [GW_OBSERVATORY_COLORS[i] for i in labels]
1196 elif not colors:
1197 colors = ['r', 'b', 'orange', 'c', 'g', 'purple']
1198 while len(colors) <= len(labels):
1199 colors += colors
1200 for num, i in enumerate(frequencies):
1201 ff = np.array(i)
1202 ss = np.array(strains[num])
1203 cond = np.ones_like(strains[num], dtype=bool)
1204 if fmin is not None:
1205 cond *= ff >= fmin
1206 if fmax is not None:
1207 cond *= ff <= fmax
1208 i = ff[cond]
1209 strains[num] = ss[cond]
1210 ax.loglog(i, strains[num], color=colors[num], label=labels[num])
1211 ax.tick_params(which="both", bottom=True, length=3, width=1)
1212 ax.set_xlabel(r"Frequency $[\mathrm{Hz}]$")
1213 ax.set_ylabel(r"Power Spectral Density [$\mathrm{strain}^{2}/\mathrm{Hz}$]")
1214 ax.legend(loc="best")
1215 fig.tight_layout()
1216 return fig
1219@no_latex_plot
1220def _calibration_envelope_plot(frequency, calibration_envelopes, ifos,
1221 colors=None, prior=[]):
1222 """Generate a plot showing the calibration envelope
1224 Parameters
1225 ----------
1226 frequency: array
1227 frequency bandwidth that you would like to use
1228 calibration_envelopes: nd list
1229 list containing the calibration envelope data for different IFOs
1230 ifos: list
1231 list of IFOs that are associated with the calibration envelopes
1232 colors: list, optional
1233 list of colors to be used to differentiate the different calibration
1234 envelopes
1235 prior: list, optional
1236 list containing the prior calibration envelope data for different IFOs
1237 """
1238 from gwpy.plot.colors import GW_OBSERVATORY_COLORS
1240 def interpolate_calibration(data):
1241 """Interpolate the calibration data using spline
1243 Parameters
1244 ----------
1245 data: np.ndarray
1246 array containing the calibration data
1247 """
1248 interp = [
1249 np.interp(frequency, data[:, 0], data[:, j], left=k, right=k)
1250 for j, k in zip(range(1, 7), [1, 0, 1, 0, 1, 0])
1251 ]
1252 amp_median = (interp[0] - 1) * 100
1253 phase_median = interp[1] * 180. / np.pi
1254 amp_lower_sigma = (interp[2] - 1) * 100
1255 phase_lower_sigma = interp[3] * 180. / np.pi
1256 amp_upper_sigma = (interp[4] - 1) * 100
1257 phase_upper_sigma = interp[5] * 180. / np.pi
1258 data_dict = {
1259 "amplitude": {
1260 "median": amp_median,
1261 "lower": amp_lower_sigma,
1262 "upper": amp_upper_sigma
1263 },
1264 "phase": {
1265 "median": phase_median,
1266 "lower": phase_lower_sigma,
1267 "upper": phase_upper_sigma
1268 }
1269 }
1270 return data_dict
1272 fig, (ax1, ax2) = subplots(2, 1, sharex=True, gca=False)
1273 if not colors and all(i in GW_OBSERVATORY_COLORS.keys() for i in ifos):
1274 colors = [GW_OBSERVATORY_COLORS[i] for i in ifos]
1275 elif not colors:
1276 colors = ['r', 'b', 'orange', 'c', 'g', 'purple']
1277 while len(colors) <= len(ifos):
1278 colors += colors
1280 for num, i in enumerate(calibration_envelopes):
1281 calibration_envelopes[num] = np.array(calibration_envelopes[num])
1283 for num, i in enumerate(calibration_envelopes):
1284 calibration_data = interpolate_calibration(i)
1285 if prior != []:
1286 prior_data = interpolate_calibration(prior[num])
1287 ax1.plot(
1288 frequency, calibration_data["amplitude"]["upper"], color=colors[num],
1289 linestyle="-", label=ifos[num]
1290 )
1291 ax1.plot(
1292 frequency, calibration_data["amplitude"]["lower"], color=colors[num],
1293 linestyle="-"
1294 )
1295 ax1.set_ylabel(r"Amplitude deviation $[\%]$", fontsize=10)
1296 ax1.legend(loc="best")
1297 ax2.plot(
1298 frequency, calibration_data["phase"]["upper"], color=colors[num],
1299 linestyle="-", label=ifos[num]
1300 )
1301 ax2.plot(
1302 frequency, calibration_data["phase"]["lower"], color=colors[num],
1303 linestyle="-"
1304 )
1305 ax2.set_ylabel(r"Phase deviation $[\degree]$", fontsize=10)
1306 if prior != []:
1307 ax1.fill_between(
1308 frequency, prior_data["amplitude"]["upper"],
1309 prior_data["amplitude"]["lower"], color=colors[num], alpha=0.2
1310 )
1311 ax2.fill_between(
1312 frequency, prior_data["phase"]["upper"],
1313 prior_data["phase"]["lower"], color=colors[num], alpha=0.2
1314 )
1316 ax1.set_xscale('log')
1317 ax2.set_xscale('log')
1318 ax2.set_xlabel(r"Frequency $[Hz]$")
1319 fig.tight_layout()
1320 return fig
1323def _strain_plot(strain, maxL_params, **kwargs):
1324 """Generate a plot showing the strain data and the maxL waveform
1326 Parameters
1327 ----------
1328 strain: gwpy.timeseries
1329 timeseries containing the strain data
1330 maxL_samples: dict
1331 dictionary of maximum likelihood parameter values
1332 """
1333 logger.debug("Generating the strain plot")
1334 from pesummary.gw.conversions import time_in_each_ifo
1335 from gwpy.timeseries import TimeSeries
1337 fig, axs = subplots(nrows=len(strain.keys()), sharex=True)
1338 time = maxL_params["geocent_time"]
1339 delta_t = 1. / 4096.
1340 minimum_frequency = kwargs.get("f_min", 5.)
1341 t_start = time - 15.0
1342 t_finish = time + 0.06
1343 time_array = np.arange(t_start, t_finish, delta_t)
1345 approx = lalsim.GetApproximantFromString(maxL_params["approximant"])
1346 mass_1 = maxL_params["mass_1"] * MSUN_SI
1347 mass_2 = maxL_params["mass_2"] * MSUN_SI
1348 luminosity_distance = maxL_params["luminosity_distance"] * PC_SI * 10**6
1349 phase = maxL_params["phase"] if "phase" in maxL_params.keys() else 0.0
1350 cartesian = [
1351 "iota", "spin_1x", "spin_1y", "spin_1z", "spin_2x", "spin_2y", "spin_2z"
1352 ]
1353 if not all(param in maxL_params.keys() for param in cartesian):
1354 if "phi_jl" in maxL_params.keys():
1355 iota, S1x, S1y, S1z, S2x, S2y, S2z = \
1356 lalsim.SimInspiralTransformPrecessingNewInitialConditions(
1357 maxL_params["theta_jn"], maxL_params["phi_jl"],
1358 maxL_params["tilt_1"], maxL_params["tilt_2"],
1359 maxL_params["phi_12"], maxL_params["a_1"],
1360 maxL_params["a_2"], mass_1, mass_2, kwargs.get("f_ref", 10.),
1361 phase
1362 )
1363 else:
1364 iota, S1x, S1y, S1z, S2x, S2y, S2z = maxL_params["iota"], 0., 0., \
1365 0., 0., 0., 0.
1366 else:
1367 iota, S1x, S1y, S1z, S2x, S2y, S2z = [
1368 maxL_params[param] for param in cartesian
1369 ]
1370 h_plus, h_cross = lalsim.SimInspiralChooseTDWaveform(
1371 mass_1, mass_2, S1x, S1y, S1z, S2x, S2y, S2z, luminosity_distance, iota,
1372 phase, 0.0, 0.0, 0.0, delta_t, minimum_frequency,
1373 kwargs.get("f_ref", 10.), None, approx)
1374 h_plus = TimeSeries(
1375 h_plus.data.data[:], dt=h_plus.deltaT, t0=h_plus.epoch
1376 )
1377 h_cross = TimeSeries(
1378 h_cross.data.data[:], dt=h_cross.deltaT, t0=h_cross.epoch
1379 )
1381 for num, key in enumerate(list(strain.keys())):
1382 ifo_time = time_in_each_ifo(key, maxL_params["ra"], maxL_params["dec"],
1383 maxL_params["geocent_time"])
1385 asd = strain[key].asd(8, 4, method="median")
1386 strain_data_frequency = strain[key].fft()
1387 asd_interp = asd.interpolate(float(np.array(strain_data_frequency.df)))
1388 asd_interp = asd_interp[:len(strain_data_frequency)]
1389 strain_data_time = (strain_data_frequency / asd_interp).ifft()
1390 strain_data_time = strain_data_time.highpass(30)
1391 strain_data_time = strain_data_time.lowpass(300)
1393 ar = __antenna_response(key, maxL_params["ra"], maxL_params["dec"],
1394 maxL_params["psi"], maxL_params["geocent_time"])
1396 h_t = ar[0] * h_plus + ar[1] * h_cross
1397 h_t_frequency = h_t.fft()
1398 asd_interp = asd.interpolate(float(np.array(h_t_frequency.df)))
1399 asd_interp = asd_interp[:len(h_t_frequency)]
1400 h_t_time = (h_t_frequency / asd_interp).ifft()
1401 h_t_time = h_t_time.highpass(30)
1402 h_t_time = h_t_time.lowpass(300)
1403 h_t_time.times = [float(np.array(i)) + ifo_time for i in h_t.times]
1405 strain_data_crop = strain_data_time.crop(ifo_time - 0.2, ifo_time + 0.06)
1406 try:
1407 h_t_time = h_t_time.crop(ifo_time - 0.2, ifo_time + 0.06)
1408 except Exception:
1409 pass
1410 max_strain = np.max(strain_data_crop).value
1412 axs[num].plot(strain_data_crop, color='grey', alpha=0.75, label="data")
1413 axs[num].plot(h_t_time, color='orange', label="template")
1414 axs[num].set_xlim([ifo_time - 0.2, ifo_time + 0.06])
1415 if not math.isnan(max_strain):
1416 axs[num].set_ylim([-max_strain * 1.5, max_strain * 1.5])
1417 axs[num].set_ylabel("Whitened %s strain" % (key), fontsize=8)
1418 axs[num].grid(False)
1419 axs[num].legend(loc="best", prop={'size': 8})
1420 axs[-1].set_xlabel("Time $[s]$", fontsize=16)
1421 fig.tight_layout()
1422 return fig
1425def _format_prob(prob):
1426 """Format the probabilities for use with _classification_plot
1427 """
1428 if prob >= 1:
1429 return '100%'
1430 elif prob <= 0:
1431 return '0%'
1432 elif prob > 0.99:
1433 return '>99%'
1434 elif prob < 0.01:
1435 return '<1%'
1436 else:
1437 return '{}%'.format(int(np.round(100 * prob)))
1440@no_latex_plot
1441def _classification_plot(classification):
1442 """Generate a bar chart showing the source classifications probabilities
1444 Parameters
1445 ----------
1446 classification: dict
1447 dictionary of source classifications
1448 """
1449 from matplotlib import rcParams
1451 original_fontsize = rcParams["font.size"]
1452 original_ylabel = rcParams["ytick.labelsize"]
1453 rcParams["font.size"] = 12
1454 rcParams["ytick.labelsize"] = 12
1455 probs, names = zip(
1456 *sorted(zip(classification.values(), classification.keys())))
1457 with matplotlib.style.context('seaborn-white'):
1458 fig, ax = figure(figsize=(2.5, 2), gca=True)
1459 ax.barh(names, probs)
1460 for i, prob in enumerate(probs):
1461 ax.annotate(_format_prob(prob), (0, i), (4, 0),
1462 textcoords='offset points', ha='left', va='center')
1463 ax.set_xlim(0, 1)
1464 ax.set_xticks([])
1465 ax.tick_params(left=False)
1466 for side in ['top', 'bottom', 'right']:
1467 ax.spines[side].set_visible(False)
1468 fig.tight_layout()
1469 rcParams["font.size"] = original_fontsize
1470 rcParams["ytick.labelsize"] = original_ylabel
1471 return fig