From 9863e1b859a9494e1841a1060d382c728be6d80a Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 19:47:11 +1100 Subject: [PATCH 01/37] Added direct access to transient classes --- redback/transient/__init__.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/redback/transient/__init__.py b/redback/transient/__init__.py index 3920ee86e..c317581c8 100644 --- a/redback/transient/__init__.py +++ b/redback/transient/__init__.py @@ -1,5 +1,12 @@ from redback.transient import afterglow, kilonova, prompt, supernova, tde, transient +from afterglow import Afterglow, LGRB, SGRB +from kilonova import Kilonova +from prompt import PromptTimeSeries +from supernova import Supernova +from tde import TDE +from transient import Transient -TRANSIENT_DICT = dict(afterglow=afterglow.Afterglow, lgrb=afterglow.LGRB, sgrb=afterglow.SGRB, - kilonova=kilonova.Kilonova, prompt=prompt.PromptTimeSeries, supernova=supernova.Supernova, - tde=tde.TDE, transient=transient.Transient) + +TRANSIENT_DICT = dict( + afterglow=Afterglow, lgrb=LGRB, sgrb=SGRB, kilonova=Kilonova, prompt=PromptTimeSeries, supernova=Supernova, + tde=TDE, transient=Transient) From cd7c4f921394b30c260eca678472723823a4364f Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 20:22:29 +1100 Subject: [PATCH 02/37] Added plotter classes and fixed some issue with reading non-existing metadata --- redback/plotting.py | 144 ++++++++++++++++++++++++++++++++- redback/transient/__init__.py | 12 +-- redback/transient/afterglow.py | 66 +++++++-------- redback/transient/transient.py | 54 +++---------- 4 files changed, 190 insertions(+), 86 deletions(-) diff --git a/redback/plotting.py b/redback/plotting.py index b92666a67..1f022f03f 100644 --- a/redback/plotting.py +++ b/redback/plotting.py @@ -1,9 +1,11 @@ import matplotlib import numpy as np import matplotlib.pyplot as plt -import redback from os.path import join +import redback +from redback.get_data.directory import afterglow_directory_structure + class MultiBandPlotter(object): @@ -111,3 +113,143 @@ def plot_multiband( plt.subplots_adjust(wspace=wspace, hspace=hspace) plt.savefig(join(self.transient.transient_dir, filename), bbox_inches="tight") return axes + + +class IntegratedFluxPlotter(object): + + def __init__(self, transient): + self.transient = transient + + def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k', **kwargs) -> matplotlib.axes.Axes: + """ + Plots the Afterglow lightcurve and returns Axes. + + Parameters + ---------- + axes : Union[matplotlib.axes.Axes, None], optional + Matplotlib axes to plot the lightcurve into. Useful for user specific modifications to the plot. + colour: str, optional + Colour of the data. + kwargs: dict + Additional keyword arguments. + + Returns + ---------- + matplotlib.axes.Axes: The axes with the plot. + """ + + if self.transient.x_err is not None: + x_err = [np.abs(self.transient.x_err[1, :]), self.transient.x_err[0, :]] + else: + x_err = None + y_err = [np.abs(self.transient.y_err[1, :]), self.transient.y_err[0, :]] + + ax = axes or plt.gca() + ax.errorbar(self.transient.x, self.transient.y, xerr=x_err, yerr=y_err, + fmt='x', c=colour, ms=1, elinewidth=2, capsize=0.) + + ax.set_xscale('log') + ax.set_yscale('log') + + ax.set_xlim(0.5 * self.transient.x[0], 2 * (self.transient.x[-1] + x_err[1][-1])) + ax.set_ylim(0.5 * min(self.transient.y), 2. * np.max(self.transient.y)) + + ax.annotate(self.transient.name, xy=(0.95, 0.9), xycoords='axes fraction', + horizontalalignment='right', size=20) + + ax.set_xlabel(r'Time since burst [s]') + ax.set_ylabel(self.transient.ylabel) + ax.tick_params(axis='x', pad=10) + + if axes is None: + plt.tight_layout() + + directory_structure = afterglow_directory_structure(grb=self.transient.name, data_mode=self.transient.data_mode) + filename = f"{self.transient.name}_lc.png" + plt.savefig(join(directory_structure.directory_path, filename)) + if axes is None: + plt.clf() + return ax + + +class LuminosityPlotter(IntegratedFluxPlotter): + pass + + +class MagnitudePlotter(object): + + def __init__(self, transient): + self.transient = transient + + def plot_data( + self, axes: matplotlib.axes.Axes = None, filters: list = None, plot_others: bool = True, + plot_save: bool = True, **plot_kwargs: dict) -> None: + """ + Plots the data. + + Parameters + ---------- + axes: matplotlib.axes.Axes, optional + Axes can be given if defaults are not satisfying + filters: list, optional + Which bands to plot. Will use default filters if None is given. + plot_others: bool, optional + Plot all bands outside filters in black without label if True. + plot_kwargs: + Additional optional plotting kwargs: + errorbar_fmt: Errorbar format ('fmt' argument in matplotlib.pyplot.errorbar) + colors: colors to be used for the bands + xlabel: Plot xlabel + ylabel: Plot ylabel + plot_label: Additional filename label appended to the default name + """ + if filters is None: + filters = self.transient.active_bands + + errorbar_fmt = plot_kwargs.get("errorbar_fmt", "x") + colors = plot_kwargs.get("colors", self.transient.get_colors(filters)) + xlabel = plot_kwargs.get("xlabel", self.transient.xlabel) + ylabel = plot_kwargs.get("ylabel", self.transient.ylabel) + plot_label = plot_kwargs.get("plot_label", "data") + + ax = axes or plt.gca() + for indices, band in zip(self.transient.list_of_band_indices, self.transient.unique_bands): + x_err = self.transient.x_err[indices] if self.transient.x_err is not None else self.transient.x_err + if band in filters: + color = colors[filters.index(band)] + label = band + elif plot_others: + color = "black" + label = None + else: + continue + ax.errorbar( + self.transient.x[indices], self.transient.y[indices], xerr=x_err, yerr=self.transient.y_err[indices], + fmt=errorbar_fmt, ms=1, color=color, elinewidth=2, capsize=0., label=label) + + ax.set_xlim(0.5 * self.transient.x[0], 1.2 * self.transient.x[-1]) + self._set_y_axis(ax) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + ax.tick_params(axis='x', pad=10) + ax.legend(ncol=2) + + if axes is None: + plt.tight_layout() + + if plot_save: + filename = f"{self.transient.name}_{self.transient.data_mode}_{plot_label}.png" + plt.savefig(join(self.transient.transient_dir, filename), bbox_inches='tight') + plt.clf() + return axes + + def _set_y_axis(self, ax): + ax.set_ylim(0.8 * min(self.transient.y), 1.2 * np.max(self.transient.y)) + ax.invert_yaxis() + + +class FluxDensityPlotter(MagnitudePlotter): + + def _set_y_axis(self, ax): + ax.set_ylim(0.5 * min(self.transient.y), 2. * np.max(self.transient.y)) + ax.set_yscale('log') diff --git a/redback/transient/__init__.py b/redback/transient/__init__.py index c317581c8..dc911a15b 100644 --- a/redback/transient/__init__.py +++ b/redback/transient/__init__.py @@ -1,10 +1,10 @@ from redback.transient import afterglow, kilonova, prompt, supernova, tde, transient -from afterglow import Afterglow, LGRB, SGRB -from kilonova import Kilonova -from prompt import PromptTimeSeries -from supernova import Supernova -from tde import TDE -from transient import Transient +from redback.transient.afterglow import Afterglow, LGRB, SGRB +from redback.transient.kilonova import Kilonova +from redback.transient.prompt import PromptTimeSeries +from redback.transient.supernova import Supernova +from redback.transient.tde import TDE +from redback.transient.transient import Transient TRANSIENT_DICT = dict( diff --git a/redback/transient/afterglow.py b/redback/transient/afterglow.py index 08ac5b3ea..8450238f8 100644 --- a/redback/transient/afterglow.py +++ b/redback/transient/afterglow.py @@ -12,7 +12,8 @@ from redback.utils import logger from redback.get_data.directory import afterglow_directory_structure -from redback.plotting import MultiBandPlotter +from redback.plotting import \ + MultiBandPlotter, IntegratedFluxPlotter, LuminosityPlotter, FluxDensityPlotter, MagnitudePlotter from redback.transient.transient import Transient dirname = os.path.dirname(__file__) @@ -235,10 +236,14 @@ def _set_data(self) -> None: """ Loads data from the meta data table and sets it to the respective attribute. """ - meta_data = pd.read_csv(self.event_table, header=0, error_bad_lines=False, delimiter='\t', dtype='str') - meta_data['BAT Photon Index (15-150 keV) (PL = simple power-law, CPL = cutoff power-law)'] = meta_data[ - 'BAT Photon Index (15-150 keV) (PL = simple power-law, CPL = cutoff power-law)'].fillna(0) - self.meta_data = meta_data + try: + meta_data = pd.read_csv(self.event_table, header=0, error_bad_lines=False, delimiter='\t', dtype='str') + meta_data['BAT Photon Index (15-150 keV) (PL = simple power-law, CPL = cutoff power-law)'] = meta_data[ + 'BAT Photon Index (15-150 keV) (PL = simple power-law, CPL = cutoff power-law)'].fillna(0) + self.meta_data = meta_data + except FileNotFoundError: + logger.warning("Meta data does not exist for this event.") + self.meta_data = None def _set_photon_index(self) -> None: """ @@ -246,11 +251,13 @@ def _set_photon_index(self) -> None: """ if not np.isnan(self.photon_index): return + if self.magnitude_data or self.flux_density_data: + self.photon_index = np.nan try: photon_index = self.meta_data.query('GRB == @self._stripped_name')[ 'BAT Photon Index (15-150 keV) (PL = simple power-law, CPL = cutoff power-law)'].values[0] self.photon_index = self.__clean_string(photon_index) - except IndexError: + except (AttributeError, IndexError): self.photon_index = np.nan def _get_redshift(self) -> None: @@ -265,7 +272,7 @@ def _get_redshift(self) -> None: self.redshift = self.__clean_string(redshift) else: self.redshift = redshift - except IndexError: + except (AttributeError, IndexError): self.redshift = np.nan def _get_redshift_for_luminosity_calculation(self) -> Union[float, None]: @@ -293,7 +300,7 @@ def _set_t90(self) -> None: if t90 == 0.: return np.nan self.t90 = self.__clean_string(t90) - except IndexError: + except (AttributeError, IndexError): self.t90 = np.nan @staticmethod @@ -368,7 +375,7 @@ def _convert_flux_to_luminosity( self.x, self.x_err, self.y, self.y_err = converter.convert_flux_to_luminosity() self._save_luminosity_data() - def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k') -> matplotlib.axes.Axes: + def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k', **kwargs: dict) -> matplotlib.axes.Axes: """ Plots the Afterglow lightcurve and returns Axes. @@ -378,41 +385,26 @@ def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k') -> mat Matplotlib axes to plot the lightcurve into. Useful for user specific modifications to the plot. colour: str, optional Colour of the data. + kwargs: dict + Additional keyword arguments to pass in the Plotter methods. Returns ---------- matplotlib.axes.Axes: The axes with the plot. """ - x_err = [np.abs(self.x_err[1, :]), self.x_err[0, :]] - y_err = [np.abs(self.y_err[1, :]), self.y_err[0, :]] - - ax = axes or plt.gca() - ax.errorbar(self.x, self.y, xerr=x_err, yerr=y_err, - fmt='x', c=colour, ms=1, elinewidth=2, capsize=0.) - - ax.set_xscale('log') - ax.set_yscale('log') - - ax.set_xlim(0.5 * self.x[0], 2 * (self.x[-1] + x_err[1][-1])) - ax.set_ylim(0.5 * min(self.y), 2. * np.max(self.y)) - - ax.annotate(self.name, xy=(0.95, 0.9), xycoords='axes fraction', - horizontalalignment='right', size=20) - - ax.set_xlabel(r'Time since burst [s]') - ax.set_ylabel(self.ylabel) - ax.tick_params(axis='x', pad=10) - - if axes is None: - plt.tight_layout() + if self.flux_data: + plotter = IntegratedFluxPlotter(transient=self) + elif self.luminosity_data: + plotter = LuminosityPlotter(transient=self) + elif self.flux_density_data: + plotter = FluxDensityPlotter(transient=self) + elif self.magnitude_data: + plotter = MagnitudePlotter(transient=self) + else: + return axes + return plotter.plot_data(axes=axes, colour=colour, **kwargs) - grb_dir, _, _ = afterglow_directory_structure(grb=self._stripped_name, data_mode=self.data_mode) - filename = f"{self.name}_lc.png" - plt.savefig(join(grb_dir, filename)) - if axes is None: - plt.clf() - return ax def plot_multiband( self, figure: matplotlib.figure.Figure = None, axes: matplotlib.axes.Axes = None, ncols: int = 2, diff --git a/redback/transient/transient.py b/redback/transient/transient.py index df1e45bb8..6beacbb02 100644 --- a/redback/transient/transient.py +++ b/redback/transient/transient.py @@ -8,7 +8,8 @@ from typing import Union import redback -from redback.plotting import MultiBandPlotter +from redback.plotting import \ + MultiBandPlotter, LuminosityPlotter, FluxDensityPlotter, IntegratedFluxPlotter, MagnitudePlotter class Transient(object): @@ -805,49 +806,18 @@ def plot_data( ylabel: Plot ylabel plot_label: Additional filename label appended to the default name """ - if filters is None: - filters = self.active_bands - - errorbar_fmt = plot_kwargs.get("errorbar_fmt", "x") - colors = plot_kwargs.get("colors", self.get_colors(filters)) - xlabel = plot_kwargs.get("xlabel", self.xlabel) - ylabel = plot_kwargs.get("ylabel", self.ylabel) - plot_label = plot_kwargs.get("plot_label", "data") - - ax = axes or plt.gca() - for indices, band in zip(self.list_of_band_indices, self.unique_bands): - x_err = self.x_err[indices] if self.x_err is not None else self.x_err - if band in filters: - color = colors[filters.index(band)] - label = band - elif plot_others: - color = "black" - label = None - else: - continue - ax.errorbar(self.x[indices], self.y[indices], xerr=x_err, yerr=self.y_err[indices], - fmt=errorbar_fmt, ms=1, color=color, elinewidth=2, capsize=0., label=label) - - ax.set_xlim(0.5 * self.x[0], 1.2 * self.x[-1]) - if self.magnitude_data: - ax.set_ylim(0.8 * min(self.y), 1.2 * np.max(self.y)) - ax.invert_yaxis() + if self.flux_data: + plotter = IntegratedFluxPlotter(transient=self) + elif self.luminosity_data: + plotter = LuminosityPlotter(transient=self) + elif self.flux_density_data: + plotter = FluxDensityPlotter(transient=self) + elif self.magnitude_data: + plotter = MagnitudePlotter(transient=self) else: - ax.set_ylim(0.5 * min(self.y), 2. * np.max(self.y)) - ax.set_yscale('log') - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel) - ax.tick_params(axis='x', pad=10) - ax.legend(ncol=2) + return axes + return plotter.plot_data(axes=axes, filters=filters, plot_others=plot_others, plot_save=plot_save, **plot_kwargs) - if axes is None: - plt.tight_layout() - - if plot_save: - filename = f"{self.name}_{self.data_mode}_{plot_label}.png" - plt.savefig(join(self.transient_dir, filename), bbox_inches='tight') - plt.clf() - return axes def plot_multiband( self, figure: matplotlib.figure.Figure = None, axes: matplotlib.axes.Axes = None, ncols: int = 2, From 314650334601184e694e9ed82bb46726eb380541 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 20:46:03 +1100 Subject: [PATCH 03/37] Added simple prior loading test --- redback/priors/cocoon.prior | 2 +- redback/priors/cone_afterglow.prior | 4 ++-- .../exponential_powerlaw_bolometric.prior | 6 +++--- redback/priors/gaussian.prior | 4 ++-- redback/priors/gaussiancore.prior | 4 ++-- redback/priors/kilonova_afterglow.prior | 2 +- redback/priors/powerlawcore.prior | 4 ++-- redback/priors/slsn.prior | 4 ++-- redback/priors/slsn_bolometric.prior | 4 ++-- redback/priors/smoothpowerlaw.prior | 4 ++-- redback/priors/sn_exponential_powerlaw.prior | 2 +- redback/priors/tophat.prior | 4 ++-- test/prior_test.py | 21 +++++++++++++++++++ 13 files changed, 43 insertions(+), 22 deletions(-) create mode 100644 test/prior_test.py diff --git a/redback/priors/cocoon.prior b/redback/priors/cocoon.prior index 98463c728..62bf039b1 100644 --- a/redback/priors/cocoon.prior +++ b/redback/priors/cocoon.prior @@ -8,6 +8,6 @@ logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') logepse = Uniform(-5, 0, 'logepse', latex_label=r'$\log_{10} \epsilon_{e}$') logepsb = Uniform(-5, 0, 'logepsb', latex_label=r'$\log_{10} \epsilon_{B}$') -ksin = Uniform(0., 1., 'ksin', latex_label=r'$\xi_{N}$') +ksin = Uniform(0., 1., 'ksin', latex_label=r'$\\xi_{N}$') g0 = Uniform(1,100, 'g0', latex_label=r'$\Gamma_{0}$') #max_energy = Constraint(minimum=1e50,maximum=6e51) diff --git a/redback/priors/cone_afterglow.prior b/redback/priors/cone_afterglow.prior index 31fbc637c..6cb0b5eb0 100644 --- a/redback/priors/cone_afterglow.prior +++ b/redback/priors/cone_afterglow.prior @@ -1,5 +1,5 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thV = Cosine('thV', latex_label=r'$\theta_{\mathrm{observer}}$') +thV = Cosine(name='thV', latex_label=r'$\theta_{\mathrm{observer}}$') logE0 = Uniform(44, 54, 'logE0', latex_label=r'$\log_{10} E_{0}$') thC = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') thW = Uniform(1, 8, 'thW', latex_label=r'$\theta_{\mathrm{truncation}}$') @@ -8,5 +8,5 @@ logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') logepse = Uniform(-5, 0, 'logepse', latex_label=r'$\log_{10} \epsilon_{e}$') logepsb = Uniform(-5, 0, 'logepsb', latex_label=r'$\log_{10} \epsilon_{B}$') -ksin = Uniform(0., 1., 'ksin', latex_label=r'$\xi_{N}$') +ksin = Uniform(0., 1., 'ksin', latex_label=r'$\\xi_{N}$') g0 = Uniform(100,2000, 'g0', latex_label=r'$\Gamma_{0}$') diff --git a/redback/priors/exponential_powerlaw_bolometric.prior b/redback/priors/exponential_powerlaw_bolometric.prior index e7c7cb771..58366511a 100644 --- a/redback/priors/exponential_powerlaw_bolometric.prior +++ b/redback/priors/exponential_powerlaw_bolometric.prior @@ -1,6 +1,6 @@ -lbol_0 = LogUniform(1e36,1e48,name='lbol_0', latex_label = r'$L_{\mathrm{bol},0}}$' -alpha_1 = Uniform(0,10,name='alpha_1',latex_label='$\alpha_{1}$') -alpha_2 = Uniform(0,10,name='alpha_2',latex_label='$\alpha_{2}$') +lbol_0 = LogUniform(1e36,1e48,name='lbol_0', latex_label = r'$L_{\mathrm{bol},0}$') +alpha_1 = Uniform(0,10,name='alpha_1',latex_label='$\\alpha_{1}$') +alpha_2 = Uniform(0,10,name='alpha_2',latex_label='$\\alpha_{2}$') tpeak_d = LogUniform(0.001,200,name='tpeak', latex_label = r'$t_{\mathrm{peak}}$~[days]') mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$') vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') diff --git a/redback/priors/gaussian.prior b/redback/priors/gaussian.prior index 0237b20f4..def6b406e 100644 --- a/redback/priors/gaussian.prior +++ b/redback/priors/gaussian.prior @@ -1,5 +1,5 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thv = Sine('thV', latex_label=r'$\theta_{\mathrm{observer}}$') +thv = Sine(name='thV', latex_label=r'$\theta_{\mathrm{observer}}$') loge0 = Uniform(44, 54, 'loge0', latex_label=r'$\log_{10} E_{0}$') thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') thw = Uniform(1, 8, 'thw', latex_label=r'$\theta_{\mathrm{truncation}}$') @@ -7,5 +7,5 @@ logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') logepse = Uniform(-5, 0, 'logepse', latex_label=r'$\log_{10} \epsilon_{e}$') logepsb = Uniform(-5, 0, 'logepsb', latex_label=r'$\log_{10} \epsilon_{B}$') -ksin = Uniform(0., 1., 'ksin', latex_label=r'$\xi_{N}$') +ksin = Uniform(0., 1., 'ksin', latex_label=r'$\\xi_{N}$') g0 = Uniform(100,2000, 'g0', latex_label=r'$\Gamma_{0}$') diff --git a/redback/priors/gaussiancore.prior b/redback/priors/gaussiancore.prior index 5e194128c..94b03fe56 100644 --- a/redback/priors/gaussiancore.prior +++ b/redback/priors/gaussiancore.prior @@ -1,5 +1,5 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thv = Sine('thv', latex_label=r'$\theta_{\mathrm{observer}}$') +thv = Sine(name='thv', latex_label=r'$\theta_{\mathrm{observer}}$') loge0 = Uniform(44, 54, 'loge0', latex_label=r'$\log_{10} E_{0}$') thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') thw = Uniform(1, 8, 'thw', latex_label=r'$\theta_{\mathrm{truncation}}$') @@ -7,5 +7,5 @@ logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') logepse = Uniform(-5, 0, 'logepse', latex_label=r'$\log_{10} \epsilon_{e}$') logepsb = Uniform(-5, 0, 'logepsb', latex_label=r'$\log_{10} \epsilon_{B}$') -ksin = Uniform(0., 1., 'ksin', latex_label=r'$\xi_{N}$') +ksin = Uniform(0., 1., 'ksin', latex_label=r'$\\xi_{N}$') g0 = Uniform(100,2000, 'g0', latex_label=r'$\Gamma_{0}$') diff --git a/redback/priors/kilonova_afterglow.prior b/redback/priors/kilonova_afterglow.prior index 353bba271..eb1ee0212 100644 --- a/redback/priors/kilonova_afterglow.prior +++ b/redback/priors/kilonova_afterglow.prior @@ -8,6 +8,6 @@ logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') logepse = Uniform(-5, 0, 'logepse', latex_label=r'$\log_{10} \epsilon_{e}$') logepsb = Uniform(-5, 0, 'logepsb', latex_label=r'$\log_{10} \epsilon_{B}$') -ksin = Uniform(0., 1., 'ksin', latex_label=r'$\xi_{N}$') +ksin = Uniform(0., 1., 'ksin', latex_label=r'$\\xi_{N}$') g0 = Uniform(1,10, 'g0', latex_label=r'$\Gamma_{0}$') #max_energy = Constraint(minimum=1e50,maximum=6e51) diff --git a/redback/priors/powerlawcore.prior b/redback/priors/powerlawcore.prior index 46c90e1b7..ec9f832a8 100644 --- a/redback/priors/powerlawcore.prior +++ b/redback/priors/powerlawcore.prior @@ -1,5 +1,5 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thv = Sine('thv', latex_label=r'$\theta_{\mathrm{observer}}$') +thv = Sine(name='thv', latex_label=r'$\theta_{\mathrm{observer}}$') loge0 = Uniform(44, 54, 'loge0', latex_label=r'$\log_{10} E_{0}$') thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') thw = Uniform(1, 8, 'thw', latex_label=r'$\theta_{\mathrm{truncation}}$') @@ -8,5 +8,5 @@ logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') logepse = Uniform(-5, 0, 'logepse', latex_label=r'$\log_{10} \epsilon_{e}$') logepsb = Uniform(-5, 0, 'logepsb', latex_label=r'$\log_{10} \epsilon_{B}$') -ksin = Uniform(0., 1., 'ksin', latex_label=r'$\xi_{N}$') +ksin = Uniform(0., 1., 'ksin', latex_label=r'$\\xi_{N}$') g0 = Uniform(100,2000, 'g0', latex_label=r'$\Gamma_{0}$') diff --git a/redback/priors/slsn.prior b/redback/priors/slsn.prior index 24f11aefd..444cb29d0 100644 --- a/redback/priors/slsn.prior +++ b/redback/priors/slsn.prior @@ -8,5 +8,5 @@ vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') kappa_gamma = Uniform(1e-4, 1e4, 'kappa_gamma', latex_label = r'$\kappa_{\gamma}$') temperature_floor = LogUniform(1e3,1e5,name = 'temperature_floor', latex_label = r'$T_{\mathrm{floor}$ [k]}') -e_rot_constraint = Constraint(name='e_rot_constraint', mininum=10, maximum=1e10) -t_nebula_min = Constraint(name='t_nebula_min', mininum=0.1,maximum=500) \ No newline at end of file +e_rot_constraint = Constraint(name='e_rot_constraint', minimum=10, maximum=1e10) +t_nebula_min = Constraint(name='t_nebula_min', minimum=0.1,maximum=500) \ No newline at end of file diff --git a/redback/priors/slsn_bolometric.prior b/redback/priors/slsn_bolometric.prior index 5c689a108..2eb9e6b19 100644 --- a/redback/priors/slsn_bolometric.prior +++ b/redback/priors/slsn_bolometric.prior @@ -6,5 +6,5 @@ mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$ vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') kappa_gamma = Uniform(1e-4, 1e4, 'kappa_gamma', latex_label = r'$\kappa_{\gamma}$') -e_rot_constraint = Constraint(name='e_rot_constraint', mininum=10, maximum=1e10) -t_nebula_min = Constraint(name='t_nebula_min', mininum=0.1,maximum=500) \ No newline at end of file +e_rot_constraint = Constraint(name='e_rot_constraint', minimum=10, maximum=1e10) +t_nebula_min = Constraint(name='t_nebula_min', minimum=0.1,maximum=500) \ No newline at end of file diff --git a/redback/priors/smoothpowerlaw.prior b/redback/priors/smoothpowerlaw.prior index 46c90e1b7..ec9f832a8 100644 --- a/redback/priors/smoothpowerlaw.prior +++ b/redback/priors/smoothpowerlaw.prior @@ -1,5 +1,5 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thv = Sine('thv', latex_label=r'$\theta_{\mathrm{observer}}$') +thv = Sine(name='thv', latex_label=r'$\theta_{\mathrm{observer}}$') loge0 = Uniform(44, 54, 'loge0', latex_label=r'$\log_{10} E_{0}$') thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') thw = Uniform(1, 8, 'thw', latex_label=r'$\theta_{\mathrm{truncation}}$') @@ -8,5 +8,5 @@ logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') logepse = Uniform(-5, 0, 'logepse', latex_label=r'$\log_{10} \epsilon_{e}$') logepsb = Uniform(-5, 0, 'logepsb', latex_label=r'$\log_{10} \epsilon_{B}$') -ksin = Uniform(0., 1., 'ksin', latex_label=r'$\xi_{N}$') +ksin = Uniform(0., 1., 'ksin', latex_label=r'$\\xi_{N}$') g0 = Uniform(100,2000, 'g0', latex_label=r'$\Gamma_{0}$') diff --git a/redback/priors/sn_exponential_powerlaw.prior b/redback/priors/sn_exponential_powerlaw.prior index 22702fa70..b2b059ae4 100644 --- a/redback/priors/sn_exponential_powerlaw.prior +++ b/redback/priors/sn_exponential_powerlaw.prior @@ -1,5 +1,5 @@ redshift = Uniform(1e-6,3,name='redshift', latex_label = r'$z$') -lbol_0 = LogUniform(1e36,1e48,name='lbol_0', latex_label = r'$L_{\mathrm{bol},0}}$' +lbol_0 = LogUniform(1e36,1e48,name='lbol_0', latex_label = r'$L_{\mathrm{bol},0}$') alpha_1 = Uniform(0,10,name='alpha_1',latex_label='$\alpha_{1}$') alpha_2 = Uniform(0,10,name='alpha_2',latex_label='$\alpha_{2}$') tpeak_d = LogUniform(0.001,200,name='tpeak', latex_label = r'$t_{\mathrm{peak}}$~[days]') diff --git a/redback/priors/tophat.prior b/redback/priors/tophat.prior index 89e13f1f6..73d1f66c2 100644 --- a/redback/priors/tophat.prior +++ b/redback/priors/tophat.prior @@ -1,10 +1,10 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thv = Sine('thv', latex_label=r'$\theta_{\mathrm{observer}}$') +thv = Sine(name='thv', latex_label=r'$\theta_{\mathrm{observer}}$') loge0 = Uniform(44, 54, 'loge0', latex_label=r'$\log_{10} E_{0}$') thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') logepse = Uniform(-5, 0, 'logepse', latex_label=r'$\log_{10} \epsilon_{e}$') logepsb = Uniform(-5, 0, 'logepsb', latex_label=r'$\log_{10} \epsilon_{B}$') -ksin = Uniform(0., 1., 'ksin', latex_label=r'$\xi_{N}$') +ksin = Uniform(0., 1., 'ksin', latex_label=r'$\\xi_{N}$') g0 = Uniform(100,2000, 'g0', latex_label=r'$\Gamma_{0}$') diff --git a/test/prior_test.py b/test/prior_test.py new file mode 100644 index 000000000..494de699b --- /dev/null +++ b/test/prior_test.py @@ -0,0 +1,21 @@ +import unittest +import bilby +from os import listdir + + +class TestLoadPriors(unittest.TestCase): + + def setUp(self) -> None: + pass + + def tearDown(self) -> None: + pass + + def test_load_priors(self): + path_to_files = "../redback/priors/" + prior_files = listdir(path_to_files) + + for f in prior_files: + print(f) + prior_dict = bilby.prior.PriorDict() + prior_dict.from_file(f"{path_to_files}{f}") From 351334097ebbe7fa4ddc9282a643c67ea5b393ca Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 20:59:45 +1100 Subject: [PATCH 04/37] Fixed broadband afterglow private data example --- ...roadband_afterglow_private_data_example.py | 19 ++++++++++--------- redback/plotting.py | 5 ++++- redback/priors/gaussiancore.prior | 2 +- redback/transient/transient.py | 4 +++- redback/transient_models/afterglow_models.py | 8 ++++---- 5 files changed, 22 insertions(+), 16 deletions(-) diff --git a/examples/broadband_afterglow_private_data_example.py b/examples/broadband_afterglow_private_data_example.py index 738faf0ba..ee3f85dd6 100644 --- a/examples/broadband_afterglow_private_data_example.py +++ b/examples/broadband_afterglow_private_data_example.py @@ -8,10 +8,10 @@ # load the data file data = pd.read_csv('example_data/grb_afterglow.csv') -time_d = data['time'] -flux_density = data['flux'] -frequency = data['frequency'] -flux_density_err = data['flux_err'] +time_d = data['time'].values +flux_density = data['flux'].values +frequency = data['frequency'].values +flux_density_err = data['flux_err'].values # we now load the afterglow transient object. We are using flux_density data here so we need to use that data mode data_mode = 'flux_density' @@ -20,14 +20,15 @@ name = '170817A' redshift = 1e-2 -afterglow = redback.transient.Afterglow(name=name, data_mode=data_mode, time=time_d, - flux_density=flux_density, flux_err=flux_density_err, frequency=frequency) +afterglow = redback.transient.Afterglow( + name=name, data_mode=data_mode, time=time_d, + flux_density=flux_density, flux_density_err=flux_density_err, frequency=frequency) # Now we have loaded the data up, we can plot it. afterglow.plot_data() # now let's actually fit it with data. We will use all the data and a gaussiancore structured jet from afterglowpy. -# Note this is not a fast example so we will make some sampling sacrifices for speed. +# Note this is not a fast example, so we will make some sampling sacrifices for speed. model = 'gaussiancore' @@ -43,11 +44,11 @@ priors['logepsb'] = -3.8 priors['ksin'] = 1. -model_kwargs = dict(frequencies=frequency, output_format='flux_density') +model_kwargs = dict(frequency=frequency, output_format='flux_density') # returns a supernova result object result = redback.fit_model(name=name, transient=afterglow, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, data_mode='flux_density', sample='rslice', nlive=200, resume=False) + prior=priors, data_mode='flux_density', nlive=200, resume=False) # plot corner result.plot_corner() diff --git a/redback/plotting.py b/redback/plotting.py index 1f022f03f..e956abf8d 100644 --- a/redback/plotting.py +++ b/redback/plotting.py @@ -239,7 +239,10 @@ def plot_data( if plot_save: filename = f"{self.transient.name}_{self.transient.data_mode}_{plot_label}.png" - plt.savefig(join(self.transient.transient_dir, filename), bbox_inches='tight') + structure = afterglow_directory_structure( + grb=self.transient.name, data_mode='', instrument='') + + plt.savefig(join(structure.directory_path, filename), bbox_inches='tight') plt.clf() return axes diff --git a/redback/priors/gaussiancore.prior b/redback/priors/gaussiancore.prior index 94b03fe56..0c542d2ac 100644 --- a/redback/priors/gaussiancore.prior +++ b/redback/priors/gaussiancore.prior @@ -1,5 +1,5 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thv = Sine(name='thv', latex_label=r'$\theta_{\mathrm{observer}}$') +thv = Sine(maximum=np.pi/2, name='thv', latex_label=r'$\theta_{\mathrm{observer}}$') loge0 = Uniform(44, 54, 'loge0', latex_label=r'$\log_{10} E_{0}$') thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') thw = Uniform(1, 8, 'thw', latex_label=r'$\theta_{\mathrm{truncation}}$') diff --git a/redback/transient/transient.py b/redback/transient/transient.py index 6beacbb02..0d4cc13eb 100644 --- a/redback/transient/transient.py +++ b/redback/transient/transient.py @@ -354,7 +354,9 @@ def active_bands(self, active_bands: Union[list, str]) -> None: self._active_bands = active_bands @property - def filtered_indices(self) -> list: + def filtered_indices(self) -> Union[list, None]: + if self.bands is None: + return list(np.arange(len(self.x))) return [b in self.active_bands for b in self.bands] def get_filtered_data(self) -> tuple: diff --git a/redback/transient_models/afterglow_models.py b/redback/transient_models/afterglow_models.py index 6f7447695..45f620571 100644 --- a/redback/transient_models/afterglow_models.py +++ b/redback/transient_models/afterglow_models.py @@ -88,11 +88,11 @@ def cone_afterglow(time, redshift, thv, loge0, thw, thc, logn0, p, logepse, loge def gaussiancore(time, redshift, thv, loge0, thc, thw, logn0, p, logepse, logepsb, ksin, g0, **kwargs): time = time * day_to_s dl = cosmo.luminosity_distance(redshift).cgs.value - spread = kwargs['spread'] - latres = kwargs['latres'] - tres = kwargs['tres'] + spread = kwargs.get('spread', True) + latres = kwargs.get('latres', 2) + tres = kwargs.get('tres', 100) jettype = jettype_dict['gaussian_w_core'] - spectype = kwargs['spectype'] + spectype = kwargs.get('spectype', 0) frequency = kwargs['frequency'] thw = thw * thc From 78f80b30b1906045931109450b78a2da77a5706a Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 21:00:05 +1100 Subject: [PATCH 05/37] reformatted --- examples/fit_your_own_model_example.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/fit_your_own_model_example.py b/examples/fit_your_own_model_example.py index ecc6ab774..bcb8e9797 100644 --- a/examples/fit_your_own_model_example.py +++ b/examples/fit_your_own_model_example.py @@ -1,6 +1,7 @@ import redback from bilby.core.prior import LogUniform, Uniform + # If your favourite model is not implemented in redback. You can still fit it using redback! # Now instead of passing a string as the model. You need to pass a python function. @@ -9,7 +10,8 @@ # time must be the first element. def my_favourite_model(time, l0, alpha): - return l0*time**alpha + return l0 * time ** alpha + model = my_favourite_model @@ -31,8 +33,8 @@ def my_favourite_model(time, l0, alpha): # You need to create your own priors for this new model. # The model has two parameters l0 and alpha. We use bilby priors for this priors = {} -priors['l0'] = LogUniform(1e40, 1e55, 'l0', latex_label = r'$l_{0}$') -priors['alpha_1'] = Uniform(-7, -1, 'alpha_1', latex_label = r'$\alpha_{1}$') +priors['l0'] = LogUniform(1e40, 1e55, 'l0', latex_label=r'$l_{0}$') +priors['alpha_1'] = Uniform(-7, -1, 'alpha_1', latex_label=r'$\alpha_{1}$') # Call redback.fit_model to run the sampler and obtain GRB result object result = redback.fit_model(name=GRB, model=model, sampler='dynesty', nlive=200, transient=afterglow, From 914d9cdb966756141682ce952c221c65f2176ece Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 22:00:35 +1100 Subject: [PATCH 06/37] Fixed issues with broadband_afterglow_private_data_example.py --- ...roadband_afterglow_private_data_example.py | 5 ++-- redback/plotting.py | 30 +++++++++++-------- redback/result.py | 16 ++++++++++ redback/sampler.py | 14 ++++----- redback/transient/transient.py | 16 ++++++++-- 5 files changed, 56 insertions(+), 25 deletions(-) diff --git a/examples/broadband_afterglow_private_data_example.py b/examples/broadband_afterglow_private_data_example.py index ee3f85dd6..c9fa8bcc0 100644 --- a/examples/broadband_afterglow_private_data_example.py +++ b/examples/broadband_afterglow_private_data_example.py @@ -26,6 +26,7 @@ # Now we have loaded the data up, we can plot it. afterglow.plot_data() +afterglow.plot_multiband() # now let's actually fit it with data. We will use all the data and a gaussiancore structured jet from afterglowpy. # Note this is not a fast example, so we will make some sampling sacrifices for speed. @@ -33,7 +34,7 @@ model = 'gaussiancore' # use default priors and 'nestle' sampler -sampler = 'nestle' +sampler = 'dynesty' priors = redback.priors.get_priors(model=model) priors['redshift'] = redshift @@ -48,7 +49,7 @@ # returns a supernova result object result = redback.fit_model(name=name, transient=afterglow, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, data_mode='flux_density', nlive=200, resume=False) + prior=priors, sample='rslice', data_mode='flux_density', nlive=200, resume=True) # plot corner result.plot_corner() diff --git a/redback/plotting.py b/redback/plotting.py index e956abf8d..70b779067 100644 --- a/redback/plotting.py +++ b/redback/plotting.py @@ -75,12 +75,11 @@ def plot_multiband( f"but {len(filters)} panels are needed.") if figsize is None: figsize = (4 * ncols, 4 * nrows) - figure, axes = plt.subplots(ncols=ncols, nrows=nrows, sharex='all', sharey='all', figsize=figsize) + figure, axes = plt.subplots(ncols=ncols, nrows=nrows, sharex='all', figsize=figsize) axes = axes.ravel() - i = 0 - for indices, band in zip(self.transient.list_of_band_indices, self.transient.unique_bands): + for i, (indices, band) in enumerate(zip(self.transient.list_of_band_indices, self.transient.unique_bands)): if band not in filters: continue @@ -88,11 +87,14 @@ def plot_multiband( color = colors[filters.index(band)] - freq = self.transient.bands_to_frequencies([band]) - if 1e10 < freq < 1e15: - label = band + if isinstance(band, str): + freq = self.transient.bands_to_frequencies([band]) + if 1e10 < freq < 1e15: + label = band + else: + label = freq else: - label = freq + label = f"{band:.2e}" axes[i].errorbar(self.transient.x[indices], self.transient.y[indices], xerr=x_err, yerr=self.transient.y_err[indices], fmt=errorbar_fmt, ms=1, color=color, elinewidth=2, capsize=0., label=label) @@ -105,13 +107,13 @@ def plot_multiband( axes[i].set_yscale("log") axes[i].legend(ncol=2) axes[i].tick_params(axis='both', which='major', pad=8) - i += 1 figure.supxlabel(xlabel, fontsize=fontsize) figure.supylabel(ylabel, fontsize=fontsize) - filename = f"{self.transient.name}_{self.transient.data_mode}_{plot_label}.png" + filename = f"{self.transient.name}_{plot_label}.png" plt.subplots_adjust(wspace=wspace, hspace=hspace) - plt.savefig(join(self.transient.transient_dir, filename), bbox_inches="tight") + directory_structure = afterglow_directory_structure(grb=self.transient.name, data_mode=self.transient.data_mode) + plt.savefig(join(directory_structure.directory_path, filename), bbox_inches="tight") return axes @@ -223,6 +225,8 @@ def plot_data( label = None else: continue + if isinstance(label, float): + label = f"{label:.2e}" ax.errorbar( self.transient.x[indices], self.transient.y[indices], xerr=x_err, yerr=self.transient.y_err[indices], fmt=errorbar_fmt, ms=1, color=color, elinewidth=2, capsize=0., label=label) @@ -232,15 +236,15 @@ def plot_data( ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.tick_params(axis='x', pad=10) - ax.legend(ncol=2) + ax.legend(ncol=2, loc='best') if axes is None: plt.tight_layout() if plot_save: - filename = f"{self.transient.name}_{self.transient.data_mode}_{plot_label}.png" + filename = f"{self.transient.name}_{plot_label}.png" structure = afterglow_directory_structure( - grb=self.transient.name, data_mode='', instrument='') + grb=self.transient.name, data_mode=self.transient.data_mode, instrument='') plt.savefig(join(structure.directory_path, filename), bbox_inches='tight') plt.clf() diff --git a/redback/result.py b/redback/result.py index 7719855f2..e38ab490a 100644 --- a/redback/result.py +++ b/redback/result.py @@ -147,6 +147,22 @@ def plot_lightcurve(self, model: Union[callable, str] = None, **kwargs: dict) -> self.transient.plot_lightcurve(model=model, posterior=self.posterior, outdir=self.outdir, model_kwargs=self.model_kwargs, **kwargs) + def plot_multiband_lightcurve(self, model: Union[callable, str] = None, **kwargs: dict) -> None: + """ + Reconstructs the transient and calls the specific `plot_multiband_lightcurve` method. + + Parameters + ---------- + model: Union[callable, str], optional + User specified model. + kwargs: dict + Any kwargs to be passed into the `plot_lightcurve` method. + """ + if model is None: + model = model_library.all_models_dict[self.model] + self.transient.plot_multiband_lightcurve( + model=model, posterior=self.posterior, outdir=self.outdir, model_kwargs=self.model_kwargs, **kwargs) + def plot_data(self, **kwargs: dict) -> None: """ Reconstructs the transient and calls the specific `plot_data` method. diff --git a/redback/sampler.py b/redback/sampler.py index ae5c07047..dde5dc4e9 100644 --- a/redback/sampler.py +++ b/redback/sampler.py @@ -97,10 +97,9 @@ def _fit_grb(transient, model, outdir=None, label=None, sampler='dynesty', nlive outdir = f"{outdir}/{model.__name__}" Path(outdir).mkdir(parents=True, exist_ok=True) - if label is None: - label = transient.data_mode - if use_photon_index_prior: - label += '_photon_index' + label = kwargs.get("label", transient.name) + if use_photon_index_prior: + label += '_photon_index' if transient.flux_density_data or transient.magnitude_data: x, x_err, y, y_err = transient.get_filtered_data() @@ -122,7 +121,7 @@ def _fit_grb(transient, model, outdir=None, label=None, sampler='dynesty', nlive return result -def _fit_kilonova(transient, model, outdir=None, label=None, sampler='dynesty', nlive=3000, prior=None, walks=1000, +def _fit_kilonova(transient, model, outdir=None, sampler='dynesty', nlive=3000, prior=None, walks=1000, resume=True, save_format='json', model_kwargs=None, **kwargs): if outdir is None: @@ -131,8 +130,7 @@ def _fit_kilonova(transient, model, outdir=None, label=None, sampler='dynesty', outdir = f"{outdir}/{model.__name__}" Path(outdir).mkdir(parents=True, exist_ok=True) - if label is None: - label = f"{transient.data_mode}" + label = kwargs.get("label", transient.name) if transient.flux_density_data or transient.magnitude_data: x, x_err, y, y_err = transient.get_filtered_data() @@ -162,7 +160,7 @@ def _fit_prompt(name, transient, model, outdir, integrated_rate_function=True, s outdir = f"{outdir}/GRB{name}/{model.__name__}" Path(outdir).mkdir(parents=True, exist_ok=True) - label = transient.data_mode + label = kwargs.get("label", transient.name) if use_photon_index_prior: label += '_photon_index' diff --git a/redback/transient/transient.py b/redback/transient/transient.py index 0d4cc13eb..3b598c742 100644 --- a/redback/transient/transient.py +++ b/redback/transient/transient.py @@ -132,8 +132,8 @@ def __init__( self.counts_err = np.sqrt(counts) if counts is not None else None self.ttes = ttes - self.bands = bands self.frequency = frequency + self.bands = bands self.system = system self.active_bands = active_bands self.data_mode = data_mode @@ -339,7 +339,7 @@ def active_bands(self) -> list: return self._active_bands @active_bands.setter - def active_bands(self, active_bands: Union[list, str]) -> None: + def active_bands(self, active_bands: Union[list, str, None]) -> None: """ Parameters @@ -359,6 +359,18 @@ def filtered_indices(self) -> Union[list, None]: return list(np.arange(len(self.x))) return [b in self.active_bands for b in self.bands] + @property + def bands(self) -> Union[list, None, np.ndarray]: + return self._bands + + @bands.setter + def bands(self, bands: Union[list, None, np.ndarray]): + if bands is None: + self._bands = self.frequency + else: + self._bands = bands + + def get_filtered_data(self) -> tuple: """ Used to filter flux density and photometry data, so we only use data that is using the active bands. From 35c2d5ae926ce4d6a7f7463d83867afc95c631e9 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 22:02:33 +1100 Subject: [PATCH 07/37] Moved a logger message --- redback/get_data/directory.py | 2 -- redback/get_data/swift.py | 9 +++++++++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/redback/get_data/directory.py b/redback/get_data/directory.py index 7cf83206a..0932a7bbc 100644 --- a/redback/get_data/directory.py +++ b/redback/get_data/directory.py @@ -40,11 +40,9 @@ def afterglow_directory_structure(grb: str, data_mode: str, instrument: str = 'B if instrument == 'XRT': raw_file_path = f'{path}_xrt_rawSwiftData.csv' processed_file_path = f'{path}_xrt.csv' - logger.warning('You are only downloading XRT data, you may not capture the tail of the prompt emission.') else: raw_file_path = f'{path}_rawSwiftData.csv' processed_file_path = f'{path}.csv' - logger.warning('You are downloading BAT and XRT data, you will need to truncate the data for some models.') return DirectoryStructure( directory_path=directory_path, raw_file_path=raw_file_path, processed_file_path=processed_file_path) diff --git a/redback/get_data/swift.py b/redback/get_data/swift.py index 4afe6cd5d..076c5eaaf 100644 --- a/redback/get_data/swift.py +++ b/redback/get_data/swift.py @@ -176,6 +176,15 @@ def get_data(self) -> pd.DataFrame: ------- pandas.DataFrame: The processed data. """ + if self.instrument == "BAT+XRT": + logger.warning( + "You are downloading BAT and XRT data, " + "you will need to truncate the data for some models.") + elif self.instrument == "XRT": + logger.warning( + "You are only downloading XRT data, you may not capture" + " the tail of the prompt emission.") + self.collect_data() return self.convert_raw_data_to_csv() From fff54cc0c993a4ae7d024b769de2625198b4735f Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 22:18:53 +1100 Subject: [PATCH 08/37] fixed fit_your_own_model_example.py --- examples/fit_your_own_model_example.py | 16 +++++---- redback/sampler.py | 2 +- redback/transient/transient.py | 46 +++++++++++++++----------- 3 files changed, 37 insertions(+), 27 deletions(-) diff --git a/examples/fit_your_own_model_example.py b/examples/fit_your_own_model_example.py index bcb8e9797..c67c621bf 100644 --- a/examples/fit_your_own_model_example.py +++ b/examples/fit_your_own_model_example.py @@ -1,3 +1,5 @@ +import bilby.core.prior + import redback from bilby.core.prior import LogUniform, Uniform @@ -9,7 +11,7 @@ # You can make up any model you like. # time must be the first element. -def my_favourite_model(time, l0, alpha): +def my_favourite_model(time, l0, alpha, **kwargs): return l0 * time ** alpha @@ -26,18 +28,20 @@ def my_favourite_model(time, l0, alpha): afterglow = redback.afterglow.SGRB.from_swift_grb(name=GRB, data_mode='flux', truncate=True, truncate_method="prompt_time_error") +afterglow.plot_data() + # uses an analytical k-correction expression to create luminosity data if not already there. # Can also use a numerical k-correction through CIAO afterglow.analytical_flux_to_luminosity() # You need to create your own priors for this new model. # The model has two parameters l0 and alpha. We use bilby priors for this -priors = {} -priors['l0'] = LogUniform(1e40, 1e55, 'l0', latex_label=r'$l_{0}$') -priors['alpha_1'] = Uniform(-7, -1, 'alpha_1', latex_label=r'$\alpha_{1}$') +priors = bilby.core.prior.PriorDict() +priors['l0'] = LogUniform(1e-10, 1e5, 'l0', latex_label=r'$l_{0}$') +priors['alpha'] = Uniform(-7, 0, 'alpha', latex_label=r'$\alpha$') # Call redback.fit_model to run the sampler and obtain GRB result object result = redback.fit_model(name=GRB, model=model, sampler='dynesty', nlive=200, transient=afterglow, - prior=priors, data_mode='luminosity', sample='rslice') + prior=priors, data_mode='luminosity', sample='rslice', resume=False) -result.plot_lightcurve(random_models=1000) +result.plot_lightcurve(random_models=1000, model=my_favourite_model) diff --git a/redback/sampler.py b/redback/sampler.py index dde5dc4e9..3077226b9 100644 --- a/redback/sampler.py +++ b/redback/sampler.py @@ -117,7 +117,7 @@ def _fit_grb(transient, model, outdir=None, label=None, sampler='dynesty', nlive outdir=outdir, plot=True, use_ratio=False, walks=walks, resume=resume, maxmcmc=10 * walks, result_class=RedbackResult, meta_data=meta_data, nthreads=4, save_bounds=False, nsteps=nlive, nwalkers=walks, save=save_format, **kwargs) - + plt.close('all') return result diff --git a/redback/transient/transient.py b/redback/transient/transient.py index 3b598c742..cbcc01c7c 100644 --- a/redback/transient/transient.py +++ b/redback/transient/transient.py @@ -132,8 +132,9 @@ def __init__( self.counts_err = np.sqrt(counts) if counts is not None else None self.ttes = ttes - self.frequency = frequency - self.bands = bands + self._frequency = None + self._bands = None + self.set_bands_and_frequency(bands=bands, frequency=frequency) self.system = system self.active_bands = active_bands self.data_mode = data_mode @@ -300,6 +301,18 @@ def ylabel(self) -> str: except KeyError: raise ValueError("No data mode specified") + def set_bands_and_frequency( + self, bands: Union[None, list, np.ndarray], frequency: Union[None, list, np.ndarray]): + if (bands is None and frequency is None) or (bands is not None and frequency is not None): + self._bands = bands + self._frequency = bands + elif bands is None and frequency is not None: + self._frequency = frequency + self._bands = self.frequency + elif bands is not None and frequency is None: + self._bands = bands + self._frequency = self.bands_to_frequencies(self.bands) + @property def frequency(self) -> np.ndarray: """ @@ -319,10 +332,15 @@ def frequency(self, frequency: np.ndarray) -> None: frequency: np.ndarray Set band frequencies if an array is given. Otherwise, convert bands to frequencies. """ - if frequency is None: - self._frequency = self.bands_to_frequencies(self.bands) - else: - self._frequency = frequency + self.set_bands_and_frequency(bands=self.bands, frequency=frequency) + + @property + def bands(self) -> Union[list, None, np.ndarray]: + return self._bands + + @bands.setter + def bands(self, bands: Union[list, None, np.ndarray]): + self.set_bands_and_frequency(bands=bands, frequency=self.frequency) @property def filtered_frequencies(self) -> np.array: @@ -359,18 +377,6 @@ def filtered_indices(self) -> Union[list, None]: return list(np.arange(len(self.x))) return [b in self.active_bands for b in self.bands] - @property - def bands(self) -> Union[list, None, np.ndarray]: - return self._bands - - @bands.setter - def bands(self, bands: Union[list, None, np.ndarray]): - if bands is None: - self._bands = self.frequency - else: - self._bands = bands - - def get_filtered_data(self) -> tuple: """ Used to filter flux density and photometry data, so we only use data that is using the active bands. @@ -520,9 +526,9 @@ def plot_lightcurve( if model_kwargs is None: model_kwargs = dict() axes = axes or plt.gca() - # axes = self.plot_data(axes=axes) + axes = self.plot_data(axes=axes) axes.set_yscale('log') - # plt.semilogy() + plt.semilogy() times = self._get_times(axes) posterior.sort_values(by='log_likelihood') From 82a088e26d71aaa625ccc057f8755d8952720539 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 23:40:15 +1100 Subject: [PATCH 09/37] Fixed more prior issues --- redback/priors/basic_magnetar_powered.prior | 2 +- redback/priors/basic_magnetar_powered_bolometric.prior | 2 +- redback/priors/cone_afterglow.prior | 6 +++--- redback/priors/gaussian.prior | 6 +++--- redback/priors/gaussiancore.prior | 6 +++--- redback/priors/magnetar_nickel.prior | 2 +- redback/priors/powerlawcore.prior | 6 +++--- redback/priors/slsn.prior | 2 +- redback/priors/slsn_bolometric.prior | 2 +- redback/priors/smoothpowerlaw.prior | 6 +++--- redback/priors/tophat.prior | 4 ++-- 11 files changed, 22 insertions(+), 22 deletions(-) diff --git a/redback/priors/basic_magnetar_powered.prior b/redback/priors/basic_magnetar_powered.prior index bbcef5a77..8c22038e6 100644 --- a/redback/priors/basic_magnetar_powered.prior +++ b/redback/priors/basic_magnetar_powered.prior @@ -2,7 +2,7 @@ redshift = Uniform(1e-6,3,name='redshift', latex_label = r'$z$') p0 = Uniform(1, 10, 'p0', latex_label = r'$P_{0]$ [ms]') bp = LogUniform(0.1,10,'bp',latex_label = r'$B_{p}$~[10$^{14}$G]') mass_ns = Uniform(1.1, 2.2, 'mass_ns', latex_label = r'$M_{\mathrm{NS}} [M_{\odot}]$') -theta_pb = Uniform(0, 3.14/2, 'theta_pb', latex_label = r'$\theta_{P-B}$') +theta_pb = Uniform(0, 3.14/2, 'theta_pb', latex_label = r'$\\theta_{P-B}$') mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$') vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') diff --git a/redback/priors/basic_magnetar_powered_bolometric.prior b/redback/priors/basic_magnetar_powered_bolometric.prior index a7ae7225c..0af025c5e 100644 --- a/redback/priors/basic_magnetar_powered_bolometric.prior +++ b/redback/priors/basic_magnetar_powered_bolometric.prior @@ -1,7 +1,7 @@ p0 = Uniform(1, 10, 'p0', latex_label = r'$P_{0]$ [ms]') bp = LogUniform(0.1,10,'bp',latex_label = r'$B_{p}$~[10$^{14}$G]') mass_ns = Uniform(1.1, 2.2, 'mass_ns', latex_label = r'$M_{\mathrm{NS}} [M_{\odot}]$') -theta_pb = Uniform(0, 3.14/2, 'theta_pb', latex_label = r'$\theta_{P-B}$') +theta_pb = Uniform(0, 3.14/2, 'theta_pb', latex_label = r'$\\theta_{P-B}$') mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$') vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') diff --git a/redback/priors/cone_afterglow.prior b/redback/priors/cone_afterglow.prior index 6cb0b5eb0..f9ea05ed1 100644 --- a/redback/priors/cone_afterglow.prior +++ b/redback/priors/cone_afterglow.prior @@ -1,8 +1,8 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thV = Cosine(name='thV', latex_label=r'$\theta_{\mathrm{observer}}$') +thV = Cosine(name='thV', latex_label=r'$\\theta_{\mathrm{observer}}$') logE0 = Uniform(44, 54, 'logE0', latex_label=r'$\log_{10} E_{0}$') -thC = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') -thW = Uniform(1, 8, 'thW', latex_label=r'$\theta_{\mathrm{truncation}}$') +thC = Uniform(0.01, 0.1, 'thc', latex_label=r'$\\theta_{\mathrm{core}}$') +thW = Uniform(1, 8, 'thW', latex_label=r'$\\theta_{\mathrm{truncation}}$') beta = Uniform(0.5, 10, 'beta', latex_label=r'$\\beta$') logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') diff --git a/redback/priors/gaussian.prior b/redback/priors/gaussian.prior index def6b406e..5d753709a 100644 --- a/redback/priors/gaussian.prior +++ b/redback/priors/gaussian.prior @@ -1,8 +1,8 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thv = Sine(name='thV', latex_label=r'$\theta_{\mathrm{observer}}$') +thv = Sine(name='thV', latex_label=r'$\\theta_{\mathrm{observer}}$') loge0 = Uniform(44, 54, 'loge0', latex_label=r'$\log_{10} E_{0}$') -thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') -thw = Uniform(1, 8, 'thw', latex_label=r'$\theta_{\mathrm{truncation}}$') +thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\\theta_{\mathrm{core}}$') +thw = Uniform(1, 8, 'thw', latex_label=r'$\\theta_{\mathrm{truncation}}$') logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') logepse = Uniform(-5, 0, 'logepse', latex_label=r'$\log_{10} \epsilon_{e}$') diff --git a/redback/priors/gaussiancore.prior b/redback/priors/gaussiancore.prior index 0c542d2ac..620bf4475 100644 --- a/redback/priors/gaussiancore.prior +++ b/redback/priors/gaussiancore.prior @@ -1,8 +1,8 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thv = Sine(maximum=np.pi/2, name='thv', latex_label=r'$\theta_{\mathrm{observer}}$') +thv = Sine(maximum=np.pi/2, name='thv', latex_label=r'$\\theta_{\mathrm{observer}}$') loge0 = Uniform(44, 54, 'loge0', latex_label=r'$\log_{10} E_{0}$') -thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') -thw = Uniform(1, 8, 'thw', latex_label=r'$\theta_{\mathrm{truncation}}$') +thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\\theta_{\mathrm{core}}$') +thw = Uniform(1, 8, 'thw', latex_label=r'$\\theta_{\mathrm{truncation}}$') logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') logepse = Uniform(-5, 0, 'logepse', latex_label=r'$\log_{10} \epsilon_{e}$') diff --git a/redback/priors/magnetar_nickel.prior b/redback/priors/magnetar_nickel.prior index 49748be29..40ce59300 100644 --- a/redback/priors/magnetar_nickel.prior +++ b/redback/priors/magnetar_nickel.prior @@ -3,7 +3,7 @@ f_nickel = LogUniform(1e-3,1,name='f_nickel', latex_label = r'$f_{\mathrm{Ni}}') p0 = Uniform(1, 10, 'p0', latex_label = r'$P_{0]$ [ms]') bp = LogUniform(0.1,10,'bp',latex_label = r'$B_{p}$~[10$^{14}$G]') mass_ns = Uniform(1.1, 2.2, 'mass_ns', latex_label = r'$M_{\mathrm{NS}} [M_{\odot}]$') -theta_pb = Uniform(0, 3.14/2, 'theta_pb', latex_label = r'$\theta_{P-B}$') +theta_pb = Uniform(0, 3.14/2, 'theta_pb', latex_label = r'$\\theta_{P-B}$') mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$') vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') diff --git a/redback/priors/powerlawcore.prior b/redback/priors/powerlawcore.prior index ec9f832a8..33559c8b2 100644 --- a/redback/priors/powerlawcore.prior +++ b/redback/priors/powerlawcore.prior @@ -1,8 +1,8 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thv = Sine(name='thv', latex_label=r'$\theta_{\mathrm{observer}}$') +thv = Sine(name='thv', latex_label=r'$\\theta_{\mathrm{observer}}$') loge0 = Uniform(44, 54, 'loge0', latex_label=r'$\log_{10} E_{0}$') -thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') -thw = Uniform(1, 8, 'thw', latex_label=r'$\theta_{\mathrm{truncation}}$') +thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\\theta_{\mathrm{core}}$') +thw = Uniform(1, 8, 'thw', latex_label=r'$\\theta_{\mathrm{truncation}}$') beta = Uniform(0.5, 10, 'beta', latex_label=r'$\\beta$') logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') diff --git a/redback/priors/slsn.prior b/redback/priors/slsn.prior index 444cb29d0..0c97c7345 100644 --- a/redback/priors/slsn.prior +++ b/redback/priors/slsn.prior @@ -2,7 +2,7 @@ redshift = Uniform(1e-6,3,name='redshift', latex_label = r'$z$') p0 = Uniform(1, 10, 'p0', latex_label = r'$P_{0]$ [ms]') bp = LogUniform(0.1,10,'bp',latex_label = r'$B_{p}$~[10$^{14}$G]') mass_ns = Uniform(1.1, 2.2, 'mass_ns', latex_label = r'$M_{\mathrm{NS}} [M_{\odot}]$') -theta_pb = Uniform(0, 3.14/2, 'theta_pb', latex_label = r'$\theta_{P-B}$') +theta_pb = Uniform(0, 3.14/2, 'theta_pb', latex_label = r'$\\theta_{P-B}$') mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$') vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') diff --git a/redback/priors/slsn_bolometric.prior b/redback/priors/slsn_bolometric.prior index 2eb9e6b19..0e71a38d3 100644 --- a/redback/priors/slsn_bolometric.prior +++ b/redback/priors/slsn_bolometric.prior @@ -1,7 +1,7 @@ p0 = Uniform(1, 10, 'p0', latex_label = r'$P_{0]$ [ms]') bp = LogUniform(0.1,10,'bp',latex_label = r'$B_{p}$~[10$^{14}$G]') mass_ns = Uniform(1.1, 2.2, 'mass_ns', latex_label = r'$M_{\mathrm{NS}} [M_{\odot}]$') -theta_pb = Uniform(0, 3.14/2, 'theta_pb', latex_label = r'$\theta_{P-B}$') +theta_pb = Uniform(0, 3.14/2, 'theta_pb', latex_label = r'$\\theta_{P-B}$') mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$') vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') diff --git a/redback/priors/smoothpowerlaw.prior b/redback/priors/smoothpowerlaw.prior index ec9f832a8..33559c8b2 100644 --- a/redback/priors/smoothpowerlaw.prior +++ b/redback/priors/smoothpowerlaw.prior @@ -1,8 +1,8 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thv = Sine(name='thv', latex_label=r'$\theta_{\mathrm{observer}}$') +thv = Sine(name='thv', latex_label=r'$\\theta_{\mathrm{observer}}$') loge0 = Uniform(44, 54, 'loge0', latex_label=r'$\log_{10} E_{0}$') -thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') -thw = Uniform(1, 8, 'thw', latex_label=r'$\theta_{\mathrm{truncation}}$') +thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\\theta_{\mathrm{core}}$') +thw = Uniform(1, 8, 'thw', latex_label=r'$\\theta_{\mathrm{truncation}}$') beta = Uniform(0.5, 10, 'beta', latex_label=r'$\\beta$') logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') diff --git a/redback/priors/tophat.prior b/redback/priors/tophat.prior index 73d1f66c2..47bd94afa 100644 --- a/redback/priors/tophat.prior +++ b/redback/priors/tophat.prior @@ -1,7 +1,7 @@ redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') -thv = Sine(name='thv', latex_label=r'$\theta_{\mathrm{observer}}$') +thv = Sine(name='thv', latex_label=r'$\\theta_{\mathrm{observer}}$') loge0 = Uniform(44, 54, 'loge0', latex_label=r'$\log_{10} E_{0}$') -thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\theta_{\mathrm{core}}$') +thc = Uniform(0.01, 0.1, 'thc', latex_label=r'$\\theta_{\mathrm{core}}$') logn0 = Uniform(-5, 2, 'logn0', latex_label=r'$\log_{10} n_{\mathrm{ism}}$') p = Uniform(2, 3, 'p', latex_label=r'$p$') logepse = Uniform(-5, 0, 'logepse', latex_label=r'$\log_{10} \epsilon_{e}$') From bc281ca3459399590fd6b2bdf8527f759c0e1d4a Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 23:40:41 +1100 Subject: [PATCH 10/37] removed some empty lines --- examples/broadband_afterglow_private_data_example.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/examples/broadband_afterglow_private_data_example.py b/examples/broadband_afterglow_private_data_example.py index c9fa8bcc0..c3ceeb1c1 100644 --- a/examples/broadband_afterglow_private_data_example.py +++ b/examples/broadband_afterglow_private_data_example.py @@ -55,5 +55,3 @@ # plot multiband lightcurve. This will plot a panel for every unique frequency result.plot_multiband_lightcurve(random_models=100) - - From ee20b17cedcfe0a2e5822432ed9a233a7029521e Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 23:41:18 +1100 Subject: [PATCH 11/37] Gave all practical transient classes directory structure attributes --- redback/transient/afterglow.py | 4 ++-- redback/transient/kilonova.py | 4 ++-- redback/transient/prompt.py | 1 + redback/transient/supernova.py | 4 ++-- redback/transient/tde.py | 4 ++-- 5 files changed, 9 insertions(+), 8 deletions(-) diff --git a/redback/transient/afterglow.py b/redback/transient/afterglow.py index 8450238f8..2add2cbd0 100644 --- a/redback/transient/afterglow.py +++ b/redback/transient/afterglow.py @@ -109,6 +109,7 @@ def __init__( self._set_photon_index() self._set_t90() self._get_redshift() + self.directory_structure = afterglow_directory_structure(grb=self.name, data_mode=self.data_mode, instrument="") @classmethod def from_swift_grb( @@ -221,7 +222,6 @@ def _save_luminosity_data(self) -> None: """ Saves luminosity data to a csv file. """ - grb_dir, _, _ = afterglow_directory_structure(grb=self._stripped_name, data_mode=self.data_mode) filename = f"{self.name}.csv" data = {"Time in restframe [s]": self.time_rest_frame, "Pos. time err in restframe [s]": self.time_rest_frame_err[0, :], @@ -230,7 +230,7 @@ def _save_luminosity_data(self) -> None: "Pos. luminosity err [10^50 erg s^{-1}]": self.Lum50_err[0, :], "Neg. luminosity err [10^50 erg s^{-1}]": self.Lum50_err[1, :]} df = pd.DataFrame(data=data) - df.to_csv(join(grb_dir, filename), index=False) + df.to_csv(join(self.directory_structure.directory_path, filename), index=False) def _set_data(self) -> None: """ diff --git a/redback/transient/kilonova.py b/redback/transient/kilonova.py index c90c3f161..64892c967 100644 --- a/redback/transient/kilonova.py +++ b/redback/transient/kilonova.py @@ -86,6 +86,6 @@ def __init__( magnitude=magnitude, magnitude_err=magnitude_err, data_mode=data_mode, name=name, bands=bands, system=system, active_bands=active_bands, use_phase_model=use_phase_model, redshift=redshift, photon_index=photon_index, **kwargs) - self.directory_structure = redback.get_data.directory.open_access_directory_structure(transient=name, - transient_type="kilonova") + self.directory_structure = redback.get_data.directory.open_access_directory_structure( + transient=name, transient_type="kilonova") self._set_data() diff --git a/redback/transient/prompt.py b/redback/transient/prompt.py index 06ee1b7c7..26d24883d 100644 --- a/redback/transient/prompt.py +++ b/redback/transient/prompt.py @@ -64,6 +64,7 @@ def __init__( self.channels = channels self.instrument = instrument self._set_data() + self.directory_structure = swift_prompt_directory_structure(grb=self.name, bin_size=self.bin_size) @classmethod def from_batse_grb_name( diff --git a/redback/transient/supernova.py b/redback/transient/supernova.py index 8fe7cc36c..99ab280dc 100644 --- a/redback/transient/supernova.py +++ b/redback/transient/supernova.py @@ -81,6 +81,6 @@ def __init__( magnitude=magnitude, magnitude_err=magnitude_err, data_mode=data_mode, name=name, use_phase_model=use_phase_model, bands=bands, system=system, active_bands=active_bands, redshift=redshift, photon_index=photon_index, **kwargs) - self.directory_structure = redback.get_data.directory.open_access_directory_structure(transient=name, - transient_type="supernova") + self.directory_structure = redback.get_data.directory.open_access_directory_structure( + transient=name, transient_type="supernova") self._set_data() diff --git a/redback/transient/tde.py b/redback/transient/tde.py index 4a54779a0..323b06487 100644 --- a/redback/transient/tde.py +++ b/redback/transient/tde.py @@ -77,6 +77,6 @@ def __init__( magnitude=magnitude, magnitude_err=magnitude_err, data_mode=data_mode, name=name, use_phase_model=use_phase_model, bands=bands, system=system, active_bands=active_bands, redshift=redshift, photon_index=photon_index, **kwargs) - self.directory_structure = redback.get_data.directory.open_access_directory_structure(transient=self.name, - transient_type="tidal_disruption_event") + self.directory_structure = redback.get_data.directory.open_access_directory_structure( + transient=self.name, transient_type="tidal_disruption_event") self._set_data() From a5e1f5dff72d2df189692480eeb590abd8d5a70f Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 23:41:51 +1100 Subject: [PATCH 12/37] Fixed an issue with multiband plotting --- redback/transient/transient.py | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/redback/transient/transient.py b/redback/transient/transient.py index cbcc01c7c..7deda4395 100644 --- a/redback/transient/transient.py +++ b/redback/transient/transient.py @@ -416,7 +416,11 @@ def unique_frequencies(self) -> np.ndarray: ------- np.ndarray: All frequencies that we get from the data, eliminating all duplicates. """ - return self.bands_to_frequencies(self.unique_bands) + try: + if isinstance(self.unique_bands[0], (float, int)): + return self.unique_bands + except (TypeError, IndexError): + return self.bands_to_frequencies(self.unique_bands) @property def list_of_band_indices(self) -> list: @@ -587,17 +591,18 @@ def plot_multiband_lightcurve( times = self._get_times(axes) times_mesh, frequency_mesh = np.meshgrid(times, self.unique_frequencies) - model_kwargs['frequency'] = frequency_mesh + new_model_kwargs = model_kwargs.copy() + new_model_kwargs['frequency'] = frequency_mesh posterior.sort_values(by='log_likelihood') max_like_params = posterior.iloc[-1] - ys = model(times_mesh, **max_like_params, **model_kwargs) + ys = model(times_mesh, **max_like_params, **new_model_kwargs) for i in range(len(self.unique_frequencies)): axes[i].plot(times_mesh[i], ys[i], color='blue', alpha=0.65, lw=2) - params = posterior.iloc[np.random.randint(len(posterior))] - ys = model(times_mesh, **params, **model_kwargs) for _ in range(random_models): - axes.plot(times, ys[i], color='red', alpha=0.05, lw=2, zorder=-1) + params = posterior.iloc[np.random.randint(len(posterior))] + ys = model(times_mesh, **params, **new_model_kwargs) + axes[i].plot(times, ys[i], color='red', alpha=0.05, lw=2, zorder=-1) if plot_save: plt.savefig(join(outdir, filename), dpi=300, bbox_inches="tight") if plot_show: @@ -616,7 +621,12 @@ def _get_times(self, axes: matplotlib.axes.Axes) -> np.ndarray: np.ndarray: Linearly or logarithmically scaled time values depending on the y scale used in the plot. """ - if axes.get_yscale == 'linear': + if isinstance(axes, np.ndarray): + ax = axes[0] + else: + ax = axes + + if ax.get_yscale() == 'linear': times = np.linspace(self.x[0], self.x[-1], 200) else: times = np.exp(np.linspace(np.log(self.x[0]), np.log(self.x[-1]), 200)) From 1699881098967417899785d88791708c04173ffa Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 3 Mar 2022 23:42:10 +1100 Subject: [PATCH 13/37] Fixed some more minor issues --- redback/plotting.py | 25 +++++++++++-------------- redback/sampler.py | 4 +--- 2 files changed, 12 insertions(+), 17 deletions(-) diff --git a/redback/plotting.py b/redback/plotting.py index 70b779067..21e37cafb 100644 --- a/redback/plotting.py +++ b/redback/plotting.py @@ -1,10 +1,10 @@ +from os.path import join + import matplotlib -import numpy as np import matplotlib.pyplot as plt -from os.path import join +import numpy as np import redback -from redback.get_data.directory import afterglow_directory_structure class MultiBandPlotter(object): @@ -79,13 +79,14 @@ def plot_multiband( axes = axes.ravel() - for i, (indices, band) in enumerate(zip(self.transient.list_of_band_indices, self.transient.unique_bands)): + i = 0 + for indices, band in zip(self.transient.list_of_band_indices, self.transient.unique_bands): if band not in filters: continue x_err = self.transient.x_err[indices] if self.transient.x_err is not None else self.transient.x_err - color = colors[filters.index(band)] + color = colors[list(filters).index(band)] if isinstance(band, str): freq = self.transient.bands_to_frequencies([band]) @@ -107,13 +108,13 @@ def plot_multiband( axes[i].set_yscale("log") axes[i].legend(ncol=2) axes[i].tick_params(axis='both', which='major', pad=8) + i += 1 figure.supxlabel(xlabel, fontsize=fontsize) figure.supylabel(ylabel, fontsize=fontsize) filename = f"{self.transient.name}_{plot_label}.png" plt.subplots_adjust(wspace=wspace, hspace=hspace) - directory_structure = afterglow_directory_structure(grb=self.transient.name, data_mode=self.transient.data_mode) - plt.savefig(join(directory_structure.directory_path, filename), bbox_inches="tight") + plt.savefig(join(self.transient.directory_structure.directory_path, filename), bbox_inches="tight") return axes @@ -166,9 +167,8 @@ def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k', **kwar if axes is None: plt.tight_layout() - directory_structure = afterglow_directory_structure(grb=self.transient.name, data_mode=self.transient.data_mode) filename = f"{self.transient.name}_lc.png" - plt.savefig(join(directory_structure.directory_path, filename)) + plt.savefig(join(self.transient.directory_structure.directory_path, filename)) if axes is None: plt.clf() return ax @@ -218,7 +218,7 @@ def plot_data( for indices, band in zip(self.transient.list_of_band_indices, self.transient.unique_bands): x_err = self.transient.x_err[indices] if self.transient.x_err is not None else self.transient.x_err if band in filters: - color = colors[filters.index(band)] + color = colors[list(filters).index(band)] label = band elif plot_others: color = "black" @@ -243,10 +243,7 @@ def plot_data( if plot_save: filename = f"{self.transient.name}_{plot_label}.png" - structure = afterglow_directory_structure( - grb=self.transient.name, data_mode=self.transient.data_mode, instrument='') - - plt.savefig(join(structure.directory_path, filename), bbox_inches='tight') + plt.savefig(join(self.transient.directory_structure.directory_path, filename), bbox_inches='tight') plt.clf() return axes diff --git a/redback/sampler.py b/redback/sampler.py index 3077226b9..f1adb2319 100644 --- a/redback/sampler.py +++ b/redback/sampler.py @@ -92,9 +92,7 @@ def _fit_grb(transient, model, outdir=None, label=None, sampler='dynesty', nlive prior['alpha_1'] = bilby.prior.Gaussian(mu=-(transient.photon_index + 1), sigma=0.1, latex_label=r'$\alpha_{1}$') if outdir is None: - outdir, _, _ = redback.get_data.directory.afterglow_directory_structure( - grb=transient.name, data_mode=transient.data_mode, instrument='') - outdir = f"{outdir}/{model.__name__}" + outdir = f"{transient.directory_structure.directory_path}/{model.__name__}" Path(outdir).mkdir(parents=True, exist_ok=True) label = kwargs.get("label", transient.name) From e2f1637ac363fd339b470490cfcc3ed5a0a9d71e Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 10 Mar 2022 19:23:30 +1100 Subject: [PATCH 14/37] some reformatting --- examples/kilonova_example.py | 2 +- redback/transient/transient.py | 9 +++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/examples/kilonova_example.py b/examples/kilonova_example.py index 9dfed3ce8..b0166b7c3 100644 --- a/examples/kilonova_example.py +++ b/examples/kilonova_example.py @@ -26,7 +26,7 @@ model_kwargs = dict(frequency=kilonova.filtered_frequencies, output_format='flux_density') result = redback.fit_model(name=kne, transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, data_mode='flux_density', sample='rslice', nlive=200, resume=False) + prior=priors, data_mode='flux_density', sample='rslice', nlive=200, resume=True) result.plot_corner() # returns a Kilonova result object result.plot_lightcurve(random_models=1000) diff --git a/redback/transient/transient.py b/redback/transient/transient.py index 7deda4395..ce69b77fd 100644 --- a/redback/transient/transient.py +++ b/redback/transient/transient.py @@ -420,7 +420,8 @@ def unique_frequencies(self) -> np.ndarray: if isinstance(self.unique_bands[0], (float, int)): return self.unique_bands except (TypeError, IndexError): - return self.bands_to_frequencies(self.unique_bands) + pass + return self.bands_to_frequencies(self.unique_bands) @property def list_of_band_indices(self) -> list: @@ -535,9 +536,13 @@ def plot_lightcurve( plt.semilogy() times = self._get_times(axes) + times_mesh, frequency_mesh = np.meshgrid(times, self.unique_frequencies) + new_model_kwargs = model_kwargs.copy() + new_model_kwargs['frequency'] = frequency_mesh + posterior.sort_values(by='log_likelihood') max_like_params = posterior.iloc[-1] - ys = model(times, **max_like_params, **model_kwargs) + ys = model(times, **max_like_params, **new_model_kwargs) axes.plot(times, ys, color='blue', alpha=0.65, lw=2) for _ in range(random_models): From 0447cd5cda0958afcef52203c557ef93eaf4e0f2 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 10 Mar 2022 19:33:29 +1100 Subject: [PATCH 15/37] Renamed bands_to_frequencies --> bands_to_frequency --- examples/SN2011kl_sample_in_t0_example.py | 2 +- examples/magnetar_boosted_example.py | 2 +- examples/supernova_example.py | 2 +- examples/tde_example.py | 2 +- redback/plotting.py | 2 +- redback/transient/kilonova.py | 4 ++-- redback/transient/supernova.py | 4 ++-- redback/transient/tde.py | 4 ++-- redback/transient/transient.py | 16 ++++++++-------- test/transient_test.py | 12 ++++++------ 10 files changed, 25 insertions(+), 25 deletions(-) diff --git a/examples/SN2011kl_sample_in_t0_example.py b/examples/SN2011kl_sample_in_t0_example.py index 1b96a7dce..c1118fe51 100644 --- a/examples/SN2011kl_sample_in_t0_example.py +++ b/examples/SN2011kl_sample_in_t0_example.py @@ -42,7 +42,7 @@ priors['av'] = Uniform(0.1, 1, name='av', latex_label=r'$a_{v}$') -model_kwargs = dict(frequencies=redback.utils.bands_to_frequencies(bands), output_format='flux_density') +model_kwargs = dict(frequencies=redback.utils.bands_to_frequency(bands), output_format='flux_density') # returns a supernova result object result = redback.fit_model(name=sne, transient=supernova, model=model, sampler=sampler, model_kwargs=model_kwargs, diff --git a/examples/magnetar_boosted_example.py b/examples/magnetar_boosted_example.py index 2187c7c1a..0b6eaa43c 100644 --- a/examples/magnetar_boosted_example.py +++ b/examples/magnetar_boosted_example.py @@ -27,7 +27,7 @@ priors['tau_sd'] = 1e3 priors['thermalisation_efficiency'] = 0.3 -model_kwargs = dict(frequencies=redback.utils.bands_to_frequencies(bands), output_format='flux_density') +model_kwargs = dict(frequencies=redback.utils.bands_to_frequency(bands), output_format='flux_density') result = redback.fit_model(name=kne, transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, path=path, prior=priors, data_mode='flux_density', sample='rslice', nlive=200) diff --git a/examples/supernova_example.py b/examples/supernova_example.py index 6b5900639..c9cb342ca 100644 --- a/examples/supernova_example.py +++ b/examples/supernova_example.py @@ -19,7 +19,7 @@ # use default priors priors = redback.priors.get_priors(model=model) priors['redshift'] = 0.677 -model_kwargs = dict(frequencies=redback.utils.bands_to_frequencies(bands), output_format='flux_density') +model_kwargs = dict(frequencies=redback.utils.bands_to_frequency(bands), output_format='flux_density') # returns a supernova result object result = redback.fit_model(name=sne, transient=supernova, model=model, sampler=sampler, model_kwargs=model_kwargs, diff --git a/examples/tde_example.py b/examples/tde_example.py index 511296d96..8e4375ffa 100644 --- a/examples/tde_example.py +++ b/examples/tde_example.py @@ -25,7 +25,7 @@ # We can do this through the metadata that was downloaded alongside the data, or if you just know it. priors['redshift'] = tidal_disruption_event.redshift -model_kwargs = dict(frequencies=redback.utils.bands_to_frequencies(bands), output_format='flux_density') +model_kwargs = dict(frequencies=redback.utils.bands_to_frequency(bands), output_format='flux_density') # returns a tde result object result = redback.fit_model(name=tde, transient=tidal_disruption_event, model=model, sampler=sampler, model_kwargs=model_kwargs, diff --git a/redback/plotting.py b/redback/plotting.py index 21e37cafb..458396277 100644 --- a/redback/plotting.py +++ b/redback/plotting.py @@ -89,7 +89,7 @@ def plot_multiband( color = colors[list(filters).index(band)] if isinstance(band, str): - freq = self.transient.bands_to_frequencies([band]) + freq = self.transient.bands_to_frequency([band]) if 1e10 < freq < 1e15: label = band else: diff --git a/redback/transient/kilonova.py b/redback/transient/kilonova.py index 64892c967..660b534ef 100644 --- a/redback/transient/kilonova.py +++ b/redback/transient/kilonova.py @@ -77,8 +77,8 @@ def __init__( Whether we are using a phase model. kwargs: dict, optional Additional callables: - bands_to_frequencies: Conversion function to convert a list of bands to frequencies. Use - redback.utils.bands_to_frequencies if not given. + bands_to_frequency: Conversion function to convert a list of bands to frequencies. Use + redback.utils.bands_to_frequency if not given. """ super().__init__(time=time, time_err=time_err, time_rest_frame=time_rest_frame, time_mjd=time_mjd, time_mjd_err=time_mjd_err, time_rest_frame_err=time_rest_frame_err, Lum50=Lum50, diff --git a/redback/transient/supernova.py b/redback/transient/supernova.py index 99ab280dc..e9fe80ee6 100644 --- a/redback/transient/supernova.py +++ b/redback/transient/supernova.py @@ -71,8 +71,8 @@ def __init__( Whether we are using a phase model. kwargs: dict, optional Additional callables: - bands_to_frequencies: Conversion function to convert a list of bands to frequencies. Use - redback.utils.bands_to_frequencies if not given. + bands_to_frequency: Conversion function to convert a list of bands to frequencies. Use + redback.utils.bands_to_frequency if not given. """ super().__init__(time=time, time_err=time_err, time_rest_frame=time_rest_frame, time_mjd=time_mjd, diff --git a/redback/transient/tde.py b/redback/transient/tde.py index 323b06487..d2debcfbd 100644 --- a/redback/transient/tde.py +++ b/redback/transient/tde.py @@ -68,8 +68,8 @@ def __init__( Whether we are using a phase model. kwargs: dict, optional Additional callables: - bands_to_frequencies: Conversion function to convert a list of bands to frequencies. Use - redback.utils.bands_to_frequencies if not given. + bands_to_frequency: Conversion function to convert a list of bands to frequencies. Use + redback.utils.bands_to_frequency if not given. """ super().__init__(time=time, time_err=time_err, time_rest_frame=time_rest_frame, time_mjd=time_mjd, time_mjd_err=time_mjd_err, time_rest_frame_err=time_rest_frame_err, Lum50=Lum50, diff --git a/redback/transient/transient.py b/redback/transient/transient.py index ce69b77fd..9b56b8b30 100644 --- a/redback/transient/transient.py +++ b/redback/transient/transient.py @@ -102,13 +102,13 @@ def __init__( List or array of active bands to be used in the analysis. Use all available bands if 'all' is given. kwargs: dict, optional Additional callables: - bands_to_frequencies: Conversion function to convert a list of bands to frequencies. Use - redback.utils.bands_to_frequencies if not given. - bin_ttes: Binning function for time-tagged event data. Use redback.utils.bands_to_frequencies if not given. + bands_to_frequency: Conversion function to convert a list of bands to frequencies. Use + redback.utils.bands_to_frequency if not given. + bin_ttes: Binning function for time-tagged event data. Use redback.utils.bands_to_frequency if not given. """ self.bin_size = bin_size self.bin_ttes = kwargs.get("bin_ttes", redback.utils.bin_ttes) - self.bands_to_frequencies = kwargs.get("bands_to_frequencies", redback.utils.bands_to_frequencies) + self.bands_to_frequency = kwargs.get("bands_to_frequency", redback.utils.bands_to_frequency) if data_mode == 'ttes': time, counts = self.bin_ttes(ttes, self.bin_size) @@ -311,7 +311,7 @@ def set_bands_and_frequency( self._bands = self.frequency elif bands is not None and frequency is None: self._bands = bands - self._frequency = self.bands_to_frequencies(self.bands) + self._frequency = self.bands_to_frequency(self.bands) @property def frequency(self) -> np.ndarray: @@ -421,7 +421,7 @@ def unique_frequencies(self) -> np.ndarray: return self.unique_bands except (TypeError, IndexError): pass - return self.bands_to_frequencies(self.unique_bands) + return self.bands_to_frequency(self.unique_bands) @property def list_of_band_indices(self) -> list: @@ -737,8 +737,8 @@ def __init__( Whether we are using a phase model. kwargs: dict, optional Additional callables: - bands_to_frequencies: Conversion function to convert a list of bands to frequencies. Use - redback.utils.bands_to_frequencies if not given. + bands_to_frequency: Conversion function to convert a list of bands to frequencies. Use + redback.utils.bands_to_frequency if not given. """ super().__init__(time=time, time_err=time_err, time_rest_frame=time_rest_frame, time_mjd=time_mjd, time_mjd_err=time_mjd_err, frequency=frequency, diff --git a/test/transient_test.py b/test/transient_test.py index ca5303a8e..12cb97b5e 100644 --- a/test/transient_test.py +++ b/test/transient_test.py @@ -303,14 +303,14 @@ def test_set_active_bands_all(self): def test_set_frequencies_from_bands(self): expected = [1, 2, 2] - bands_to_frequencies = MagicMock(return_value=expected) + bands_to_frequency = MagicMock(return_value=expected) self.transient = redback.transient.transient.OpticalTransient( time=self.time, time_err=self.time_err, flux_density=self.y, flux_density_err=self.y_err, redshift=self.redshift, data_mode=self.data_mode, name=self.name, photon_index=self.photon_index, use_phase_model=self.use_phase_model, bands=self.bands, - active_bands=self.active_bands, bands_to_frequencies=bands_to_frequencies) + active_bands=self.active_bands, bands_to_frequency=bands_to_frequency) self.assertTrue(np.array_equal(expected, self.transient.frequency)) - bands_to_frequencies.assert_called_once() + bands_to_frequency.assert_called_once() def test_set_frequencies_default(self): frequency = np.array([1, 2, 2]) @@ -471,14 +471,14 @@ def test_set_active_bands_all(self): def test_set_frequencies_from_bands(self): expected = [1, 2, 2] - bands_to_frequencies = MagicMock(return_value=expected) + bands_to_frequency = MagicMock(return_value=expected) self.sgrb = redback.transient.afterglow.SGRB( time=self.time, time_err=self.time_err, flux_density=self.y, flux_density_err=self.y_err, data_mode=self.data_mode, name=self.name, use_phase_model=self.use_phase_model, bands=self.bands, - active_bands=self.active_bands, bands_to_frequencies=bands_to_frequencies) + active_bands=self.active_bands, bands_to_frequency=bands_to_frequency) self.assertTrue(np.array_equal(expected, self.sgrb.frequency)) - bands_to_frequencies.assert_called_once() + bands_to_frequency.assert_called_once() def test_set_frequencies_default(self): frequency = np.array([1, 2, 2]) From 6bc9c87dbd1cbbdc1db1baed149e8bb032d02433 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 10 Mar 2022 20:33:07 +1100 Subject: [PATCH 16/37] Sorted out integrated flux plotting --- examples/magnetar_example.py | 3 +- redback/plotting.py | 129 ++++++++++++++++++++++++++++++--- redback/transient/afterglow.py | 50 ++++++++++++- 3 files changed, 171 insertions(+), 11 deletions(-) diff --git a/examples/magnetar_example.py b/examples/magnetar_example.py index ce6770be5..acf7a623d 100644 --- a/examples/magnetar_example.py +++ b/examples/magnetar_example.py @@ -16,6 +16,7 @@ # uses an analytical k-correction expression to create luminosity data if not already there. # Can also use a numerical k-correction through CIAO afterglow.analytical_flux_to_luminosity() +afterglow.plot_data() # use default priors priors = redback.priors.get_priors(model=model, data_mode='luminosity') @@ -34,6 +35,6 @@ # Call redback.fit_model to run the sampler and obtain GRB result object result = redback.fit_model(name=GRB, model=model, sampler='dynesty', nlive=200, transient=afterglow, - prior=priors, data_mode='luminosity', sample='rslice') + prior=priors, sample='rslice', resume=True) result.plot_lightcurve(random_models=1000) diff --git a/redback/plotting.py b/redback/plotting.py index 458396277..62526c5bc 100644 --- a/redback/plotting.py +++ b/redback/plotting.py @@ -3,6 +3,7 @@ import matplotlib import matplotlib.pyplot as plt import numpy as np +import pandas as pd import redback @@ -123,7 +124,25 @@ class IntegratedFluxPlotter(object): def __init__(self, transient): self.transient = transient - def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k', **kwargs) -> matplotlib.axes.Axes: + @property + def xlim_low(self): + return 0.5 * self.transient.x[0] + + @property + def xlim_high(self): + if self.x_err is None: + return 2 * self.transient.x[-1] + return 2 * (self.transient.x[-1] + self.x_err[1][-1]) + + @property + def x_err(self): + if self.transient.x_err is not None: + return [np.abs(self.transient.x_err[1, :]), self.transient.x_err[0, :]] + else: + return None + + def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k', plot_save: bool = True, + plot_show: bool = True, **kwargs) -> matplotlib.axes.Axes: """ Plots the Afterglow lightcurve and returns Axes. @@ -141,20 +160,16 @@ def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k', **kwar matplotlib.axes.Axes: The axes with the plot. """ - if self.transient.x_err is not None: - x_err = [np.abs(self.transient.x_err[1, :]), self.transient.x_err[0, :]] - else: - x_err = None y_err = [np.abs(self.transient.y_err[1, :]), self.transient.y_err[0, :]] ax = axes or plt.gca() - ax.errorbar(self.transient.x, self.transient.y, xerr=x_err, yerr=y_err, + ax.errorbar(self.transient.x, self.transient.y, xerr=self.x_err, yerr=y_err, fmt='x', c=colour, ms=1, elinewidth=2, capsize=0.) ax.set_xscale('log') ax.set_yscale('log') - ax.set_xlim(0.5 * self.transient.x[0], 2 * (self.transient.x[-1] + x_err[1][-1])) + ax.set_xlim(self.xlim_low, self.xlim_high) ax.set_ylim(0.5 * min(self.transient.y), 2. * np.max(self.transient.y)) ax.annotate(self.transient.name, xy=(0.95, 0.9), xycoords='axes fraction', @@ -167,12 +182,97 @@ def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k', **kwar if axes is None: plt.tight_layout() - filename = f"{self.transient.name}_lc.png" - plt.savefig(join(self.transient.directory_structure.directory_path, filename)) + filename = f"{self.transient.name}_data.png" + if plot_save: + plt.tight_layout() + plt.savefig(join(self.transient.directory_structure.directory_path, filename)) + if plot_show: + plt.tight_layout() + plt.show() if axes is None: plt.clf() return ax + def plot_lightcurve( + self, model: callable, filename: str = None, axes: matplotlib.axes.Axes = None, plot_save: bool = True, + plot_show: bool = True, random_models: int = 100, posterior: pd.DataFrame = None, outdir: str = '.', + model_kwargs: dict = None, **kwargs: object) -> None: + """ + + Parameters + ---------- + model: callable + The model used to plot the lightcurve. + filename: str, optional + The output filename. Otherwise, use default which starts with the name + attribute and ends with *lightcurve.png. + axes: matplotlib.axes.Axes, optional + Axes to plot in if given. + plot_save: bool, optional + Whether to save the plot. + plot_show: bool, optional + Whether to show the plot. + random_models: int, optional + Number of random posterior samples plotted faintly. Default is 100. + posterior: pd.DataFrame, optional + Posterior distribution to which to draw samples from. Is optional but must be given. + outdir: str, optional + Out directory in which to save the plot. Default is the current working directory. + model_kwargs: dict + Additional keyword arguments to be passed into the model. + kwargs: dict + No current function. + """ + if filename is None: + filename = f"{self.transient.name}_lightcurve.png" + if model_kwargs is None: + model_kwargs = dict() + axes = axes or plt.gca() + axes = self.plot_data(axes=axes, plot_save=False, plot_show=False) + axes.set_yscale('log') + plt.semilogy() + times = self._get_times(axes) + + posterior.sort_values(by='log_likelihood') + max_like_params = posterior.iloc[-1] + ys = model(times, **max_like_params, **model_kwargs) + axes.plot(times, ys, color='blue', alpha=0.65, lw=2) + + for _ in range(random_models): + params = posterior.iloc[np.random.randint(len(posterior))] + ys = model(times, **params, **model_kwargs) + axes.plot(times, ys, color='red', alpha=0.05, lw=2, zorder=-1, **kwargs) + + if plot_save: + plt.savefig(join(outdir, filename), dpi=300, bbox_inches="tight") + if plot_show: + plt.tight_layout() + plt.show() + plt.clf() + + def _get_times(self, axes: matplotlib.axes.Axes) -> np.ndarray: + """ + + Parameters + ---------- + axes: matplotlib.axes.Axes + The axes used in the plotting procedure. + Returns + ------- + np.ndarray: Linearly or logarithmically scaled time values depending on the y scale used in the plot. + + """ + if isinstance(axes, np.ndarray): + ax = axes[0] + else: + ax = axes + + if ax.get_yscale() == 'linear': + times = np.linspace(self.xlim_low, self.xlim_high, 200) + else: + times = np.exp(np.linspace(np.log(self.xlim_low), np.log(self.xlim_high), 200)) + return times + class LuminosityPlotter(IntegratedFluxPlotter): pass @@ -251,6 +351,17 @@ def _set_y_axis(self, ax): ax.set_ylim(0.8 * min(self.transient.y), 1.2 * np.max(self.transient.y)) ax.invert_yaxis() + def plot_lightcurve(self): + pass + + def plot_multiband(self): + pass + + def plot_multiband_lightcurve(self): + pass + + + class FluxDensityPlotter(MagnitudePlotter): diff --git a/redback/transient/afterglow.py b/redback/transient/afterglow.py index 2add2cbd0..81a3cf421 100644 --- a/redback/transient/afterglow.py +++ b/redback/transient/afterglow.py @@ -148,6 +148,35 @@ def from_swift_grb( def _stripped_name(self) -> str: return self.name.lstrip('GRB') + @property + def data_mode(self) -> str: + """ + + Returns + ------- + str: The currently active data mode (one in `Transient.DATA_MODES`) + """ + return self._data_mode + + @data_mode.setter + def data_mode(self, data_mode: str) -> None: + """ + + Parameters + ------- + data_mode: str + One of the data modes in `Transient.DATA_MODES` + """ + if data_mode in self.DATA_MODES or data_mode is None: + self._data_mode = data_mode + try: + self.directory_structure = afterglow_directory_structure( + grb=self.name, data_mode=self.data_mode, instrument="") + except AttributeError: + pass + else: + raise ValueError("Unknown data mode.") + def load_and_truncate_data( self, truncate: bool = True, truncate_method: str = 'prompt_time_error', data_mode: str = 'flux') -> None: """ @@ -367,14 +396,33 @@ def _convert_flux_to_luminosity( redshift = self._get_redshift_for_luminosity_calculation() if redshift is None: return + self.data_mode = "luminosity" converter = self.FluxToLuminosityConverter( redshift=redshift, photon_index=self.photon_index, time=self.time, time_err=self.time_err, flux=self.flux, flux_err=self.flux_err, counts_to_flux_absorbed=counts_to_flux_absorbed, counts_to_flux_unabsorbed=counts_to_flux_unabsorbed, conversion_method=conversion_method) - self.data_mode = "luminosity" self.x, self.x_err, self.y, self.y_err = converter.convert_flux_to_luminosity() self._save_luminosity_data() + def plot_lightcurve( + self, model: callable, filename: str = None, axes: matplotlib.axes.Axes = None, plot_save: bool = True, + plot_show: bool = True, random_models: int = 100, posterior: pd.DataFrame = None, outdir: str = '.', + model_kwargs: dict = None, **kwargs: object) -> None: + if self.flux_data: + plotter = IntegratedFluxPlotter(transient=self) + elif self.luminosity_data: + plotter = LuminosityPlotter(transient=self) + elif self.flux_density_data: + plotter = FluxDensityPlotter(transient=self) + elif self.magnitude_data: + plotter = MagnitudePlotter(transient=self) + else: + return axes + return plotter.plot_lightcurve( + model=model, filename=filename, axes=axes, plot_save=plot_save, + plot_show=plot_show, random_models=random_models, posterior=posterior, + outdir=outdir, model_kwargs=model_kwargs, **kwargs) + def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k', **kwargs: dict) -> matplotlib.axes.Axes: """ Plots the Afterglow lightcurve and returns Axes. From b0f4db9e09e4313ee7dcb0275f2346bf59e23436 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 10 Mar 2022 21:00:07 +1100 Subject: [PATCH 17/37] Rationalised plotting --- redback/plotting.py | 309 +++++++++++++++++++-------------- redback/transient/afterglow.py | 27 ++- redback/transient/transient.py | 268 +++++++++------------------- 3 files changed, 289 insertions(+), 315 deletions(-) diff --git a/redback/plotting.py b/redback/plotting.py index 62526c5bc..08dbc6dd0 100644 --- a/redback/plotting.py +++ b/redback/plotting.py @@ -8,118 +8,33 @@ import redback -class MultiBandPlotter(object): +class Plotter(object): - def __init__(self, transient): - self.transient = transient - - def plot_multiband( - self, figure: matplotlib.figure.Figure = None, axes: matplotlib.axes.Axes = None, ncols: int = 2, - nrows: int = None, figsize: tuple = None, filters: list = None, **plot_kwargs: dict) -> \ - matplotlib.axes.Axes: + def _get_times(self, axes: matplotlib.axes.Axes) -> np.ndarray: """ Parameters ---------- - figure: matplotlib.figure.Figure, optional - Figure can be given if defaults are not satisfying - axes: matplotlib.axes.Axes, optional - Axes can be given if defaults are not satisfying - ncols: int, optional - Number of columns to use on the plot. Default is 2. - nrows: int, optional - Number of rows to use on the plot. If None are given this will - be inferred from ncols and the number of filters. - figsize: tuple, optional - Size of the figure. A default based on ncols and nrows will be used if None is given. - filters: list, optional - Which bands to plot. Will use default filters if None is given. - plot_kwargs: - Additional optional plotting kwargs: - wspace: Extra argument for matplotlib.pyplot.subplots_adjust - hspace: Extra argument for matplotlib.pyplot.subplots_adjust - fontsize: Label fontsize - errorbar_fmt: Errorbar format ('fmt' argument in matplotlib.pyplot.errorbar) - colors: colors to be used for the bands - xlabel: Plot xlabel - ylabel: Plot ylabel - plot_label: Addional filename label appended to the default name - + axes: matplotlib.axes.Axes + The axes used in the plotting procedure. Returns ------- + np.ndarray: Linearly or logarithmically scaled time values depending on the y scale used in the plot. """ - if self.transient.luminosity_data or self.transient.flux_data: - redback.utils.logger.warning(f"Can't plot multiband for {self.transient.data_mode} data.") - return - - if filters is None: - filters = self.transient.active_bands - elif str(filters) == 'default': - filters = self.transient.default_filters - - wspace = plot_kwargs.get("wspace", 0.15) - hspace = plot_kwargs.get("hspace", 0.04) - fontsize = plot_kwargs.get("fontsize", 30) - errorbar_fmt = plot_kwargs.get("errorbar_fmt", "x") - colors = plot_kwargs.get("colors", self.transient.get_colors(filters)) - xlabel = plot_kwargs.get("xlabel", self.transient.xlabel) - ylabel = plot_kwargs.get("ylabel", self.transient.ylabel) - plot_label = plot_kwargs.get("plot_label", "multiband_data") - - if figure is None or axes is None: - if nrows is None: - nrows = int(np.ceil(len(filters) / 2)) - npanels = ncols * nrows - if npanels < len(filters): - raise ValueError(f"Insufficient number of panels. {npanels} panels were given " - f"but {len(filters)} panels are needed.") - if figsize is None: - figsize = (4 * ncols, 4 * nrows) - figure, axes = plt.subplots(ncols=ncols, nrows=nrows, sharex='all', figsize=figsize) - - axes = axes.ravel() - - i = 0 - for indices, band in zip(self.transient.list_of_band_indices, self.transient.unique_bands): - if band not in filters: - continue - - x_err = self.transient.x_err[indices] if self.transient.x_err is not None else self.transient.x_err - - color = colors[list(filters).index(band)] - - if isinstance(band, str): - freq = self.transient.bands_to_frequency([band]) - if 1e10 < freq < 1e15: - label = band - else: - label = freq - else: - label = f"{band:.2e}" - axes[i].errorbar(self.transient.x[indices], self.transient.y[indices], xerr=x_err, yerr=self.transient.y_err[indices], - fmt=errorbar_fmt, ms=1, color=color, elinewidth=2, capsize=0., label=label) - - axes[i].set_xlim(0.5 * self.transient.x[indices][0], 1.2 * self.transient.x[indices][-1]) - if self.transient.magnitude_data: - axes[i].set_ylim(0.8 * min(self.transient.y[indices]), 1.2 * np.max(self.transient.y[indices])) - axes[i].invert_yaxis() - else: - axes[i].set_ylim(0.5 * min(self.transient.y[indices]), 2. * np.max(self.transient.y[indices])) - axes[i].set_yscale("log") - axes[i].legend(ncol=2) - axes[i].tick_params(axis='both', which='major', pad=8) - i += 1 + if isinstance(axes, np.ndarray): + ax = axes[0] + else: + ax = axes - figure.supxlabel(xlabel, fontsize=fontsize) - figure.supylabel(ylabel, fontsize=fontsize) - filename = f"{self.transient.name}_{plot_label}.png" - plt.subplots_adjust(wspace=wspace, hspace=hspace) - plt.savefig(join(self.transient.directory_structure.directory_path, filename), bbox_inches="tight") - return axes + if ax.get_yscale() == 'linear': + times = np.linspace(self.xlim_low, self.xlim_high, 200) + else: + times = np.exp(np.linspace(np.log(self.xlim_low), np.log(self.xlim_high), 200)) + return times -class IntegratedFluxPlotter(object): +class IntegratedFluxPlotter(Plotter): def __init__(self, transient): self.transient = transient @@ -250,35 +165,12 @@ def plot_lightcurve( plt.show() plt.clf() - def _get_times(self, axes: matplotlib.axes.Axes) -> np.ndarray: - """ - - Parameters - ---------- - axes: matplotlib.axes.Axes - The axes used in the plotting procedure. - Returns - ------- - np.ndarray: Linearly or logarithmically scaled time values depending on the y scale used in the plot. - - """ - if isinstance(axes, np.ndarray): - ax = axes[0] - else: - ax = axes - - if ax.get_yscale() == 'linear': - times = np.linspace(self.xlim_low, self.xlim_high, 200) - else: - times = np.exp(np.linspace(np.log(self.xlim_low), np.log(self.xlim_high), 200)) - return times - class LuminosityPlotter(IntegratedFluxPlotter): pass -class MagnitudePlotter(object): +class MagnitudePlotter(Plotter): def __init__(self, transient): self.transient = transient @@ -354,12 +246,173 @@ def _set_y_axis(self, ax): def plot_lightcurve(self): pass - def plot_multiband(self): - pass + def plot_multiband( + self, figure: matplotlib.figure.Figure = None, axes: matplotlib.axes.Axes = None, ncols: int = 2, + nrows: int = None, figsize: tuple = None, filters: list = None, **plot_kwargs: dict) -> \ + matplotlib.axes.Axes: + """ - def plot_multiband_lightcurve(self): - pass + Parameters + ---------- + figure: matplotlib.figure.Figure, optional + Figure can be given if defaults are not satisfying + axes: matplotlib.axes.Axes, optional + Axes can be given if defaults are not satisfying + ncols: int, optional + Number of columns to use on the plot. Default is 2. + nrows: int, optional + Number of rows to use on the plot. If None are given this will + be inferred from ncols and the number of filters. + figsize: tuple, optional + Size of the figure. A default based on ncols and nrows will be used if None is given. + filters: list, optional + Which bands to plot. Will use default filters if None is given. + plot_kwargs: + Additional optional plotting kwargs: + wspace: Extra argument for matplotlib.pyplot.subplots_adjust + hspace: Extra argument for matplotlib.pyplot.subplots_adjust + fontsize: Label fontsize + errorbar_fmt: Errorbar format ('fmt' argument in matplotlib.pyplot.errorbar) + colors: colors to be used for the bands + xlabel: Plot xlabel + ylabel: Plot ylabel + plot_label: Addional filename label appended to the default name + + Returns + ------- + """ + if self.transient.luminosity_data or self.transient.flux_data: + redback.utils.logger.warning(f"Can't plot multiband for {self.transient.data_mode} data.") + return + + if filters is None: + filters = self.transient.active_bands + elif str(filters) == 'default': + filters = self.transient.default_filters + + wspace = plot_kwargs.get("wspace", 0.15) + hspace = plot_kwargs.get("hspace", 0.04) + fontsize = plot_kwargs.get("fontsize", 30) + errorbar_fmt = plot_kwargs.get("errorbar_fmt", "x") + colors = plot_kwargs.get("colors", self.transient.get_colors(filters)) + xlabel = plot_kwargs.get("xlabel", self.transient.xlabel) + ylabel = plot_kwargs.get("ylabel", self.transient.ylabel) + plot_label = plot_kwargs.get("plot_label", "multiband_data") + + if figure is None or axes is None: + if nrows is None: + nrows = int(np.ceil(len(filters) / 2)) + npanels = ncols * nrows + if npanels < len(filters): + raise ValueError(f"Insufficient number of panels. {npanels} panels were given " + f"but {len(filters)} panels are needed.") + if figsize is None: + figsize = (6 * nrows, 4 * ncols) + figure, axes = plt.subplots(ncols=ncols, nrows=nrows, sharex='all', figsize=figsize) + + axes = axes.ravel() + + i = 0 + for indices, band in zip(self.transient.list_of_band_indices, self.transient.unique_bands): + if band not in filters: + continue + + x_err = self.transient.x_err[indices] if self.transient.x_err is not None else self.transient.x_err + + color = colors[list(filters).index(band)] + + if isinstance(band, str): + freq = self.transient.bands_to_frequency([band]) + if 1e10 < freq < 1e15: + label = band + else: + label = freq + else: + label = f"{band:.2e}" + axes[i].errorbar(self.transient.x[indices], self.transient.y[indices], xerr=x_err, yerr=self.transient.y_err[indices], + fmt=errorbar_fmt, ms=1, color=color, elinewidth=2, capsize=0., label=label) + + axes[i].set_xlim(0.5 * self.transient.x[indices][0], 1.2 * self.transient.x[indices][-1]) + if self.transient.magnitude_data: + axes[i].set_ylim(0.8 * min(self.transient.y[indices]), 1.2 * np.max(self.transient.y[indices])) + axes[i].invert_yaxis() + else: + axes[i].set_ylim(0.5 * min(self.transient.y[indices]), 2. * np.max(self.transient.y[indices])) + axes[i].set_yscale("log") + axes[i].legend(ncol=2) + axes[i].tick_params(axis='both', which='major', pad=8) + i += 1 + + figure.supxlabel(xlabel, fontsize=fontsize) + figure.supylabel(ylabel, fontsize=fontsize) + filename = f"{self.transient.name}_{plot_label}.png" + plt.subplots_adjust(wspace=wspace, hspace=hspace) + plt.savefig(join(self.transient.directory_structure.directory_path, filename), bbox_inches="tight") + return axes + + def plot_multiband_lightcurve( + self, model: callable, filename: str = None, axes: matplotlib.axes.Axes = None, plot_save: bool = True, + plot_show: bool = True, random_models: int = 100, posterior: pd.DataFrame = None, outdir: str = '.', + model_kwargs: dict = None, **kwargs: object) -> None: + """ + + Parameters + ---------- + model: callable + The model used to plot the lightcurve + filename: str, optional + The output filename. Otherwise, use default which starts with the name + attribute and ends with *lightcurve.png. + axes: matplotlib.axes.Axes, optional + Axes to plot in if given. + plot_save: bool, optional + Whether to save the plot. + plot_show: bool, optional + Whether to show the plot. + random_models: int, optional + Number of random posterior samples plotted faintly. Default is 100. + posterior: pd.DataFrame, optional + Posterior distribution to which to draw samples from. Is optional but must be given. + outdir: str, optional + Out directory in which to save the plot. Default is the current working directory. + model_kwargs: dict + Additional keyword arguments to be passed into the model. + kwargs: dict + No current function. + ------- + + """ + if self.transient.luminosity_data or self.transient.flux_data: + redback.utils.logger.warning( + f"Plotting multiband lightcurve not possible for {self.transient.data_mode}. Returning.") + return + + if filename is None: + filename = f"{self.transient.name}_multiband_lightcurve.png" + axes = axes or plt.gca() + axes = self.plot_multiband(axes=axes) + + times = self._get_times(axes) + + times_mesh, frequency_mesh = np.meshgrid(times, self.transient.unique_frequencies) + new_model_kwargs = model_kwargs.copy() + new_model_kwargs['frequency'] = frequency_mesh + posterior.sort_values(by='log_likelihood') + max_like_params = posterior.iloc[-1] + ys = model(times_mesh, **max_like_params, **new_model_kwargs) + + for i in range(len(self.transient.unique_frequencies)): + axes[i].plot(times_mesh[i], ys[i], color='blue', alpha=0.65, lw=2) + for _ in range(random_models): + params = posterior.iloc[np.random.randint(len(posterior))] + ys = model(times_mesh, **params, **new_model_kwargs) + axes[i].plot(times, ys[i], color='red', alpha=0.05, lw=2, zorder=-1) + if plot_save: + plt.savefig(join(outdir, filename), dpi=300, bbox_inches="tight") + if plot_show: + plt.show() + plt.clf() diff --git a/redback/transient/afterglow.py b/redback/transient/afterglow.py index 81a3cf421..5d46da30a 100644 --- a/redback/transient/afterglow.py +++ b/redback/transient/afterglow.py @@ -13,7 +13,7 @@ from redback.utils import logger from redback.get_data.directory import afterglow_directory_structure from redback.plotting import \ - MultiBandPlotter, IntegratedFluxPlotter, LuminosityPlotter, FluxDensityPlotter, MagnitudePlotter + IntegratedFluxPlotter, LuminosityPlotter, FluxDensityPlotter, MagnitudePlotter from redback.transient.transient import Transient dirname = os.path.dirname(__file__) @@ -493,10 +493,27 @@ def plot_multiband( if self.data_mode not in ['flux_density', 'magnitude']: raise ValueError( f'You cannot plot multiband data with {self.data_mode} data mode . Why are you doing this?') - mbd = MultiBandPlotter(transient=self) - return mbd.plot_multiband( - figure=figure, axes=axes, ncols=ncols, - nrows=nrows, figsize=figsize, filters=filters, **plot_kwargs) + if self.magnitude_data: + plotter = MagnitudePlotter(transient=self) + elif self.flux_density_data: + plotter = FluxDensityPlotter(transient=self) + else: + return + return plotter.plot_multiband( + figure=figure, axes=axes, ncols=ncols, nrows=nrows, figsize=figsize, filters=filters, **plot_kwargs) + + def plot_multiband_lightcurve( + self, model: callable, filename: str = None, axes: matplotlib.axes.Axes = None, plot_save: bool = True, + plot_show: bool = True, random_models: int = 100, posterior: pd.DataFrame = None, outdir: str = '.', + model_kwargs: dict = None, **kwargs: object) -> None: + + if self.data_mode not in ['flux_density', 'magnitude']: + raise ValueError( + f'You cannot plot multiband data with {self.data_mode} data mode . Why are you doing this?') + return super(Afterglow, self).plot_multiband_lightcurve( + model=model, filename=filename, axes=axes, plot_save=plot_save, plot_show=plot_show, + random_models=random_models, posterior=posterior, outdir=outdir, model_kwargs=model_kwargs, **kwargs) + class SGRB(Afterglow): pass diff --git a/redback/transient/transient.py b/redback/transient/transient.py index 9b56b8b30..10bd98d17 100644 --- a/redback/transient/transient.py +++ b/redback/transient/transient.py @@ -9,7 +9,7 @@ import redback from redback.plotting import \ - MultiBandPlotter, LuminosityPlotter, FluxDensityPlotter, IntegratedFluxPlotter, MagnitudePlotter + LuminosityPlotter, FluxDensityPlotter, IntegratedFluxPlotter, MagnitudePlotter class Transient(object): @@ -458,43 +458,83 @@ def get_colors(filters: Union[np.ndarray, list]) -> matplotlib.colors.Colormap: """ return matplotlib.cm.rainbow(np.linspace(0, 1, len(filters))) - def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k') -> matplotlib.axes.Axes: + def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k', **kwargs: dict) -> matplotlib.axes.Axes: """ - Base function for data plotting. + Plots the Afterglow lightcurve and returns Axes. Parameters ---------- - axes: matplotlib.axes.Axes, optional - Axes can be given if defaults are not satisfying. + axes : Union[matplotlib.axes.Axes, None], optional + Matplotlib axes to plot the lightcurve into. Useful for user specific modifications to the plot. colour: str, optional - Colour for plotted data. + Colour of the data. + kwargs: dict + Additional keyword arguments to pass in the Plotter methods. Returns - ------- - matplotlib.axes.Axes: The user can make additional modifications to the axes. - + ---------- + matplotlib.axes.Axes: The axes with the plot. """ - fig, axes = plt.subplots() - return axes - def plot_multiband(self, axes: matplotlib.axes.Axes = None, colour: str = 'k') -> matplotlib.axes.Axes: + if self.flux_data: + plotter = IntegratedFluxPlotter(transient=self) + elif self.luminosity_data: + plotter = LuminosityPlotter(transient=self) + elif self.flux_density_data: + plotter = FluxDensityPlotter(transient=self) + elif self.magnitude_data: + plotter = MagnitudePlotter(transient=self) + else: + return axes + return plotter.plot_data(axes=axes, colour=colour, **kwargs) + + def plot_multiband( + self, figure: matplotlib.figure.Figure = None, axes: matplotlib.axes.Axes = None, ncols: int = 2, + nrows: int = None, figsize: tuple = None, filters: list = None, **plot_kwargs: dict) \ + -> matplotlib.axes.Axes: """ - Base function for multiband data plotting. Parameters ---------- - axes: matplotlib.axes.Axes - Axes can be given if defaults are not satisfying. - colour: str, optional - Colour for plotted data. + figure: matplotlib.figure.Figure, optional + Figure can be given if defaults are not satisfying + axes: matplotlib.axes.Axes, optional + Axes can be given if defaults are not satisfying + ncols: int, optional + Number of columns to use on the plot. Default is 2. + nrows: int, optional + Number of rows to use on the plot. If None are given this will + be inferred from ncols and the number of filters. + figsize: tuple, optional + Size of the figure. A default based on ncols and nrows will be used if None is given. + filters: list, optional + Which bands to plot. Will use default filters if None is given. + plot_kwargs: + Additional optional plotting kwargs: + wspace: Extra argument for matplotlib.pyplot.subplots_adjust + hspace: Extra argument for matplotlib.pyplot.subplots_adjust + fontsize: Label fontsize + errorbar_fmt: Errorbar format ('fmt' argument in matplotlib.pyplot.errorbar) + colors: colors to be used for the bands + xlabel: Plot xlabel + ylabel: Plot ylabel + plot_label: Addional filename label appended to the default name Returns ------- - matplotlib.axes.Axes: The user can make additional modifications to the axes. """ - fig, axes = plt.subplots() - return axes + if self.data_mode not in ['flux_density', 'magnitude']: + raise ValueError( + f'You cannot plot multiband data with {self.data_mode} data mode . Why are you doing this?') + if self.magnitude_data: + plotter = MagnitudePlotter(transient=self) + elif self.flux_density_data: + plotter = FluxDensityPlotter(transient=self) + else: + return + return plotter.plot_multiband( + figure=figure, axes=axes, ncols=ncols, nrows=nrows, figsize=figsize, filters=filters, **plot_kwargs) def plot_lightcurve( self, model: callable, filename: str = None, axes: matplotlib.axes.Axes = None, plot_save: bool = True, @@ -526,31 +566,21 @@ def plot_lightcurve( kwargs: dict No current function. """ - if filename is None: - filename = f"{self.data_mode}_lightcurve.png" - if model_kwargs is None: - model_kwargs = dict() - axes = axes or plt.gca() - axes = self.plot_data(axes=axes) - axes.set_yscale('log') - plt.semilogy() - times = self._get_times(axes) - - times_mesh, frequency_mesh = np.meshgrid(times, self.unique_frequencies) - new_model_kwargs = model_kwargs.copy() - new_model_kwargs['frequency'] = frequency_mesh - - posterior.sort_values(by='log_likelihood') - max_like_params = posterior.iloc[-1] - ys = model(times, **max_like_params, **new_model_kwargs) - axes.plot(times, ys, color='blue', alpha=0.65, lw=2) - - for _ in range(random_models): - params = posterior.iloc[np.random.randint(len(posterior))] - ys = model(times, **params, **model_kwargs) - axes.plot(times, ys, color='red', alpha=0.05, lw=2, zorder=-1) - plt.savefig(join(outdir, filename), dpi=300, bbox_inches="tight") - plt.clf() + if self.flux_data: + plotter = IntegratedFluxPlotter(transient=self) + elif self.luminosity_data: + plotter = LuminosityPlotter(transient=self) + elif self.flux_density_data: + plotter = FluxDensityPlotter(transient=self) + elif self.magnitude_data: + plotter = MagnitudePlotter(transient=self) + else: + return axes + return plotter.plot_lightcurve( + model=model, filename=filename, axes=axes, plot_save=plot_save, + plot_show=plot_show, random_models=random_models, posterior=posterior, + outdir=outdir, model_kwargs=model_kwargs, **kwargs) + def plot_multiband_lightcurve( self, model: callable, filename: str = None, axes: matplotlib.axes.Axes = None, plot_save: bool = True, @@ -584,58 +614,18 @@ def plot_multiband_lightcurve( ------- """ - if self.luminosity_data or self.flux_data: - redback.utils.logger.warning(f"Plotting multiband lightcurve not possible for {self.data_mode}. Returning.") - return - - if filename is None: - filename = f"{self.name}_multiband_lightcurve.png" - axes = axes or plt.gca() - axes = self.plot_multiband(axes=axes) - - times = self._get_times(axes) - - times_mesh, frequency_mesh = np.meshgrid(times, self.unique_frequencies) - new_model_kwargs = model_kwargs.copy() - new_model_kwargs['frequency'] = frequency_mesh - posterior.sort_values(by='log_likelihood') - max_like_params = posterior.iloc[-1] - ys = model(times_mesh, **max_like_params, **new_model_kwargs) - - for i in range(len(self.unique_frequencies)): - axes[i].plot(times_mesh[i], ys[i], color='blue', alpha=0.65, lw=2) - for _ in range(random_models): - params = posterior.iloc[np.random.randint(len(posterior))] - ys = model(times_mesh, **params, **new_model_kwargs) - axes[i].plot(times, ys[i], color='red', alpha=0.05, lw=2, zorder=-1) - if plot_save: - plt.savefig(join(outdir, filename), dpi=300, bbox_inches="tight") - if plot_show: - plt.show() - plt.clf() - - def _get_times(self, axes: matplotlib.axes.Axes) -> np.ndarray: - """ - - Parameters - ---------- - axes: matplotlib.axes.Axes - The axes used in the plotting procedure. - Returns - ------- - np.ndarray: Linearly or logarithmically scaled time values depending on the y scale used in the plot. - - """ - if isinstance(axes, np.ndarray): - ax = axes[0] - else: - ax = axes - - if ax.get_yscale() == 'linear': - times = np.linspace(self.x[0], self.x[-1], 200) + if self.data_mode not in ['flux_density', 'magnitude']: + raise ValueError( + f'You cannot plot multiband data with {self.data_mode} data mode . Why are you doing this?') + if self.magnitude_data: + plotter = MagnitudePlotter(transient=self) + elif self.flux_density_data: + plotter = FluxDensityPlotter(transient=self) else: - times = np.exp(np.linspace(np.log(self.x[0]), np.log(self.x[-1]), 200)) - return times + return + return plotter.plot_multiband_lightcurve(model=model, filename=filename, axes=axes, plot_save=plot_save, + plot_show=plot_show, random_models=random_models, posterior=posterior, outdir=outdir, + model_kwargs=model_kwargs, **kwargs) class OpticalTransient(Transient): @@ -818,89 +808,3 @@ def _get_transient_dir(self) -> str: transient_dir, _, _ = redback.get_data.directory.open_access_directory_structure(transient=self.name, transient_type=self.__class__.__name__.lower()) return transient_dir - - def plot_data( - self, axes: matplotlib.axes.Axes = None, filters: list = None, plot_others: bool = True, - plot_save: bool = True, **plot_kwargs: dict) -> None: - """ - Plots the data. - - Parameters - ---------- - axes: matplotlib.axes.Axes, optional - Axes can be given if defaults are not satisfying - filters: list, optional - Which bands to plot. Will use default filters if None is given. - plot_others: bool, optional - Plot all bands outside filters in black without label if True. - plot_kwargs: - Additional optional plotting kwargs: - errorbar_fmt: Errorbar format ('fmt' argument in matplotlib.pyplot.errorbar) - colors: colors to be used for the bands - xlabel: Plot xlabel - ylabel: Plot ylabel - plot_label: Additional filename label appended to the default name - """ - if self.flux_data: - plotter = IntegratedFluxPlotter(transient=self) - elif self.luminosity_data: - plotter = LuminosityPlotter(transient=self) - elif self.flux_density_data: - plotter = FluxDensityPlotter(transient=self) - elif self.magnitude_data: - plotter = MagnitudePlotter(transient=self) - else: - return axes - return plotter.plot_data(axes=axes, filters=filters, plot_others=plot_others, plot_save=plot_save, **plot_kwargs) - - - def plot_multiband( - self, figure: matplotlib.figure.Figure = None, axes: matplotlib.axes.Axes = None, ncols: int = 2, - nrows: int = None, figsize: tuple = None, filters: list = None, **plot_kwargs: dict) \ - -> matplotlib.axes.Axes: - """ - - Parameters - ---------- - figure: matplotlib.figure.Figure, optional - Figure can be given if defaults are not satisfying - axes: matplotlib.axes.Axes, optional - Axes can be given if defaults are not satisfying - ncols: int, optional - Number of columns to use on the plot. Default is 2. - nrows: int, optional - Number of rows to use on the plot. If None are given this will - be inferred from ncols and the number of filters. - figsize: tuple, optional - Size of the figure. A default based on ncols and nrows will be used if None is given. - filters: list, optional - Which bands to plot. Will use default filters if None is given. - plot_kwargs: - Additional optional plotting kwargs: - wspace: Extra argument for matplotlib.pyplot.subplots_adjust - hspace: Extra argument for matplotlib.pyplot.subplots_adjust - fontsize: Label fontsize - errorbar_fmt: Errorbar format ('fmt' argument in matplotlib.pyplot.errorbar) - colors: colors to be used for the bands - xlabel: Plot xlabel - ylabel: Plot ylabel - plot_label: Addional filename label appended to the default name - - Returns - ------- - - """ - mbd = MultiBandPlotter(transient=self) - return mbd.plot_multiband( - figure=figure, axes=axes, ncols=ncols, - nrows=nrows, figsize=figsize, filters=filters, **plot_kwargs) - - def plot_lightcurve(self, model, filename=None, axes=None, plot_save=True, plot_show=True, random_models=100, - posterior=None, outdir='.', model_kwargs=None, **kwargs): - - axes = axes or plt.gca() - axes = self.plot_data(axes=axes, plot_save=False) - - super(OpticalTransient, self).plot_lightcurve( - model=model, filename=filename, axes=axes, plot_save=plot_save, plot_show=plot_show, - random_models=random_models, posterior=posterior, outdir=outdir, model_kwargs=model_kwargs, **kwargs) From f151cc557fa16cb22216831f75a202bd97f96e0d Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 10 Mar 2022 21:09:20 +1100 Subject: [PATCH 18/37] Some more rationalisation --- redback/plotting.py | 42 +++++++++++++++++++++++++++++------------- 1 file changed, 29 insertions(+), 13 deletions(-) diff --git a/redback/plotting.py b/redback/plotting.py index 08dbc6dd0..4bdc834f0 100644 --- a/redback/plotting.py +++ b/redback/plotting.py @@ -10,6 +10,9 @@ class Plotter(object): + def __init__(self, transient): + self.transient = transient + def _get_times(self, axes: matplotlib.axes.Axes) -> np.ndarray: """ @@ -36,9 +39,6 @@ def _get_times(self, axes: matplotlib.axes.Axes) -> np.ndarray: class IntegratedFluxPlotter(Plotter): - def __init__(self, transient): - self.transient = transient - @property def xlim_low(self): return 0.5 * self.transient.x[0] @@ -49,6 +49,14 @@ def xlim_high(self): return 2 * self.transient.x[-1] return 2 * (self.transient.x[-1] + self.x_err[1][-1]) + @property + def ylim_low(self): + return 0.5 * min(self.transient.y) + + @property + def ylim_high(self): + return 2. * np.max(self.transient.y) + @property def x_err(self): if self.transient.x_err is not None: @@ -56,6 +64,14 @@ def x_err(self): else: return None + @property + def y_err(self): + return [np.abs(self.transient.y_err[1, :]), self.transient.y_err[0, :]] + + @property + def xlabel(self): + return r"Time since burst [s]" + def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k', plot_save: bool = True, plot_show: bool = True, **kwargs) -> matplotlib.axes.Axes: """ @@ -74,23 +90,27 @@ def plot_data(self, axes: matplotlib.axes.Axes = None, colour: str = 'k', plot_s ---------- matplotlib.axes.Axes: The axes with the plot. """ + xy = kwargs.get("xy", (0.95, 0.9)) + xycoords = kwargs.get("xycoords", "axes fraction") + horizontalalignment = kwargs.get("horizontalalignment", "right") + annotation_size = kwargs.get("annotation_size", 20) - y_err = [np.abs(self.transient.y_err[1, :]), self.transient.y_err[0, :]] ax = axes or plt.gca() - ax.errorbar(self.transient.x, self.transient.y, xerr=self.x_err, yerr=y_err, + ax.errorbar(self.transient.x, self.transient.y, xerr=self.x_err, yerr=self.y_err, fmt='x', c=colour, ms=1, elinewidth=2, capsize=0.) ax.set_xscale('log') ax.set_yscale('log') ax.set_xlim(self.xlim_low, self.xlim_high) - ax.set_ylim(0.5 * min(self.transient.y), 2. * np.max(self.transient.y)) + ax.set_ylim(self.ylim_low, self.ylim_high) - ax.annotate(self.transient.name, xy=(0.95, 0.9), xycoords='axes fraction', - horizontalalignment='right', size=20) + ax.annotate( + self.transient.name, xy=xy, xycoords=xycoords, + horizontalalignment=horizontalalignment, size=annotation_size) - ax.set_xlabel(r'Time since burst [s]') + ax.set_xlabel(self.xlabel) ax.set_ylabel(self.transient.ylabel) ax.tick_params(axis='x', pad=10) @@ -172,9 +192,6 @@ class LuminosityPlotter(IntegratedFluxPlotter): class MagnitudePlotter(Plotter): - def __init__(self, transient): - self.transient = transient - def plot_data( self, axes: matplotlib.axes.Axes = None, filters: list = None, plot_others: bool = True, plot_save: bool = True, **plot_kwargs: dict) -> None: @@ -415,7 +432,6 @@ def plot_multiband_lightcurve( plt.clf() - class FluxDensityPlotter(MagnitudePlotter): def _set_y_axis(self, ax): From ac80860856040d9a78bfeaf6c346dd6ed413cb96 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 10 Mar 2022 21:23:52 +1100 Subject: [PATCH 19/37] Removed datamode kwargs from fit_model --- examples/SN2011kl_sample_in_t0_example.py | 2 +- examples/broadband_afterglow_private_data_example.py | 2 +- examples/fit_your_own_model_example.py | 2 +- examples/kilonova_example.py | 2 +- examples/magnetar_boosted_example.py | 2 +- examples/prompt_example.py | 2 +- examples/supernova_example.py | 2 +- examples/tde_example.py | 4 ++-- 8 files changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/SN2011kl_sample_in_t0_example.py b/examples/SN2011kl_sample_in_t0_example.py index c1118fe51..b0043271f 100644 --- a/examples/SN2011kl_sample_in_t0_example.py +++ b/examples/SN2011kl_sample_in_t0_example.py @@ -46,7 +46,7 @@ # returns a supernova result object result = redback.fit_model(name=sne, transient=supernova, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, data_mode='flux_density', sample='rslice', nlive=200, resume=False) + prior=priors, sample='rslice', nlive=200, resume=False) result.plot_corner() diff --git a/examples/broadband_afterglow_private_data_example.py b/examples/broadband_afterglow_private_data_example.py index c3ceeb1c1..8de528b1d 100644 --- a/examples/broadband_afterglow_private_data_example.py +++ b/examples/broadband_afterglow_private_data_example.py @@ -49,7 +49,7 @@ # returns a supernova result object result = redback.fit_model(name=name, transient=afterglow, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, sample='rslice', data_mode='flux_density', nlive=200, resume=True) + prior=priors, sample='rslice', nlive=200, resume=True) # plot corner result.plot_corner() diff --git a/examples/fit_your_own_model_example.py b/examples/fit_your_own_model_example.py index c67c621bf..f3419f461 100644 --- a/examples/fit_your_own_model_example.py +++ b/examples/fit_your_own_model_example.py @@ -42,6 +42,6 @@ def my_favourite_model(time, l0, alpha, **kwargs): # Call redback.fit_model to run the sampler and obtain GRB result object result = redback.fit_model(name=GRB, model=model, sampler='dynesty', nlive=200, transient=afterglow, - prior=priors, data_mode='luminosity', sample='rslice', resume=False) + prior=priors, sample='rslice', resume=False) result.plot_lightcurve(random_models=1000, model=my_favourite_model) diff --git a/examples/kilonova_example.py b/examples/kilonova_example.py index b0166b7c3..0773696bf 100644 --- a/examples/kilonova_example.py +++ b/examples/kilonova_example.py @@ -26,7 +26,7 @@ model_kwargs = dict(frequency=kilonova.filtered_frequencies, output_format='flux_density') result = redback.fit_model(name=kne, transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, data_mode='flux_density', sample='rslice', nlive=200, resume=True) + prior=priors, sample='rslice', nlive=200, resume=True) result.plot_corner() # returns a Kilonova result object result.plot_lightcurve(random_models=1000) diff --git a/examples/magnetar_boosted_example.py b/examples/magnetar_boosted_example.py index 0b6eaa43c..27b4a2f77 100644 --- a/examples/magnetar_boosted_example.py +++ b/examples/magnetar_boosted_example.py @@ -30,7 +30,7 @@ model_kwargs = dict(frequencies=redback.utils.bands_to_frequency(bands), output_format='flux_density') result = redback.fit_model(name=kne, transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, - path=path, prior=priors, data_mode='flux_density', sample='rslice', nlive=200) + path=path, prior=priors, sample='rslice', nlive=200) result.plot_corner() # returns a Kilonova result object result.plot_lightcurve(random_models=1000) diff --git a/examples/prompt_example.py b/examples/prompt_example.py index 86c94ea68..17bf97470 100644 --- a/examples/prompt_example.py +++ b/examples/prompt_example.py @@ -20,7 +20,7 @@ y=prompt.counts, yerr=prompt.counts_err, dt=prompt.bin_size) result = redback.fit_model(source_type='prompt', name=name, model=model, transient=prompt, nlive=500, - sampler=sampler, prior=priors, data_mode='counts', outdir="GRB_results", sample='rslice') + sampler=sampler, prior=priors, outdir="GRB_results", sample='rslice') # returns a GRB prompt result object result.plot_lightcurve(random_models=1000) result.plot_corner() diff --git a/examples/supernova_example.py b/examples/supernova_example.py index c9cb342ca..45a7ceaca 100644 --- a/examples/supernova_example.py +++ b/examples/supernova_example.py @@ -23,6 +23,6 @@ # returns a supernova result object result = redback.fit_model(name=sne, transient=supernova, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, data_mode='flux_density', sample='rslice', nlive=200, resume=False) + prior=priors, sample='rslice', nlive=200, resume=False) result.plot_corner() result.plot_lightcurve(random_models=1000) diff --git a/examples/tde_example.py b/examples/tde_example.py index 8e4375ffa..aed5eacec 100644 --- a/examples/tde_example.py +++ b/examples/tde_example.py @@ -28,8 +28,8 @@ model_kwargs = dict(frequencies=redback.utils.bands_to_frequency(bands), output_format='flux_density') # returns a tde result object -result = redback.fit_model(name=tde, transient=tidal_disruption_event, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, data_mode='flux_density', sample='rslice', nlive=200, resume=False) +result = redback.fit_model(name=tde, transient=tidal_disruption_event, model=model, sampler=sampler, + model_kwargs=model_kwargs, prior=priors, sample='rslice', nlive=200, resume=False) result.plot_corner() From 15bc84bd83919a6d61269e37cd78d7d74ba15144 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 10 Mar 2022 21:26:23 +1100 Subject: [PATCH 20/37] plot data now after flux to luminosity conversion --- examples/fit_your_own_model_example.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fit_your_own_model_example.py b/examples/fit_your_own_model_example.py index f3419f461..997dff3ff 100644 --- a/examples/fit_your_own_model_example.py +++ b/examples/fit_your_own_model_example.py @@ -28,11 +28,11 @@ def my_favourite_model(time, l0, alpha, **kwargs): afterglow = redback.afterglow.SGRB.from_swift_grb(name=GRB, data_mode='flux', truncate=True, truncate_method="prompt_time_error") -afterglow.plot_data() # uses an analytical k-correction expression to create luminosity data if not already there. # Can also use a numerical k-correction through CIAO afterglow.analytical_flux_to_luminosity() +afterglow.plot_data() # You need to create your own priors for this new model. # The model has two parameters l0 and alpha. We use bilby priors for this From d076fdf6dacf6f252f3c263339e3292b48db6e27 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 10 Mar 2022 21:28:06 +1100 Subject: [PATCH 21/37] removed path --- examples/magnetar_boosted_example.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/magnetar_boosted_example.py b/examples/magnetar_boosted_example.py index 27b4a2f77..75ecdba8a 100644 --- a/examples/magnetar_boosted_example.py +++ b/examples/magnetar_boosted_example.py @@ -9,7 +9,6 @@ model = 'mergernova' kne = 'at2017gfo' -path = 'KNDir' # gets the magnitude data for AT2017gfo, the KN associated with GW170817 data = redback.get_data.get_kilonova_data_from_open_transient_catalog_data(transient=kne) # creates a GRBDir with GRB @@ -30,7 +29,7 @@ model_kwargs = dict(frequencies=redback.utils.bands_to_frequency(bands), output_format='flux_density') result = redback.fit_model(name=kne, transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, - path=path, prior=priors, sample='rslice', nlive=200) + prior=priors, sample='rslice', nlive=200) result.plot_corner() # returns a Kilonova result object result.plot_lightcurve(random_models=1000) From cbe1ea8f7e64eb77716e64e9adb9d88602b7f605 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 10 Mar 2022 21:30:28 +1100 Subject: [PATCH 22/37] frequencies --> frequency --- examples/SN2011kl_sample_in_t0_example.py | 2 +- examples/magnetar_boosted_example.py | 2 +- examples/supernova_example.py | 2 +- examples/tde_example.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/SN2011kl_sample_in_t0_example.py b/examples/SN2011kl_sample_in_t0_example.py index b0043271f..d546947c5 100644 --- a/examples/SN2011kl_sample_in_t0_example.py +++ b/examples/SN2011kl_sample_in_t0_example.py @@ -42,7 +42,7 @@ priors['av'] = Uniform(0.1, 1, name='av', latex_label=r'$a_{v}$') -model_kwargs = dict(frequencies=redback.utils.bands_to_frequency(bands), output_format='flux_density') +model_kwargs = dict(frequency=redback.utils.bands_to_frequency(bands), output_format='flux_density') # returns a supernova result object result = redback.fit_model(name=sne, transient=supernova, model=model, sampler=sampler, model_kwargs=model_kwargs, diff --git a/examples/magnetar_boosted_example.py b/examples/magnetar_boosted_example.py index 75ecdba8a..8b4bf1d93 100644 --- a/examples/magnetar_boosted_example.py +++ b/examples/magnetar_boosted_example.py @@ -26,7 +26,7 @@ priors['tau_sd'] = 1e3 priors['thermalisation_efficiency'] = 0.3 -model_kwargs = dict(frequencies=redback.utils.bands_to_frequency(bands), output_format='flux_density') +model_kwargs = dict(frequency=redback.utils.bands_to_frequency(bands), output_format='flux_density') result = redback.fit_model(name=kne, transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, prior=priors, sample='rslice', nlive=200) diff --git a/examples/supernova_example.py b/examples/supernova_example.py index 45a7ceaca..e87d254ef 100644 --- a/examples/supernova_example.py +++ b/examples/supernova_example.py @@ -19,7 +19,7 @@ # use default priors priors = redback.priors.get_priors(model=model) priors['redshift'] = 0.677 -model_kwargs = dict(frequencies=redback.utils.bands_to_frequency(bands), output_format='flux_density') +model_kwargs = dict(frequency=redback.utils.bands_to_frequency(bands), output_format='flux_density') # returns a supernova result object result = redback.fit_model(name=sne, transient=supernova, model=model, sampler=sampler, model_kwargs=model_kwargs, diff --git a/examples/tde_example.py b/examples/tde_example.py index aed5eacec..d3a8a6673 100644 --- a/examples/tde_example.py +++ b/examples/tde_example.py @@ -25,7 +25,7 @@ # We can do this through the metadata that was downloaded alongside the data, or if you just know it. priors['redshift'] = tidal_disruption_event.redshift -model_kwargs = dict(frequencies=redback.utils.bands_to_frequency(bands), output_format='flux_density') +model_kwargs = dict(frequency=redback.utils.bands_to_frequency(bands), output_format='flux_density') # returns a tde result object result = redback.fit_model(name=tde, transient=tidal_disruption_event, model=model, sampler=sampler, From 4c45364e701a5d7dfaaa7e5da20b73ecb0482fb0 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 10 Mar 2022 22:51:31 +1100 Subject: [PATCH 23/37] Fixed kilonova example (again) --- ...roadband_afterglow_private_data_example.py | 2 +- examples/kilonova_example.py | 15 +-- redback/plotting.py | 99 +++++++++++++++---- 3 files changed, 91 insertions(+), 25 deletions(-) diff --git a/examples/broadband_afterglow_private_data_example.py b/examples/broadband_afterglow_private_data_example.py index 8de528b1d..e92f8dbfa 100644 --- a/examples/broadband_afterglow_private_data_example.py +++ b/examples/broadband_afterglow_private_data_example.py @@ -49,7 +49,7 @@ # returns a supernova result object result = redback.fit_model(name=name, transient=afterglow, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, sample='rslice', nlive=200, resume=True) + prior=priors, sample='rslice', nlive=50, dlogz=10, resume=True) # plot corner result.plot_corner() diff --git a/examples/kilonova_example.py b/examples/kilonova_example.py index 0773696bf..36d13dd7e 100644 --- a/examples/kilonova_example.py +++ b/examples/kilonova_example.py @@ -15,18 +15,19 @@ # creates a GRBDir with GRB kilonova = redback.kilonova.Kilonova.from_open_access_catalogue( name=kne, data_mode="flux_density", active_bands=np.array(["g", "i"])) -# kilonova.flux_density_data = True +kilonova.plot_data(plot_show=False) fig, axes = plt.subplots(3, 2, sharex=True, sharey=True, figsize=(12, 8)) -kilonova.plot_multiband(figure=fig, axes=axes, filters=["g", "r", "i", "z", "y", "J"]) +kilonova.plot_multiband(figure=fig, axes=axes, filters=["g", "r", "i", "z", "y", "J"], plot_show=False) # use default priors priors = redback.priors.get_priors(model=model) priors['redshift'] = 1e-2 model_kwargs = dict(frequency=kilonova.filtered_frequencies, output_format='flux_density') - -result = redback.fit_model(name=kne, transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, sample='rslice', nlive=200, resume=True) -result.plot_corner() +# result = redback.fit_model(name=kne, transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, +# prior=priors, sample='rslice', nlive=200, resume=True) +result = redback.result.read_in_result("kilonova/one_component_kilonova_model/at2017gfo_result.json") +# result.plot_corner() # returns a Kilonova result object -result.plot_lightcurve(random_models=1000) +result.plot_lightcurve(plot_show=False) +result.plot_multiband_lightcurve(filters=["g", "r", "i", "z", "y", "J"], plot_show=False) diff --git a/redback/plotting.py b/redback/plotting.py index 4bdc834f0..dd72493d2 100644 --- a/redback/plotting.py +++ b/redback/plotting.py @@ -36,9 +36,6 @@ def _get_times(self, axes: matplotlib.axes.Axes) -> np.ndarray: times = np.exp(np.linspace(np.log(self.xlim_low), np.log(self.xlim_high), 200)) return times - -class IntegratedFluxPlotter(Plotter): - @property def xlim_low(self): return 0.5 * self.transient.x[0] @@ -68,6 +65,9 @@ def x_err(self): def y_err(self): return [np.abs(self.transient.y_err[1, :]), self.transient.y_err[0, :]] + +class IntegratedFluxPlotter(Plotter): + @property def xlabel(self): return r"Time since burst [s]" @@ -183,7 +183,7 @@ def plot_lightcurve( if plot_show: plt.tight_layout() plt.show() - plt.clf() + return axes class LuminosityPlotter(IntegratedFluxPlotter): @@ -193,8 +193,8 @@ class LuminosityPlotter(IntegratedFluxPlotter): class MagnitudePlotter(Plotter): def plot_data( - self, axes: matplotlib.axes.Axes = None, filters: list = None, plot_others: bool = True, - plot_save: bool = True, **plot_kwargs: dict) -> None: + self, axes: matplotlib.axes.Axes = None, filters: list = None, plot_others: bool = False, + plot_save: bool = True, plot_show: bool = True, **plot_kwargs: dict) -> None: """ Plots the data. @@ -253,19 +253,79 @@ def plot_data( if plot_save: filename = f"{self.transient.name}_{plot_label}.png" plt.savefig(join(self.transient.directory_structure.directory_path, filename), bbox_inches='tight') - plt.clf() + if plot_show: + plt.tight_layout() + plt.show() return axes def _set_y_axis(self, ax): ax.set_ylim(0.8 * min(self.transient.y), 1.2 * np.max(self.transient.y)) ax.invert_yaxis() - def plot_lightcurve(self): - pass + def plot_lightcurve( + self, model: callable, filename: str = None, axes: matplotlib.axes.Axes = None, plot_save: bool = True, + plot_show: bool = True, random_models: int = 100, posterior: pd.DataFrame = None, outdir: str = '.', + model_kwargs: dict = None, **kwargs: object) -> None: + """ + + Parameters + ---------- + model: callable + The model used to plot the lightcurve. + filename: str, optional + The output filename. Otherwise, use default which starts with the name + attribute and ends with *lightcurve.png. + axes: matplotlib.axes.Axes, optional + Axes to plot in if given. + plot_save: bool, optional + Whether to save the plot. + plot_show: bool, optional + Whether to show the plot. + random_models: int, optional + Number of random posterior samples plotted faintly. Default is 100. + posterior: pd.DataFrame, optional + Posterior distribution to which to draw samples from. Is optional but must be given. + outdir: str, optional + Out directory in which to save the plot. Default is the current working directory. + model_kwargs: dict + Additional keyword arguments to be passed into the model. + kwargs: dict + No current function. + """ + if filename is None: + filename = f"{self.transient.name}_lightcurve.png" + if model_kwargs is None: + model_kwargs = dict() + axes = axes or plt.gca() + axes = self.plot_data(axes=axes, plot_save=False, plot_show=False) + axes.set_yscale('log') + # plt.semilogy() + times = self._get_times(axes) + + posterior.sort_values(by='log_likelihood') + max_like_params = posterior.iloc[-1] + random_params_list = [posterior.iloc[np.random.randint(len(posterior))] for _ in range(random_models)] + + for band, color in zip(self.transient.active_bands, self.transient.get_colors(self.transient.active_bands)): + frequency = redback.utils.bands_to_frequency([band]) + model_kwargs["frequency"] = np.ones(len(times)) * frequency + ys = model(times, **max_like_params, **model_kwargs) + axes.plot(times, ys, color=color, alpha=0.65, lw=2) + + for params in random_params_list: + ys = model(times, **params, **model_kwargs) + axes.plot(times, ys, color='red', alpha=0.05, lw=2, zorder=-1) + if plot_save: + plt.savefig(join(outdir, filename), dpi=300, bbox_inches="tight") + if plot_show: + plt.tight_layout() + plt.show() + return axes def plot_multiband( self, figure: matplotlib.figure.Figure = None, axes: matplotlib.axes.Axes = None, ncols: int = 2, - nrows: int = None, figsize: tuple = None, filters: list = None, **plot_kwargs: dict) -> \ + nrows: int = None, figsize: tuple = None, filters: list = None, plot_save: bool=True, plot_show: bool=True, + **plot_kwargs: dict) -> \ matplotlib.axes.Axes: """ @@ -325,7 +385,7 @@ def plot_multiband( raise ValueError(f"Insufficient number of panels. {npanels} panels were given " f"but {len(filters)} panels are needed.") if figsize is None: - figsize = (6 * nrows, 4 * ncols) + figsize = (4 + 4 * ncols, 2 + 2 * nrows) figure, axes = plt.subplots(ncols=ncols, nrows=nrows, sharex='all', figsize=figsize) axes = axes.ravel() @@ -365,7 +425,11 @@ def plot_multiband( figure.supylabel(ylabel, fontsize=fontsize) filename = f"{self.transient.name}_{plot_label}.png" plt.subplots_adjust(wspace=wspace, hspace=hspace) - plt.savefig(join(self.transient.directory_structure.directory_path, filename), bbox_inches="tight") + if plot_save: + plt.savefig(join(self.transient.directory_structure.directory_path, filename), bbox_inches="tight") + if plot_show: + plt.tight_layout() + plt.show() return axes def plot_multiband_lightcurve( @@ -408,19 +472,20 @@ def plot_multiband_lightcurve( if filename is None: filename = f"{self.transient.name}_multiband_lightcurve.png" axes = axes or plt.gca() - axes = self.plot_multiband(axes=axes) + filters = kwargs.get("filters", self.transient.active_bands) + axes = self.plot_multiband(axes=axes, plot_save=False, plot_show=False, filters=filters) times = self._get_times(axes) - times_mesh, frequency_mesh = np.meshgrid(times, self.transient.unique_frequencies) + times_mesh, frequency_mesh = np.meshgrid(times, redback.utils.bands_to_frequency(filters)) new_model_kwargs = model_kwargs.copy() new_model_kwargs['frequency'] = frequency_mesh posterior.sort_values(by='log_likelihood') max_like_params = posterior.iloc[-1] ys = model(times_mesh, **max_like_params, **new_model_kwargs) - for i in range(len(self.transient.unique_frequencies)): - axes[i].plot(times_mesh[i], ys[i], color='blue', alpha=0.65, lw=2) + for i in range(len(filters)): + axes[i].plot(times, ys[i], color='blue', alpha=0.65, lw=2) for _ in range(random_models): params = posterior.iloc[np.random.randint(len(posterior))] ys = model(times_mesh, **params, **new_model_kwargs) @@ -429,7 +494,7 @@ def plot_multiband_lightcurve( plt.savefig(join(outdir, filename), dpi=300, bbox_inches="tight") if plot_show: plt.show() - plt.clf() + return axes class FluxDensityPlotter(MagnitudePlotter): From fccaa55ef82cbb4d40e3fe0a212e75bbdbd4aa47 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 10 Mar 2022 22:59:09 +1100 Subject: [PATCH 24/37] Remove apostrophe's --- redback/get_data/open_data.py | 1 + 1 file changed, 1 insertion(+) diff --git a/redback/get_data/open_data.py b/redback/get_data/open_data.py index b88b632c1..63daf5ed3 100644 --- a/redback/get_data/open_data.py +++ b/redback/get_data/open_data.py @@ -139,6 +139,7 @@ def convert_raw_data_to_csv(self) -> Union[pd.DataFrame, None]: magnitude_error=data['e_magnitude'].values, reference_flux=3631, magnitude_system='AB') + data['band'] = [b.replace("'", "") for b in data["band"]] metadata = pd.read_csv(f"{self.directory_path}{self.transient}_metadata.csv") metadata.replace(r'^\s+$', np.nan, regex=True) time_of_event = self.get_time_of_event(data=data, metadata=metadata) From e1b51ff893535149c6785ec3d5e60aff72179a16 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Thu, 10 Mar 2022 23:31:53 +1100 Subject: [PATCH 25/37] Fixed some priors --- redback/priors/arnett.prior | 2 +- redback/priors/arnett_bolometric.prior | 2 +- redback/priors/basic_magnetar_powered.prior | 4 ++-- .../basic_magnetar_powered_bolometric.prior | 2 +- redback/priors/magnetar_nickel.prior | 2 +- test/prior_test.py | 19 ++++++++++++++----- 6 files changed, 20 insertions(+), 11 deletions(-) diff --git a/redback/priors/arnett.prior b/redback/priors/arnett.prior index 44e8ca322..1dd25117a 100644 --- a/redback/priors/arnett.prior +++ b/redback/priors/arnett.prior @@ -1,5 +1,5 @@ redshift = Uniform(1e-6,3,name='redshift', latex_label = r'$z$') -f_nickel = LogUniform(1e-3,1,name='f_nickel', latex_label = r'$f_{\mathrm{Ni}}') +f_nickel = LogUniform(1e-3,1,name='f_nickel', latex_label = r'$f_{\mathrm{Ni}}$') mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$') vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') diff --git a/redback/priors/arnett_bolometric.prior b/redback/priors/arnett_bolometric.prior index 4a38085f7..639ec8109 100644 --- a/redback/priors/arnett_bolometric.prior +++ b/redback/priors/arnett_bolometric.prior @@ -1,4 +1,4 @@ -f_nickel = LogUniform(1e-3,1,name='f_nickel', latex_label = r'$f_{\mathrm{Ni}}') +f_nickel = LogUniform(1e-3,1,name='f_nickel', latex_label = r'$f_{\mathrm{Ni}}$') mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$') vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') diff --git a/redback/priors/basic_magnetar_powered.prior b/redback/priors/basic_magnetar_powered.prior index 8c22038e6..7b545707f 100644 --- a/redback/priors/basic_magnetar_powered.prior +++ b/redback/priors/basic_magnetar_powered.prior @@ -1,5 +1,5 @@ redshift = Uniform(1e-6,3,name='redshift', latex_label = r'$z$') -p0 = Uniform(1, 10, 'p0', latex_label = r'$P_{0]$ [ms]') +p0 = Uniform(1, 10, 'p0', latex_label = r'$P_{0}$ [ms]') bp = LogUniform(0.1,10,'bp',latex_label = r'$B_{p}$~[10$^{14}$G]') mass_ns = Uniform(1.1, 2.2, 'mass_ns', latex_label = r'$M_{\mathrm{NS}} [M_{\odot}]$') theta_pb = Uniform(0, 3.14/2, 'theta_pb', latex_label = r'$\\theta_{P-B}$') @@ -7,4 +7,4 @@ mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$ vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') kappa_gamma = Uniform(1e-4, 1e4, 'kappa_gamma', latex_label = r'$\kappa_{\gamma}$') -temperature_floor = LogUniform(1e3,1e5,name = 'temperature_floor', latex_label = r'$T_{\mathrm{floor}$ [k]}') \ No newline at end of file +temperature_floor = LogUniform(1e3,1e5,name = 'temperature_floor', latex_label = r'$T_{\mathrm{floor}}$ [k]') \ No newline at end of file diff --git a/redback/priors/basic_magnetar_powered_bolometric.prior b/redback/priors/basic_magnetar_powered_bolometric.prior index 0af025c5e..8114adf0e 100644 --- a/redback/priors/basic_magnetar_powered_bolometric.prior +++ b/redback/priors/basic_magnetar_powered_bolometric.prior @@ -1,4 +1,4 @@ -p0 = Uniform(1, 10, 'p0', latex_label = r'$P_{0]$ [ms]') +p0 = Uniform(1, 10, 'p0', latex_label = r'$P_{0}$ [ms]') bp = LogUniform(0.1,10,'bp',latex_label = r'$B_{p}$~[10$^{14}$G]') mass_ns = Uniform(1.1, 2.2, 'mass_ns', latex_label = r'$M_{\mathrm{NS}} [M_{\odot}]$') theta_pb = Uniform(0, 3.14/2, 'theta_pb', latex_label = r'$\\theta_{P-B}$') diff --git a/redback/priors/magnetar_nickel.prior b/redback/priors/magnetar_nickel.prior index 40ce59300..371fc3c8e 100644 --- a/redback/priors/magnetar_nickel.prior +++ b/redback/priors/magnetar_nickel.prior @@ -1,5 +1,5 @@ redshift = Uniform(1e-6,3,name='redshift', latex_label = r'$z$') -f_nickel = LogUniform(1e-3,1,name='f_nickel', latex_label = r'$f_{\mathrm{Ni}}') +f_nickel = LogUniform(1e-3,1,name='f_nickel', latex_label = r'$f_{\mathrm{Ni}}$') p0 = Uniform(1, 10, 'p0', latex_label = r'$P_{0]$ [ms]') bp = LogUniform(0.1,10,'bp',latex_label = r'$B_{p}$~[10$^{14}$G]') mass_ns = Uniform(1.1, 2.2, 'mass_ns', latex_label = r'$M_{\mathrm{NS}} [M_{\odot}]$') diff --git a/test/prior_test.py b/test/prior_test.py index 494de699b..a46779d35 100644 --- a/test/prior_test.py +++ b/test/prior_test.py @@ -6,16 +6,25 @@ class TestLoadPriors(unittest.TestCase): def setUp(self) -> None: - pass + self.path_to_files = "../redback/priors/" + self.prior_files = listdir(self.path_to_files) def tearDown(self) -> None: pass def test_load_priors(self): - path_to_files = "../redback/priors/" - prior_files = listdir(path_to_files) + for f in self.prior_files: + prior_dict = bilby.prior.PriorDict() + prior_dict.from_file(f"{self.path_to_files}{f}") - for f in prior_files: + def test_dollar_signs(self): + for f in self.prior_files: + print() print(f) prior_dict = bilby.prior.PriorDict() - prior_dict.from_file(f"{path_to_files}{f}") + prior_dict.from_file(f"{self.path_to_files}{f}") + for k, p in prior_dict.items(): + print(k) + occurences = p.latex_label.count("$") + assert occurences % 2 == 0 + From 3b20c240eebe819f3f37c0e385e66a71e9661062 Mon Sep 17 00:00:00 2001 From: Moritz Huebner Date: Fri, 11 Mar 2022 00:25:40 +1100 Subject: [PATCH 26/37] Fixed examples in various ways --- examples/SN2011kl_sample_in_t0_example.py | 21 ++++++++----------- ...roadband_afterglow_private_data_example.py | 2 +- examples/fit_your_own_model_example.py | 2 +- examples/kilonova_example.py | 2 +- examples/magnetar_boosted_example.py | 9 ++++---- examples/magnetar_example.py | 2 +- examples/prompt_example.py | 2 +- examples/supernova_example.py | 12 +++++------ examples/tde_example.py | 13 ++++++------ redback/plotting.py | 5 ++++- redback/priors/arnett.prior | 2 +- redback/priors/tde_analytical.prior | 8 +++++++ redback/sampler.py | 10 ++++----- redback/sed.py | 15 +++++++------ redback/transient/transient.py | 8 +++++-- redback/transient_models/extinction_models.py | 7 ++++--- redback/transient_models/supernova_models.py | 15 +++++++------ redback/transient_models/tde_models.py | 19 +++++++++-------- 18 files changed, 87 insertions(+), 67 deletions(-) create mode 100644 redback/priors/tde_analytical.prior diff --git a/examples/SN2011kl_sample_in_t0_example.py b/examples/SN2011kl_sample_in_t0_example.py index d546947c5..39e957a94 100644 --- a/examples/SN2011kl_sample_in_t0_example.py +++ b/examples/SN2011kl_sample_in_t0_example.py @@ -1,5 +1,3 @@ -import matplotlib.pyplot as plt - import redback from bilby.core.prior import Uniform, Gaussian @@ -9,18 +7,17 @@ # we want to sample in t0 and with extinction so we use model = 't0_supernova_extinction' # we want to specify a base model for the actual physics -base_model = "exponential_powerlaw" +base_model = "arnett" data = redback.get_data.get_supernova_data_from_open_transient_catalog_data(sne) # load the data into a supernova transient object which does all the processing and can be used to make plots etc # we set the data_mode to flux density to use/fit flux density. We could use 'magnitude' or 'luminosity' or flux as well. # However, for optical transients we recommend fitting in flux_density. -supernova = redback.supernova.Supernova.from_open_access_catalogue(name=sne, data_mode='flux_density') +supernova = redback.supernova.Supernova.from_open_access_catalogue(name=sne, data_mode='flux_density', use_phase_model=True) # lets make a plot of the data -fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(12, 8)) -supernova.plot_multiband(figure=fig, axes=axes, filters=["J", "H", "g", "i"]) +supernova.plot_multiband(filters=["J", "H", "g", "i"]) # we are now only going to fit g and i bands. We set this using the transient object and the active bands attribute bands = ["g", "i"] @@ -28,6 +25,7 @@ # use default priors priors = redback.priors.get_priors(model=model) +priors.update(redback.priors.get_priors(model=base_model)) # we know the redshift for SN2011kl so we just fix the prior for the redshift to the known value priors['redshift'] = 0.677 @@ -35,19 +33,18 @@ # we use bilby prior objects for this # let's set t0 as a Gaussian around the first observation with sigma = 0.5 - There are many other priors to choose from. # This also demonstates how one can change the latex labels etc -priors['t0'] = Gaussian(data['time [mjd]'], sigma=0.5, name='t0', latex_label=r'$T_{\rm{0}}$') +priors['t0'] = Gaussian(data['time'][0], sigma=0.5, name='t0', latex_label=r'$T_{\rm{0}}$') # We also need a prior on A_v i.e., the total mag extinction. # Just use uniform prior for now. priors['av'] = Uniform(0.1, 1, name='av', latex_label=r'$a_{v}$') - -model_kwargs = dict(frequency=redback.utils.bands_to_frequency(bands), output_format='flux_density') +model_kwargs = dict(frequency=supernova.filtered_frequencies, output_format='flux_density', base_model=base_model) # returns a supernova result object -result = redback.fit_model(name=sne, transient=supernova, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, sample='rslice', nlive=200, resume=False) +result = redback.fit_model(transient=supernova, model=model, sampler=sampler, model_kwargs=model_kwargs, + prior=priors, sample='rslice', nlive=200, resume=True) result.plot_corner() -result.plot_lightcurve(random_models=1000) +result.plot_lightcurve(random_models=100) diff --git a/examples/broadband_afterglow_private_data_example.py b/examples/broadband_afterglow_private_data_example.py index e92f8dbfa..853f194db 100644 --- a/examples/broadband_afterglow_private_data_example.py +++ b/examples/broadband_afterglow_private_data_example.py @@ -48,7 +48,7 @@ model_kwargs = dict(frequency=frequency, output_format='flux_density') # returns a supernova result object -result = redback.fit_model(name=name, transient=afterglow, model=model, sampler=sampler, model_kwargs=model_kwargs, +result = redback.fit_model(transient=afterglow, model=model, sampler=sampler, model_kwargs=model_kwargs, prior=priors, sample='rslice', nlive=50, dlogz=10, resume=True) # plot corner result.plot_corner() diff --git a/examples/fit_your_own_model_example.py b/examples/fit_your_own_model_example.py index 997dff3ff..09a33d59c 100644 --- a/examples/fit_your_own_model_example.py +++ b/examples/fit_your_own_model_example.py @@ -41,7 +41,7 @@ def my_favourite_model(time, l0, alpha, **kwargs): priors['alpha'] = Uniform(-7, 0, 'alpha', latex_label=r'$\alpha$') # Call redback.fit_model to run the sampler and obtain GRB result object -result = redback.fit_model(name=GRB, model=model, sampler='dynesty', nlive=200, transient=afterglow, +result = redback.fit_model(model=model, sampler='dynesty', nlive=200, transient=afterglow, prior=priors, sample='rslice', resume=False) result.plot_lightcurve(random_models=1000, model=my_favourite_model) diff --git a/examples/kilonova_example.py b/examples/kilonova_example.py index 36d13dd7e..49d8a6f7f 100644 --- a/examples/kilonova_example.py +++ b/examples/kilonova_example.py @@ -24,7 +24,7 @@ priors['redshift'] = 1e-2 model_kwargs = dict(frequency=kilonova.filtered_frequencies, output_format='flux_density') -# result = redback.fit_model(name=kne, transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, +# result = redback.fit_model(transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, # prior=priors, sample='rslice', nlive=200, resume=True) result = redback.result.read_in_result("kilonova/one_component_kilonova_model/at2017gfo_result.json") # result.plot_corner() diff --git a/examples/magnetar_boosted_example.py b/examples/magnetar_boosted_example.py index 8b4bf1d93..ab4331739 100644 --- a/examples/magnetar_boosted_example.py +++ b/examples/magnetar_boosted_example.py @@ -28,8 +28,9 @@ model_kwargs = dict(frequency=redback.utils.bands_to_frequency(bands), output_format='flux_density') -result = redback.fit_model(name=kne, transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, sample='rslice', nlive=200) -result.plot_corner() +result = redback.fit_model(transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, + prior=priors, sample='rslice', nlive=200, resume=True) +# result.plot_corner() # returns a Kilonova result object -result.plot_lightcurve(random_models=1000) +# result.plot_lightcurve(random_models=1000) +result.plot_multiband_lightcurve() diff --git a/examples/magnetar_example.py b/examples/magnetar_example.py index acf7a623d..1d43eb2da 100644 --- a/examples/magnetar_example.py +++ b/examples/magnetar_example.py @@ -34,7 +34,7 @@ # Call redback.fit_model to run the sampler and obtain GRB result object -result = redback.fit_model(name=GRB, model=model, sampler='dynesty', nlive=200, transient=afterglow, +result = redback.fit_model(model=model, sampler='dynesty', nlive=200, transient=afterglow, prior=priors, sample='rslice', resume=True) result.plot_lightcurve(random_models=1000) diff --git a/examples/prompt_example.py b/examples/prompt_example.py index 17bf97470..a93ce504e 100644 --- a/examples/prompt_example.py +++ b/examples/prompt_example.py @@ -19,7 +19,7 @@ priors = redback.priors.get_priors(model=model, data_mode='counts', times=prompt.time, y=prompt.counts, yerr=prompt.counts_err, dt=prompt.bin_size) -result = redback.fit_model(source_type='prompt', name=name, model=model, transient=prompt, nlive=500, +result = redback.fit_model(source_type='prompt', model=model, transient=prompt, nlive=500, sampler=sampler, prior=priors, outdir="GRB_results", sample='rslice') # returns a GRB prompt result object result.plot_lightcurve(random_models=1000) diff --git a/examples/supernova_example.py b/examples/supernova_example.py index e87d254ef..00e6706db 100644 --- a/examples/supernova_example.py +++ b/examples/supernova_example.py @@ -12,17 +12,17 @@ data = redback.get_data.get_supernova_data_from_open_transient_catalog_data(sne) supernova = redback.supernova.Supernova.from_open_access_catalogue(name=sne, data_mode='flux_density') -fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(12, 8)) bands = ["g", "i"] -supernova.plot_multiband(figure=fig, axes=axes, filters=["J", "H", "g", "i"]) +supernova.plot_multiband(filters=["J", "H", "g", "i"]) # use default priors priors = redback.priors.get_priors(model=model) priors['redshift'] = 0.677 -model_kwargs = dict(frequency=redback.utils.bands_to_frequency(bands), output_format='flux_density') +model_kwargs = dict(frequency=supernova.filtered_frequencies, output_format='flux_density') # returns a supernova result object -result = redback.fit_model(name=sne, transient=supernova, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, sample='rslice', nlive=200, resume=False) +result = redback.fit_model(transient=supernova, model=model, sampler=sampler, model_kwargs=model_kwargs, + prior=priors, sample='rslice', nlive=200, resume=True) result.plot_corner() -result.plot_lightcurve(random_models=1000) +result.plot_lightcurve(random_models=100) +result.plot_multiband_lightcurve() diff --git a/examples/tde_example.py b/examples/tde_example.py index d3a8a6673..80cd87732 100644 --- a/examples/tde_example.py +++ b/examples/tde_example.py @@ -8,11 +8,10 @@ data = redback.get_data.get_tidal_disruption_event_data_from_open_transient_catalog_data(tde) -tidal_disruption_event = redback.tde.TDE.from_open_access_catalogue(tidal_disruption_event, data_mode='flux_density') +tidal_disruption_event = redback.tde.TDE.from_open_access_catalogue(tde, data_mode='flux_density') # lets make a plot of the data -fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(12, 8)) -tidal_disruption_event.plot_multiband(figure=fig, axes=axes, filters=["J", "H", "g", "i"]) +tidal_disruption_event.plot_multiband(filters=["V", "r", "g", "i"]) # we are now only going to fit u and r bands. We set this using the transient object and the active bands attribute. # By default all data is used, this is just for speed in the example or if you only trust some of the data/physics. @@ -23,14 +22,14 @@ priors = redback.priors.get_priors(model=model) # we know the redshift for PS18kh so we just fix the prior for the redshift to the known value. # We can do this through the metadata that was downloaded alongside the data, or if you just know it. -priors['redshift'] = tidal_disruption_event.redshift +# priors['redshift'] = tidal_disruption_event.redshift -model_kwargs = dict(frequency=redback.utils.bands_to_frequency(bands), output_format='flux_density') +model_kwargs = dict(frequency=tidal_disruption_event.filtered_frequencies, output_format='flux_density') # returns a tde result object -result = redback.fit_model(name=tde, transient=tidal_disruption_event, model=model, sampler=sampler, +result = redback.fit_model(transient=tidal_disruption_event, model=model, sampler=sampler, model_kwargs=model_kwargs, prior=priors, sample='rslice', nlive=200, resume=False) result.plot_corner() -result.plot_lightcurve(random_models=1000) +result.plot_lightcurve(random_models=100) diff --git a/redback/plotting.py b/redback/plotting.py index dd72493d2..48247dbe2 100644 --- a/redback/plotting.py +++ b/redback/plotting.py @@ -38,7 +38,10 @@ def _get_times(self, axes: matplotlib.axes.Axes) -> np.ndarray: @property def xlim_low(self): - return 0.5 * self.transient.x[0] + xlim_low = 0.5 * self.transient.x[0] + if xlim_low == 0: + xlim_low += 1e-3 + return xlim_low @property def xlim_high(self): diff --git a/redback/priors/arnett.prior b/redback/priors/arnett.prior index 1dd25117a..86a351681 100644 --- a/redback/priors/arnett.prior +++ b/redback/priors/arnett.prior @@ -4,4 +4,4 @@ mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$ vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') kappa_gamma = Uniform(1e-4, 1e4, 'kappa_gamma', latex_label = r'$\kappa_{\gamma}$') -temperature_floor = LogUniform(1e3,1e5,name = 'temperature_floor', latex_label = r'$T_{\mathrm{floor}$ [k]}') \ No newline at end of file +temperature_floor = LogUniform(1e3,1e5,name = 'temperature_floor', latex_label = r'$T_{\mathrm{floor}}$ [k]') \ No newline at end of file diff --git a/redback/priors/tde_analytical.prior b/redback/priors/tde_analytical.prior new file mode 100644 index 000000000..362f91575 --- /dev/null +++ b/redback/priors/tde_analytical.prior @@ -0,0 +1,8 @@ +mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$') +vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') +kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') +kappa_gamma = Uniform(1e-4, 1e4, 'kappa_gamma', latex_label = r'$\kappa_{\gamma}$') +temperature_floor = LogUniform(1e3,1e5,name = 'temperature_floor', latex_label = r'$T_{\mathrm{floor}$ [k]}') +redshift = Uniform(0.01, 3, 'redshift', latex_label=r'$z$') +l0 = LogUniform(1e51, 1e58, "l0", latex_label="$l_0$") +t_0 = LogUniform(1e-4, 1e2, "t_0", latex_label="$t_0$") \ No newline at end of file diff --git a/redback/sampler.py b/redback/sampler.py index 609d3842d..459e77a8c 100644 --- a/redback/sampler.py +++ b/redback/sampler.py @@ -18,7 +18,7 @@ dirname = os.path.dirname(__file__) -def fit_model(name, transient, model, outdir=None, sampler='dynesty', nlive=2000, prior=None, +def fit_model(transient, model, outdir=None, sampler='dynesty', nlive=2000, prior=None, walks=200, truncate=True, use_photon_index_prior=False, truncate_method='prompt_time_error', resume=True, save_format='json', model_kwargs=None, **kwargs): """ @@ -61,17 +61,17 @@ def fit_model(name, transient, model, outdir=None, sampler='dynesty', nlive=2000 truncate=truncate, use_photon_index_prior=use_photon_index_prior, truncate_method=truncate_method, **kwargs) elif isinstance(transient, PromptTimeSeries): - return _fit_prompt(name=name, transient=transient, model=model, outdir=outdir, sampler=sampler, nlive=nlive, + return _fit_prompt(transient=transient, model=model, outdir=outdir, sampler=sampler, nlive=nlive, prior=prior, walks=walks, use_photon_index_prior=use_photon_index_prior, resume=resume, save_format=save_format, model_kwargs=model_kwargs, **kwargs) elif isinstance(transient, Supernova): - return _fit_supernova(name=name, transient=transient, model=model, outdir=outdir, sampler=sampler, nlive=nlive, + return _fit_kilonova(transient=transient, model=model, outdir=outdir, sampler=sampler, nlive=nlive, prior=prior, walks=walks, truncate=truncate, use_photon_index_prior=use_photon_index_prior, truncate_method=truncate_method, resume=resume, save_format=save_format, model_kwargs=model_kwargs, **kwargs) elif isinstance(transient, TDE): - return _fit_tde(name=name, transient=transient, model=model, outdir=outdir, sampler=sampler, nlive=nlive, + return _fit_kilonova(transient=transient, model=model, outdir=outdir, sampler=sampler, nlive=nlive, prior=prior, walks=walks, truncate=truncate, use_photon_index_prior=use_photon_index_prior, truncate_method=truncate_method, resume=resume, save_format=save_format, model_kwargs=model_kwargs, **kwargs) @@ -147,7 +147,7 @@ def _fit_kilonova(transient, model, outdir=None, sampler='dynesty', nlive=3000, return result -def _fit_prompt(name, transient, model, outdir, integrated_rate_function=True, sampler='dynesty', nlive=3000, +def _fit_prompt(transient, model, outdir, integrated_rate_function=True, sampler='dynesty', nlive=3000, prior=None, walks=1000, use_photon_index_prior=False, resume=True, save_format='json', model_kwargs=None, **kwargs): diff --git a/redback/sed.py b/redback/sed.py index 3245454c7..7c45f37b3 100644 --- a/redback/sed.py +++ b/redback/sed.py @@ -62,7 +62,7 @@ def calculate_flux_density(self): mask = wavelength < cutoff_wavelength sed = np.zeros(len(self.time)) - sed[mask] = flux_const * (self.r_photosphere[mask]**2 / cutoff_wavelength[mask] / wavelength ** 4) \ + sed[mask] = flux_const * (self.r_photosphere[mask]**2 / cutoff_wavelength / wavelength[mask] ** 4) \ / np.expm1(x_const / wavelength[mask] / self.temperature[mask]) sed[~mask] = flux_const * (self.r_photosphere[~mask]**2 / wavelength[~mask]**5) \ / np.expm1(x_const / wavelength[~mask] / self.temperature[~mask]) @@ -82,10 +82,13 @@ def calculate_flux_density(self): tp3 = tp**3 nxcs = x_const * np.array(range(1, 11)) - f_blue_reds = np.sum((np.exp(-nxcs / (cutoff_wavelength * tp)) * ( - nxcs ** 2 + 2 * (nxcs * cutoff_wavelength * tp + cutoff_wavelength**2 * tp2)) / (nxcs ** 3 * cutoff_wavelength**3)) + ( - (6 * tp3 - np.exp(-nxcs / (cutoff_wavelength * tp)) * (nxcs ** 3 + 3 * nxcs ** 2 * cutoff_wavelength * tp + 6 * ( - nxcs * cutoff_wavelength**2 * tp2 + cutoff_wavelength**3 *tp3)) / cutoff_wavelength**3) / (nxcs ** 4), 1)) + f_blue_reds = \ + np.sum( + np.exp(-nxcs / (cutoff_wavelength * tp)) * (nxcs ** 2 + 2 * (nxcs * cutoff_wavelength * tp + cutoff_wavelength**2 * tp2)) / (nxcs ** 3 * cutoff_wavelength**3) + + ( + (6 * tp3 - np.exp(-nxcs / (cutoff_wavelength * tp)) * (nxcs ** 3 + 3 * nxcs ** 2 * cutoff_wavelength * tp + 6 * (nxcs * cutoff_wavelength**2 * tp2 + cutoff_wavelength**3 *tp3)) / cutoff_wavelength**3) / (nxcs ** 4) + ), 1 + ) norms /= f_blue_reds @@ -95,7 +98,7 @@ def calculate_flux_density(self): self.sed = sed # sed units are erg/s/Angstrom - need to turn them into flux density compatible units - units = uu.erg / uu.s / uu.hz / uu.cm**2. + units = uu.erg / uu.s / uu.Hz / uu.cm**2. sed = sed / (4 * np.pi * self.luminosity_distance ** 2) * lambda_to_nu(1.) # add units diff --git a/redback/transient/transient.py b/redback/transient/transient.py index 10bd98d17..7e900c63b 100644 --- a/redback/transient/transient.py +++ b/redback/transient/transient.py @@ -763,8 +763,12 @@ def from_open_access_catalogue( ------- OpticalTransient: A class instance. """ - directory_structure = redback.get_data.directory.open_access_directory_structure(transient=name, - transient_type=cls.__name__.lower()) + if cls.__name__ == "TDE": + transient_type = "tidal_disruption_event" + else: + transient_type = cls.__name__.lower() + directory_structure = redback.get_data.directory.open_access_directory_structure( + transient=name, transient_type=transient_type) time_days, time_mjd, flux_density, flux_density_err, magnitude, magnitude_err, bands, system = \ cls.load_data(processed_file_path=directory_structure.processed_file_path, data_mode="all") return cls(name=name, data_mode=data_mode, time=time_days, time_err=None, time_mjd=time_mjd, diff --git a/redback/transient_models/extinction_models.py b/redback/transient_models/extinction_models.py index 25cdcf699..0224d1606 100644 --- a/redback/transient_models/extinction_models.py +++ b/redback/transient_models/extinction_models.py @@ -1,6 +1,7 @@ from inspect import isfunction import numpy as np +import redback.utils from redback.transient_models.fireball_models import predeceleration from redback.utils import logger, calc_ABmag_from_flux_density, citation_wrapper import astropy.units as uu @@ -67,8 +68,8 @@ def _perform_extinction(flux_density, frequency, av, r_v): """ import extinction # noqa # convert to angstrom - frequency = (frequency * uu.Hz).to(uu.Angstrom).value - mag_extinction = extinction.fitzpatrick99(frequency, av, r_v=r_v) + angstroms = redback.utils.nu_to_lambda(frequency) + mag_extinction = extinction.fitzpatrick99(angstroms, av, r_v=r_v) flux_density = extinction.apply(mag_extinction, flux_density) return flux_density @@ -92,7 +93,7 @@ def _evaluate_extinction_model(time, av, model_type, **kwargs): r_v = kwargs.get('r_v', 3.1) flux_density = _perform_extinction(flux_density=flux_density, frequency=frequency, av=av, r_v=r_v) if kwargs['output_format'] == 'flux_density': - return flux_density.value + return flux_density elif kwargs['output_format'] == 'magnitude': return calc_ABmag_from_flux_density(flux_density).value diff --git a/redback/transient_models/supernova_models.py b/redback/transient_models/supernova_models.py index fd48c1f77..77954d645 100644 --- a/redback/transient_models/supernova_models.py +++ b/redback/transient_models/supernova_models.py @@ -122,9 +122,7 @@ def arnett_bolometric(time, f_nickel, mej, interaction_process=ip.Diffusion, **k return lbol @citation_wrapper('https://ui.adsabs.harvard.edu/abs/1982ApJ...253..785A/abstract') -def arnett(time, redshift, f_nickel, mej, interaction_process=ip.Diffusion, - photosphere=photosphere.TemperatureFloor, - sed=sed.Blackbody, **kwargs): +def arnett(time, redshift, f_nickel, mej, **kwargs): """ :param time: time in days :param redshift: source redshift @@ -139,13 +137,18 @@ def arnett(time, redshift, f_nickel, mej, interaction_process=ip.Diffusion, e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor :return: flux_density or magnitude depending on output_format kwarg """ + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + _photosphere = kwargs.get("photosphere", photosphere.TemperatureFloor) + _sed = kwargs.get("sed", sed.Blackbody) + + frequency = kwargs['frequency'] frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value - lbol = arnett_bolometric(time=time, f_nickel=f_nickel, mej=mej, interaction_process=interaction_process, **kwargs) - photo = photosphere(time=time, luminosity=lbol, **kwargs) - sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + lbol = arnett_bolometric(time=time, f_nickel=f_nickel, mej=mej, interaction_process=_interaction_process, **kwargs) + photo = _photosphere(time=time, luminosity=lbol, **kwargs) + sed_1 = _sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, frequency=frequency, luminosity_distance=dl) flux_density = sed_1.flux_density diff --git a/redback/transient_models/tde_models.py b/redback/transient_models/tde_models.py index 35bf5ab31..f039aecdf 100644 --- a/redback/transient_models/tde_models.py +++ b/redback/transient_models/tde_models.py @@ -15,7 +15,7 @@ def _analytic_fallback(time, l0, t_0): """ mask = time - t_0 > 0. lbol = np.zeros(len(time)) - lbol[mask] = l0 / (time * 86400)**(5./3.) + lbol[mask] = l0 / (time[mask] * 86400)**(5./3.) lbol[~mask] = l0 / (t_0 * 86400)**(5./3.) return lbol @@ -43,10 +43,7 @@ def tde_analytical_bolometric(time, l0, t_0,interaction_process = ip.Diffusion, return lbol @citation_wrapper('redback') -def tde_analytical(time, redshift, l0, t_0, interaction_process=ip.Diffusion, - photosphere=photosphere.TemperatureFloor, - sed=sed.CutoffBlackbody, - **kwargs): +def tde_analytical(time, redshift, l0, t_0, **kwargs): """ :param time: rest frame time in days :param l0: bolometric luminosity at 1 second in cgs @@ -59,15 +56,19 @@ def tde_analytical(time, redshift, l0, t_0, interaction_process=ip.Diffusion, e.g., for Diffusion TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor :return: flux_density or magnitude depending on output_format kwarg """ + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + _photosphere = kwargs.get("photosphere", photosphere.TemperatureFloor) + _sed = kwargs.get("sed", sed.CutoffBlackbody) + frequency = kwargs['frequency'] cutoff_wavelength = kwargs.get('cutoff_wavelength', 3000) frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value - lbol = tde_analytical_bolometric(time=time, l0=l0, t_0=t_0, interaction_process=interaction_process) + lbol = tde_analytical_bolometric(time=time, l0=l0, t_0=t_0, interaction_process=_interaction_process, **kwargs) - photo = photosphere(time=time, luminosity=lbol, **kwargs) - sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, - frequency=frequency, luminosity_distance=dl, cutoff_wavelength=cutoff_wavelength) + photo = _photosphere(time=time, luminosity=lbol, **kwargs) + sed_1 = _sed(time=time, temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + frequency=frequency, luminosity_distance=dl, cutoff_wavelength=cutoff_wavelength, luminosity=lbol) flux_density = sed_1.flux_density From b5a74523c191956969b6b02e9eaa3f9665993c1f Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Fri, 11 Mar 2022 13:47:44 +0100 Subject: [PATCH 27/37] some docs and general fixes to models --- redback/transient_models/fireball_models.py | 14 +- .../integrated_flux_afterglow_models.py | 19 +- redback/transient_models/kilonova_models.py | 20 +- .../magnetar_boosted_ejecta_models.py | 8 +- redback/transient_models/magnetar_models.py | 268 ++++++++++-------- 5 files changed, 192 insertions(+), 137 deletions(-) diff --git a/redback/transient_models/fireball_models.py b/redback/transient_models/fireball_models.py index 6605b9b7a..59de84db0 100644 --- a/redback/transient_models/fireball_models.py +++ b/redback/transient_models/fireball_models.py @@ -1,25 +1,25 @@ -import numpy as np - +from redback.utils import citation_wrapper +@citation_wrapper('redback') def predeceleration(time, aa, mm, t0, **kwargs): """ - :param time: - :param aa: + :param time: time array in seconds + :param aa: amplitude term for powerlaw :param mm: deceleration powerlaw gradient; typically 3 but depends on physics :param t0: time GRB went off. - :param kwargs: + :param kwargs: None :return: deceleration powerlaw; units are arbitrary and dependent on a_1. """ return aa * (time - t0)**mm - +@citation_wrapper('redback') def one_component_fireball_model(time, a_1, alpha_1, **kwargs): """ :param time: time array for power law :param a_1: power law decay amplitude :param alpha_1: power law decay exponent :param kwargs: - :return: + :return: powerlaw; units are arbitrary and dependent on a_1. """ return a_1 * time ** alpha_1 diff --git a/redback/transient_models/integrated_flux_afterglow_models.py b/redback/transient_models/integrated_flux_afterglow_models.py index 5deb2d59b..3d870ee08 100644 --- a/redback/transient_models/integrated_flux_afterglow_models.py +++ b/redback/transient_models/integrated_flux_afterglow_models.py @@ -13,6 +13,13 @@ @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2021arXiv210510108S/abstract') def integrated_flux_afterglowpy_base_model(time, **kwargs): + """ + Synchrotron afterglow with integrated flux + + :param time: time in days + :param kwargs:all kwargs required by model + frequency: a list of two frequencies to integrate over. + :return: integrated flux + """ from ..model_library import modules_dict # import model library in function to avoid circular dependency base_model = kwargs['base_model'] @@ -42,11 +49,15 @@ def integrated_flux_afterglowpy_base_model(time, **kwargs): @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2021arXiv210510108S/abstract') def integrated_flux_rate_model(time, **kwargs): """ - :param time: - :param kwargs: - :return: rate + Synchrotron afterglow with approximate calculation of the counts + + :param time: time in days + :param kwargs:all kwargs required by model + frequency: an array of two frequencies to integrate over. + + prefactor an array of values same size as time array + or float which calculates the effective Ei/area for the specific time bin. + :return: counts """ - prefactor = kwargs.get('prefactor', 0) + prefactor = kwargs.get('prefactor', 1) dt = kwargs.get('dt', 1) background_rate = kwargs.get('background_rate', 0) integrated_flux = integrated_flux_afterglowpy_base_model(time, **kwargs) diff --git a/redback/transient_models/kilonova_models.py b/redback/transient_models/kilonova_models.py index b2c5a62c9..fea364d5f 100644 --- a/redback/transient_models/kilonova_models.py +++ b/redback/transient_models/kilonova_models.py @@ -49,8 +49,8 @@ def two_layer_stratified_kilonova(time, redshift, mass, vej_1, vej_2, kappa, bet :return: flux_density or magnitude """ velocity_array = np.array([vej_1, vej_2]) - flux_density = _kilonova_hr(time, redshift, mass, velocity_array, kappa, beta, **kwargs) - return flux_density + output = _kilonova_hr(time, redshift, mass, velocity_array, kappa, beta, **kwargs) + return output def _kilonova_hr(time, redshift, mass, velocity_array, kappa_array, beta, **kwargs): @@ -70,7 +70,7 @@ def _kilonova_hr(time, redshift, mass, velocity_array, kappa_array, beta, **kwar """ frequency = kwargs['frequency'] # convert to source frame time and frequency - time = time * 86400 + time = time * day_to_s frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value _, temperature, r_photosphere = _kilonova_hr_sourceframe(time, mass, velocity_array, kappa_array, beta) @@ -140,7 +140,7 @@ def three_component_kilonova_model(time, redshift, mej_1, vej_1, temperature_flo """ frequency = kwargs['frequency'] # convert to source frame time and frequency - time = time * 86400 + time = time * day_to_s frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value @@ -190,7 +190,7 @@ def two_component_kilonova_model(time, redshift, mej_1, vej_1, temperature_floor """ frequency = kwargs['frequency'] # convert to source frame time and frequency - time = time * 86400 + time = time * day_to_s frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value @@ -257,7 +257,7 @@ def one_component_kilonova_model(time, redshift, mej, vej, kappa, **kwargs): output_format :return: flux_density or magnitude """ - time = time * 86400 + time = time * day_to_s frequency = kwargs['frequency'] time_temp = np.geomspace(1e-4, 1e7, 300) _, temperature, r_photosphere = _one_component_kilonova_model(time_temp, mej, vej, kappa, **kwargs) @@ -290,7 +290,7 @@ def _one_component_kilonova_model(time, mej, vej, kappa, **kwargs): :param kwargs: temperature_floor :return: bolometric_luminosity, temperature, r_photosphere """ - tdays = time/86400 + tdays = time/day_to_s # set up kilonova physics av, bv, dv = interpolated_barnes_and_kasen_thermalisation_efficiency(mej, vej) @@ -335,6 +335,7 @@ def metzger_kilonova_model(time, redshift, mej, vej, beta, kappa_r, **kwargs): frequency (frequency to calculate - Must be same length as time array or a single number) :return: flux_density or magnitude """ + time = time * day_to_s frequency = kwargs['frequency'] time_temp = np.geomspace(1e-4, 1e7, 300) bolometric_luminosity, temperature, r_photosphere = _metzger_kilonova_model(time_temp, mej, vej, beta, @@ -372,7 +373,7 @@ def _metzger_kilonova_model(time, mej, vej, beta, kappa_r, **kwargs): neutron_precursor_switch = kwargs.get('neutron_precursor_switch', True) time = time - tdays = time/86400 + tdays = time/day_to_s time_len = len(time) mass_len = 250 @@ -506,7 +507,7 @@ def _generate_single_lightcurve_at_times(model, times, **parameters): model: str The `gwemlightcurve` model, e.g. 'DiUj2017' times: array_like - Times at which we interpolate the `gwemlightcurves` values + Times at which we interpolate the `gwemlightcurves` values, in days parameters: dict Function parameters for the given model. Returns @@ -530,7 +531,6 @@ def _generate_single_lightcurve_at_times(model, times, **parameters): def _gwemlightcurve_interface_factory(model): """ Generates `bilby`-compatible functions from `gwemlightcurve` models. - This is currently very inefficient as there is an interpolation step. Parameters ---------- diff --git a/redback/transient_models/magnetar_boosted_ejecta_models.py b/redback/transient_models/magnetar_boosted_ejecta_models.py index 7e08fb0c4..fc59e3433 100644 --- a/redback/transient_models/magnetar_boosted_ejecta_models.py +++ b/redback/transient_models/magnetar_boosted_ejecta_models.py @@ -38,7 +38,7 @@ def metzger_magnetar_boosted_kilonova_model(time, redshift, mej, vej, beta, kapp temp_func = interp1d(time_temp, y=temperature) rad_func = interp1d(time_temp, y=r_photosphere) # convert to source frame time and frequency - time = time * 86400 + time = time * day_to_s frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) temp = temp_func(time) @@ -73,7 +73,7 @@ def _metzger_magnetar_boosted_kilonova_model(time, mej, vej, beta, kappa_r, l0, time = time - tdays = time/86400 + tdays = time/day_to_s time_len = len(time) mass_len = 250 @@ -197,7 +197,7 @@ def _metzger_magnetar_boosted_kilonova_model(time, mej, vej, beta, kappa_r, l0, ejecta_albedo = kwargs.get('ejecta_albedo', 0.5) pair_cascade_fraction = kwargs.get('pair_cascade_fraction', 0.01) tlife_t = (0.6/(1 - ejecta_albedo))*(pair_cascade_fraction/0.1)**0.5 * (lsd/1.0e45)**0.5 \ - * (v0/(0.3*speed_of_light))**(0.5) * (time/86400)**(-0.5) + * (v0/(0.3*speed_of_light))**(0.5) * (time/day_to_s)**(-0.5) bolometric_luminosity = bolometric_luminosity / (1.0 + tlife_t) temperature = (bolometric_luminosity / (4.0 * np.pi * (r_photosphere) ** (2.0) * sigma_sb)) ** (0.25) @@ -365,7 +365,7 @@ def mergernova(time, redshift, mej, beta, ejecta_radius, kappa, n_ism, l0, tau_s rad_func = interp1d(time_temp, y=radius) d_func = interp1d(time_temp, y=doppler_factor) # convert to source frame time and frequency - time = time * 86400 + time = time * day_to_s frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) temp = temp_func(time) diff --git a/redback/transient_models/magnetar_models.py b/redback/transient_models/magnetar_models.py index a2720956b..aff5c6450 100644 --- a/redback/transient_models/magnetar_models.py +++ b/redback/transient_models/magnetar_models.py @@ -27,7 +27,13 @@ def _evolving_magnetar_only(time, mu0, muinf, p0, sinalpha0, tm, II, **kwargs): """ Millisecond magnetar model with evolution of inclination angle - :param time: time + :param time: time in seconds + :param mu0: Initial magnetic moment [10^33 G cm^3] + :param muinf: Magnetic moment when field relaxes [10^33 G cm^3] + :param p0: initial spin period + :param sinalpha0: initial sin(alpha0) where alpha is the angle between B and P axes. + :param tm: magnetic field decay timescale in days + :param II: Moment of inertia in cgs :param kwargs: key word argument for handling plotting :return: luminosity (depending on scaling) as a function of time. """ @@ -53,14 +59,21 @@ def _evolving_magnetar_only(time, mu0, muinf, p0, sinalpha0, tm, II, **kwargs): return luminosity / 1e50 @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2019ApJ...886....5S/abstract') -def evolving_magnetar(time, A_1, alpha_1, mu0, muinf, p0, sinalpha0, tm, II, **kwargs): +def evolving_magnetar(time, a_1, alpha_1, mu0, muinf, p0, sinalpha0, tm, II, **kwargs): """ - :param time: time - :param L0: luminosity parameter - :param tau: spin-down damping timescale - :param nn: braking index + Millisecond magnetar model with evolution of inclination angle + + :param time: time in seconds + :param A_1: amplitude of curvature effect power law + :param alpha_1: index of curvature effect power law + :param mu0: Initial magnetic moment [10^33 G cm^3] + :param muinf: Magnetic moment when field relaxes [10^33 G cm^3] + :param p0: initial spin period + :param sinalpha0: initial sin(alpha0) where alpha is the angle between B and P axes. + :param tm: magnetic field decay timescale in days + :param II: Moment of inertia in cgs :param kwargs: key word argument for handling plotting - :return: luminosity or flux (depending on scaling) as a function of time. + :return: luminosity (depending on scaling) as a function of time. """ pl = one_component_fireball_model(time=time, a_1=A_1, alpha_1=alpha_1) magnetar = _evolving_magnetar_only(time=time, mu0=mu0, muinf=muinf, @@ -70,12 +83,12 @@ def evolving_magnetar(time, A_1, alpha_1, mu0, muinf, p0, sinalpha0, tm, II, **k @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2017ApJ...843L...1L/abstract') def _magnetar_only(time, l0, tau, nn, **kwargs): """ - :param time: time - :param l0: luminosity parameter + :param time: time in seconds + :param l0: initial luminosity parameter :param tau: spin-down damping timescale :param nn: braking index - :param kwargs: key word argument for handling plotting - :return: luminosity or flux (depending on scaling) as a function of time. + :param kwargs: key word arguments for handling plotting/other functionality + :return: luminosity or flux (depending on scaling of l0) as a function of time. """ lum = l0 * (1. + time / tau) ** ((1. + nn) / (1. - nn)) return lum @@ -83,15 +96,15 @@ def _magnetar_only(time, l0, tau, nn, **kwargs): @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018PhRvD..98d3011S/abstract') def gw_magnetar(time, a_1, alpha_1, fgw0, tau, nn, log_ii, **kwargs): """ - :param time: - :param a_1: - :param alpha_1: + :param time: time in seconds + :param A_1: amplitude of curvature effect power law + :param alpha_1: index of curvature effect power law :param fgw0: initial gravitational-wave frequency - :param tau: - :param nn: + :param tau: spin-down damping timescale + :param nn: braking index :param log_ii: log10 moment of inertia :param eta: fixed to 0.1, its a fudge factor for the efficiency - :param kwargs: + :param kwargs: key word arguments for handling plotting/other functionality :return: luminosity """ eta = 0.1 @@ -108,14 +121,16 @@ def gw_magnetar(time, a_1, alpha_1, fgw0, tau, nn, log_ii, **kwargs): @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2017ApJ...843L...1L/abstract') def full_magnetar(time, a_1, alpha_1, l0, tau, nn, **kwargs): """ - :param time: - :param a_1: - :param alpha_1: - :param l0: - :param tau: - :param nn: - :param kwargs: - :return: + Generalised millisecond magnetar with curvature effect power law + + :param time: time in seconds + :param A_1: amplitude of curvature effect power law + :param alpha_1: index of curvature effect power law + :param l0: initial luminosity parameter + :param tau: spin-down damping timescale + :param nn: braking index + :param kwargs: key word arguments for handling plotting/other functionality + :return: luminosity or flux (depending on scaling of l0) as a function of time. """ pl = one_component_fireball_model(time=time, a_1=a_1, alpha_1=alpha_1) mag = _magnetar_only(time=time, l0=l0, tau=tau, nn=nn) @@ -124,15 +139,17 @@ def full_magnetar(time, a_1, alpha_1, l0, tau, nn, **kwargs): @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020PhRvD.101f3021S/abstract') def collapsing_magnetar(time, a_1, alpha_1, l0, tau, nn, tcol, **kwargs): """ - :param time: - :param a_1: - :param alpha_1: - :param l0: - :param tau: - :param nn: - :param tcol: - :param kwargs: - :return: + Generalised millisecond magnetar with curvature effect power law and a collapse time + + :param time: time in seconds + :param A_1: amplitude of curvature effect power law + :param alpha_1: index of curvature effect power law + :param l0: initial luminosity parameter + :param tau: spin-down damping timescale + :param nn: braking index + :param tcol: collapse time in seconds + :param kwargs: key word arguments for handling plotting/other functionality + :return: luminosity or flux (depending on scaling of l0) as a function of time. """ pl = one_component_fireball_model(time, a_1, alpha_1) mag = np.heaviside(tcol - time, 1e-50) * _magnetar_only(time, l0, tau, nn) @@ -151,6 +168,8 @@ def general_magnetar(time, a_1, alpha_1, :param delta_time_one: time between start and end of prompt emission :param alpha_2: Reparameterized braking index n :param delta_time_two: time between end of prompt emission and end of magnetar model plateau phase, (tau) + :param kwargs: key word arguments for handling plotting/other functionality + :return: luminosity or flux (depending on scaling of l0) as a function of time. """ time_one = delta_time_one @@ -172,13 +191,14 @@ def general_magnetar(time, a_1, alpha_1, def _integral_general(time, t0, kappa, tau, nn, **kwargs): """ General integral for radiative losses model - :param time: - :param t0: - :param kappa: - :param tau: - :param nn: - :param kwargs: - :return: + + :param time: time in seconds + :param t0: time for radiative losses to start in seconds + :param kappa: radiative efficiency + :param tau: spin down damping timescale in seconds + :param nn: braking index + :param kwargs: key word arguments for handling plotting/other functionality + :return: luminosity """ first_term, second_term = _get_integral_terms(time=time, t0=t0, kappa=kappa, tau=tau, nn=nn) return first_term - second_term @@ -186,21 +206,32 @@ def _integral_general(time, t0, kappa, tau, nn, **kwargs): @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020MNRAS.499.5986S/abstract') def _integral_general_collapsing(time, t0, kappa, tau, nn, tcol, **kwargs): """ - General collapsing integral for radiative losses model - :param time: - :param t0: - :param kappa: - :param tau: - :param nn: - :param tcol: - :param kwargs: - :return: + General integral for radiative losses model + + :param time: time in seconds + :param t0: time for radiative losses to start in seconds + :param kappa: radiative efficiency + :param tau: spin down damping timescale in seconds + :param nn: braking index + :param tcol: collapse time in seconds + :param kwargs: key word arguments for handling plotting/other functionality + :return: luminosity """ first_term, second_term = _get_integral_terms(time=time, t0=t0, kappa=kappa, tau=tau, nn=nn) return np.heaviside(tcol - time, 1e-50) * (first_term - second_term) @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020MNRAS.499.5986S/abstract') def _get_integral_terms(time, t0, kappa, tau, nn): + """ + Get integral terms for radiative losses model + + :param time: time in seconds + :param t0: time for radiative losses to start in seconds + :param kappa: radiative efficiency + :param tau: spin down damping timescale in seconds + :param nn: braking index + :return: first and second terms + """ alpha = (1 + nn) / (-1 + nn) pft = ss.hyp2f1(1 + kappa, alpha, 2 + kappa, -time / tau) pst = ss.hyp2f1(1 + kappa, alpha, 2 + kappa, -t0 / tau) @@ -210,6 +241,15 @@ def _get_integral_terms(time, t0, kappa, tau, nn): @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020MNRAS.499.5986S/abstract') def _integral_mdr(time, t0, kappa, a, **kwargs): + """ + Calculate integral for vacuum dipole radiation + + :param time: time in seconds + :param t0: time for radiative losses to start in seconds + :param kappa: radiative efficiency + :param a: 1/tau (spin down damping timescale) + :return: integral + """ z_f = (1 + a * time) ** (-1) z_int = (1 + a * t0) ** (-1) divisor_i = a ** (1 + kappa) * (kappa - 1) * (1 + a * t0) ** (1 - kappa) @@ -223,16 +263,16 @@ def piecewise_radiative_losses(time, a_1, alpha_1, l0, tau, nn, kappa, t0, **kwa """ Assumes smoothness and continuity between the prompt and magnetar term by fixing e0 variable - :param time: - :param a_1: - :param alpha_1: - :param l0: - :param tau: - :param nn: - :param kappa: - :param t0: - :param kwargs: - :return: + :param time: time in seconds + :param A_1: amplitude of curvature effect power law + :param alpha_1: index of curvature effect power law + :param l0: initial luminosity parameter + :param tau: spin-down damping timescale + :param nn: braking index + :param kappa: radiative efficiency + :param t0: time for radiative losses to start in seconds + :param kwargs: key word arguments for handling plotting/other functionality + :return: luminosity """ pl_time = np.where(time <= t0) magnetar_time = np.where(time > t0) @@ -252,17 +292,17 @@ def radiative_losses(time, a_1, alpha_1, l0, tau, nn, kappa, t0, log_e0, **kwarg """ radiative losses model with a step function, indicating the magnetar term turns on at T0 - :param time: - :param a_1: - :param alpha_1: - :param l0: - :param tau: - :param nn: - :param kappa: - :param t0: - :param log_e0: - :param kwargs: - :return: + :param time: time in seconds + :param A_1: amplitude of curvature effect power law + :param alpha_1: index of curvature effect power law + :param l0: initial luminosity parameter + :param tau: spin-down damping timescale + :param nn: braking index + :param kappa: radiative efficiency + :param t0: time for radiative losses to start in seconds + :param log_e0: log10 E0 to connect curvature effect energy with transition point energy, captures flares + :param kwargs: key word arguments for handling plotting/other functionality + :return: luminosity """ e0 = 10 ** log_e0 pl = one_component_fireball_model(time, a_1, alpha_1) @@ -279,14 +319,16 @@ def radiative_only(time, l0, tau, nn, kappa, t0, log_e0, **kwargs): """ radiative losses model only - :param time: - :param l0: - :param tau: - :param nn: - :param kappa: - :param t0: - :param log_e0: - :param kwargs: + :param time: time in seconds + :param l0: initial luminosity parameter + :param tau: spin-down damping timescale + :param nn: braking index + :param kappa: radiative efficiency + :param t0: time for radiative losses to start in seconds + :param log_e0: log10 E0 to connect curvature effect energy smoothly with transition point energy + :param kwargs: key word arguments for handling plotting/other functionality + :return: luminosity + :return: """ e0 = 10 ** log_e0 @@ -303,17 +345,18 @@ def radiative_losses_smoothness(time, a_1, alpha_1, l0, tau, nn, kappa, t0, log_ """ radiative losses model with a step function, indicating the magnetar term turns on at T0 - :param time: - :param a_1: - :param alpha_1: - :param l0: - :param tau: - :param nn: - :param kappa: - :param t0: - :param log_e0: - :param kwargs: - :return: + :param time: time in seconds + :param A_1: amplitude of curvature effect power law + :param alpha_1: index of curvature effect power law + :param l0: initial luminosity parameter + :param tau: spin-down damping timescale + :param nn: braking index + :param kappa: radiative efficiency + :param t0: time for radiative losses to start in seconds + :param log_e0: log10 E0 to connect curvature effect energy smoothly with transition point energy + :param kwargs: key word arguments for handling plotting/other functionality + :return: luminosity + """ pl = one_component_fireball_model(time, a_1, alpha_1) e0 = 10 ** log_e0 @@ -333,16 +376,16 @@ def radiative_losses_mdr(time, a_1, alpha_1, l0, tau, kappa, log_e0, t0, **kwarg """ radiative losses model for vacuum dipole radiation - :param time: - :param a_1: - :param alpha_1: - :param l0: - :param tau: - :param kappa: - :param t0: - :param log_e0: - :param kwargs: - :return: + :param time: time in seconds + :param A_1: amplitude of curvature effect power law + :param alpha_1: index of curvature effect power law + :param l0: initial luminosity parameter + :param tau: spin-down damping timescale + :param kappa: radiative efficiency + :param t0: time for radiative losses to start in seconds + :param log_e0: log10 E0 to connect curvature effect energy smoothly with transition point energy + :param kwargs: key word arguments for handling plotting/other functionality + :return: luminosity """ a = 1. / tau e0 = 10 ** log_e0 @@ -360,18 +403,19 @@ def collapsing_radiative_losses(time, a_1, alpha_1, l0, tau, nn, tcol, kappa, t0 """ radiative losses model with collapse time - :param time: - :param a_1: - :param alpha_1: - :param l0: - :param tau: - :param nn: - :param tcol: - :param kappa: - :param t0: - :param log_e0: - :param kwargs: - :return: + :param time: time in seconds + :param A_1: amplitude of curvature effect power law + :param alpha_1: index of curvature effect power law + :param l0: initial luminosity parameter + :param tau: spin-down damping timescale + :param nn: braking index + :param tcol: collapse time in seconds + :param kappa: radiative efficiency + :param t0: time for radiative losses to start in seconds + :param log_e0: log10 E0 to connect curvature effect energy smoothly with transition point energy + :param kwargs: key word arguments for handling plotting/other functionality + :return: luminosity + """ e0 = 10 ** log_e0 pl = one_component_fireball_model(time, a_1, alpha_1) From c2c1d443246258fa57dc9c42a12eea577ffbc45c Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Fri, 11 Mar 2022 13:48:27 +0100 Subject: [PATCH 28/37] update prior to not have A_1 in caps --- redback/priors/evolving_magnetar.prior | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/redback/priors/evolving_magnetar.prior b/redback/priors/evolving_magnetar.prior index 700c30cdf..61a054271 100644 --- a/redback/priors/evolving_magnetar.prior +++ b/redback/priors/evolving_magnetar.prior @@ -1,4 +1,4 @@ -A_1 = LogUniform(1e-20, 1e20, 'a_1', latex_label = r'$A_{1}$') +a_1 = LogUniform(1e-20, 1e20, 'a_1', latex_label = r'$a_{1}$') alpha_1 = Uniform(-10, -0.5, 'alpha_1', latex_label = r'$\\alpha_{1}$') p0 = Uniform(0.7e-3, 0.1, 'p0', latex_label = r'$P_{0} [s]$') mu0 = Uniform(1e-3, 10, 'mu0', latex_label = r'$\mu_{0} [10^{33} G cm^{3}]$') From c85a38d63572cea12b5b4656f49a28885f145324 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Fri, 11 Mar 2022 14:16:57 +0100 Subject: [PATCH 29/37] Put kwargs.get for default dependency injections and rewrite some docs for sn models --- redback/transient_models/supernova_models.py | 250 +++++++++++-------- 1 file changed, 139 insertions(+), 111 deletions(-) diff --git a/redback/transient_models/supernova_models.py b/redback/transient_models/supernova_models.py index 77954d645..86de784ef 100644 --- a/redback/transient_models/supernova_models.py +++ b/redback/transient_models/supernova_models.py @@ -25,32 +25,29 @@ def thermal_synchrotron(): pass @citation_wrapper('redback') -def exponential_powerlaw_bolometric(time, lbol_0, alpha_1, alpha_2, tpeak_d, interaction_process = ip.Diffusion, - **kwargs): +def exponential_powerlaw_bolometric(time, lbol_0, alpha_1, alpha_2, tpeak_d, **kwargs): """ :param time: rest frame time in days :param lbol_0: bolometric luminosity scale in cgs :param alpha_1: first exponent :param alpha_2: second exponent :param tpeak_d: peak time in days - :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process - e.g., for Diffusion: kappa, kappa_gamma, mej (solar masses), vej (km/s), temperature_floor + e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor + :param interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :return: bolometric_luminosity """ + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) lbol = exponential_powerlaw(time, a_1=lbol_0, alpha_1=alpha_1, alpha_2=alpha_2, tpeak=tpeak_d, **kwargs) - if interaction_process is not None: - interaction_class = interaction_process(time=time, luminosity=lbol, **kwargs) + if _interaction_process is not None: + interaction_class = _interaction_process(time=time, luminosity=lbol, **kwargs) lbol = interaction_class.new_luminosity return lbol @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018ApJS..236....6G/abstract') -def sn_exponential_powerlaw(time, redshift, lbol_0, alpha_1, alpha_2, tpeak_d, - interaction_process = ip.Diffusion, - photosphere=photosphere.TemperatureFloor, - sed=sed.Blackbody,**kwargs): +def sn_exponential_powerlaw(time, redshift, lbol_0, alpha_1, alpha_2, tpeak_d, **kwargs): """ :param time: observer frame time in days :param redshift: source redshift @@ -58,24 +55,28 @@ def sn_exponential_powerlaw(time, redshift, lbol_0, alpha_1, alpha_2, tpeak_d, :param alpha_1: first exponent :param alpha_2: second exponent :param tpeak_d: peak time in days + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, mej (solar masses), vej (km/s), floor temperature :param interaction_process: Default is Diffusion. Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. - :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used - e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, mej (solar masses), vej (km/s), floor temperature :return: flux_density or magnitude depending on output_format kwarg """ + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + _photosphere = kwargs.get("photosphere", photosphere.TemperatureFloor) + _sed = kwargs.get("sed", sed.Blackbody) + frequency = kwargs['frequency'] frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value lbol = exponential_powerlaw_bolometric(time=time, lbol_0=lbol_0, alpha_1=alpha_1,alpha_2=alpha_2, tpeak_d=tpeak_d, - interaction_process=interaction_process, **kwargs) - photo = photosphere(time=time, luminosity=lbol, **kwargs) - sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + interaction_process=_interaction_process, **kwargs) + photo = _photosphere(time=time, luminosity=lbol, **kwargs) + sed_1 = _sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, frequency=frequency, luminosity_distance=dl) flux_density = sed_1.flux_density @@ -103,21 +104,21 @@ def _nickelcobalt_engine(time, f_nickel, mej, **kwargs): return lbol @citation_wrapper('https://ui.adsabs.harvard.edu/abs/1982ApJ...253..785A/abstract') -def arnett_bolometric(time, f_nickel, mej, interaction_process=ip.Diffusion, **kwargs): +def arnett_bolometric(time, f_nickel, mej, **kwargs): """ :param time: time in days :param f_nickel: fraction of nickel mass :param mej: total ejecta mass in solar masses - :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process - e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor + interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. + e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor :return: bolometric_luminosity """ - + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) lbol = _nickelcobalt_engine(time=time, f_nickel=f_nickel, mej=mej) - if interaction_process is not None: - interaction_class = interaction_process(time=time, luminosity=lbol, mej=mej, **kwargs) + if _interaction_process is not None: + interaction_class = _interaction_process(time=time, luminosity=lbol, mej=mej, **kwargs) lbol = interaction_class.new_luminosity return lbol @@ -128,13 +129,13 @@ def arnett(time, redshift, f_nickel, mej, **kwargs): :param redshift: source redshift :param f_nickel: fraction of nickel mass :param mej: total ejecta mass in solar masses + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor :param interaction_process: Default is Diffusion. Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. - :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used - e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor :return: flux_density or magnitude depending on output_format kwarg """ _interaction_process = kwargs.get("interaction_process", ip.Diffusion) @@ -174,28 +175,29 @@ def _basic_magnetar(time, p0, bp, mass_ns, theta_pb, **kwargs): return luminosity @citation_wrapper('redback') -def basic_magnetar_powered_bolometric(time, p0, bp, mass_ns, theta_pb,interaction_process=ip.Diffusion, **kwargs): +def basic_magnetar_powered_bolometric(time, p0, bp, mass_ns, theta_pb, **kwargs): """ :param time: time in days in source frame :param p0: initial spin period :param bp: polar magnetic field strength in Gauss :param mass_ns: mass of neutron star in solar masses :param theta_pb: angle between spin and magnetic field axes - :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor + :param interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :return: bolometric_luminosity """ + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + lbol = _basic_magnetar(time=time*day_to_s, p0=p0, bp=bp, mass_ns=mass_ns, theta_pb=theta_pb) - if interaction_process is not None: - interaction_class = interaction_process(time=time, luminosity=lbol, **kwargs) + if _interaction_process is not None: + interaction_class = _interaction_process(time=time, luminosity=lbol, **kwargs) lbol = interaction_class.new_luminosity return lbol @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2017ApJ...850...55N/abstract') -def basic_magnetar_powered(time, redshift, p0, bp, mass_ns, theta_pb, interaction_process=ip.Diffusion, - photosphere=photosphere.TemperatureFloor, sed=sed.Blackbody,**kwargs): +def basic_magnetar_powered(time, redshift, p0, bp, mass_ns, theta_pb,**kwargs): """ :param time: time in days in observer frame :param redshift: source redshift @@ -203,24 +205,28 @@ def basic_magnetar_powered(time, redshift, p0, bp, mass_ns, theta_pb, interactio :param bp: polar magnetic field strength in Gauss :param mass_ns: mass of neutron star in solar masses :param theta_pb: angle between spin and magnetic field axes + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor :param interaction_process: Default is Diffusion. Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. - :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used - e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor :return: flux_density or magnitude depending on output_format kwarg """ + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + _photosphere = kwargs.get("photosphere", photosphere.TemperatureFloor) + _sed = kwargs.get("sed", sed.Blackbody) + frequency = kwargs['frequency'] frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value lbol = basic_magnetar_powered_bolometric(time=time, p0=p0,bp=bp, mass_ns=mass_ns, theta_pb=theta_pb, - interaction_process=interaction_process, **kwargs) - photo = photosphere(time=time, luminosity=lbol, **kwargs) + interaction_process=_interaction_process, **kwargs) + photo = _photosphere(time=time, luminosity=lbol, **kwargs) - sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + sed_1 = _sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, frequency=frequency, luminosity_distance=dl) flux_density = sed_1.flux_density @@ -231,7 +237,7 @@ def basic_magnetar_powered(time, redshift, p0, bp, mass_ns, theta_pb, interactio return flux_density.to(uu.ABmag).value @citation_wrapper('redback') -def slsn_bolometric(time, p0, bp, mass_ns, theta_pb,interaction_process=ip.Diffusion, **kwargs): +def slsn_bolometric(time, p0, bp, mass_ns, theta_pb,**kwargs): """ Same as basic magnetar_powered but with constraint on rotational_energy/kinetic_energy and nebula phase @@ -240,19 +246,19 @@ def slsn_bolometric(time, p0, bp, mass_ns, theta_pb,interaction_process=ip.Diffu :param bp: polar magnetic field strength in Gauss :param mass_ns: mass of neutron star in solar masses :param theta_pb: angle between spin and magnetic field axes - :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor + :param interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :return: bolometric_luminosity """ + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) return basic_magnetar_powered_bolometric(time=time, p0=p0, bp=bp, mass_ns=mass_ns, - theta_pb=theta_pb, interaction_process=interaction_process, **kwargs) + theta_pb=theta_pb, interaction_process=_interaction_process, **kwargs) @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2017ApJ...850...55N/abstract') -def slsn(time, redshift, p0, bp, mass_ns, theta_pb, interaction_process=ip.Diffusion, - photosphere=photosphere.TemperatureFloor, sed=sed.CutoffBlackbody,**kwargs): +def slsn(time, redshift, p0, bp, mass_ns, theta_pb,**kwargs): """ Same as basic magnetar_powered but with constraint on rotational_energy/kinetic_energy and nebula phase @@ -262,24 +268,28 @@ def slsn(time, redshift, p0, bp, mass_ns, theta_pb, interaction_process=ip.Diffu :param bp: polar magnetic field strength in Gauss :param mass_ns: mass of neutron star in solar masses :param theta_pb: angle between spin and magnetic field axes + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor + and CutoffBlackbody: cutoff_wavelength, default is 3000 Angstrom :param interaction_process: Default is Diffusion. Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is CutoffBlackbody. - :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used - e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor - and CutoffBlackbody: cutoff_wavelength, default is 3000 Angstrom :return: flux_density or magnitude depending on output_format kwarg """ frequency = kwargs['frequency'] + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + _photosphere = kwargs.get("photosphere", photosphere.TemperatureFloor) + _sed = kwargs.get("sed", sed.CutoffBlackbody) cutoff_wavelength = kwargs.get('cutoff_wavelength', 3000) frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value + lbol = slsn_bolometric(time=time, p0=p0, bp=bp, mass_ns=mass_ns, theta_pb=theta_pb, - interaction_process=interaction_process) - photo = photosphere(time=time, luminosity=lbol, **kwargs) - sed_1 = sed(time=time, luminosity=lbol, temperature=photo.photosphere_temperature, + interaction_process=_interaction_process) + photo = _photosphere(time=time, luminosity=lbol, **kwargs) + sed_1 = _sed(time=time, luminosity=lbol, temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere,frequency=frequency, luminosity_distance=dl, cutoff_wavelength=cutoff_wavelength) @@ -291,8 +301,7 @@ def slsn(time, redshift, p0, bp, mass_ns, theta_pb, interaction_process=ip.Diffu return flux_density.to(uu.ABmag).value @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018ApJS..236....6G/abstract') -def magnetar_nickel(time, redshift, f_nickel, mej, p0, bp, mass_ns, theta_pb, interaction_process=ip.Diffusion, - photosphere=photosphere.TemperatureFloor, sed=sed.Blackbody, **kwargs): +def magnetar_nickel(time, redshift, f_nickel, mej, p0, bp, mass_ns, theta_pb, **kwargs): """ :param time: time in days in observer frame :param f_nickel: fraction of nickel mass @@ -302,15 +311,19 @@ def magnetar_nickel(time, redshift, f_nickel, mej, p0, bp, mass_ns, theta_pb, in :param bp: polar magnetic field strength in Gauss :param mass_ns: mass of neutron star in solar masses :param theta_pb: angle between spin and magnetic field axes + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor :param interaction_process: Default is Diffusion. Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. - :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used - e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor :return: flux_density or magnitude depending on output_format kwarg """ + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + _photosphere = kwargs.get("photosphere", photosphere.TemperatureFloor) + _sed = kwargs.get("sed", sed.Blackbody) + frequency = kwargs['frequency'] frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value @@ -319,13 +332,13 @@ def magnetar_nickel(time, redshift, f_nickel, mej, p0, bp, mass_ns, theta_pb, in lbol_arnett = _nickelcobalt_engine(time=time, f_nickel=f_nickel, mej=mej) lbol = lbol_mag + lbol_arnett - if interaction_process is not None: - interaction_class = interaction_process(time=time, luminosity=lbol, **kwargs) + if _interaction_process is not None: + interaction_class = _interaction_process(time=time, luminosity=lbol, **kwargs) lbol = interaction_class.new_luminosity - photo = photosphere(time=time, luminosity=lbol, **kwargs) + photo = _photosphere(time=time, luminosity=lbol, **kwargs) - sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + sed_1 = _sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, frequency=frequency, luminosity_distance=dl) flux_density = sed_1.flux_density @@ -336,19 +349,18 @@ def magnetar_nickel(time, redshift, f_nickel, mej, p0, bp, mass_ns, theta_pb, in return flux_density.to(uu.ABmag).value @citation_wrapper('redback') -def homologous_expansion_supernova_model_bolometric(time, mej, ek, interaction_process=ip.Diffusion, - **kwargs): +def homologous_expansion_supernova_model_bolometric(time, mej, ek, **kwargs): """ Assumes homologous expansion to transform kinetic energy to ejecta velocity :param time: time in days in source frame :param mej: ejecta mass in solar masses :param ek: kinetic energy in ergs - :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor 'base model' from homologous_expansion_models list + :param interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :return: bolometric_luminosity """ from redback.model_library import modules_dict # import model library in function to avoid circular dependency @@ -366,25 +378,25 @@ def homologous_expansion_supernova_model_bolometric(time, mej, ek, interaction_p kwargs['vej'] = v_ejecta kwargs['mej'] = mej - lbol = function(time, interaction_process=interaction_process, **kwargs) + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + lbol = function(time, interaction_process=_interaction_process, **kwargs) return lbol @citation_wrapper('redback') -def thin_shell_supernova_model_bolometric(time, mej, ek, interaction_process=ip.Diffusion, - **kwargs): +def thin_shell_supernova_model_bolometric(time, mej, ek, **kwargs): """ Assumes thin shell ejecta to transform kinetic energy into ejecta velocity :param time: time in days in source frame :param mej: ejecta mass in solar masses :param ek: kinetic energy in ergs - :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor, 'base model' from homologous_expansion_models list + :param interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :return: bolometric_luminosity """ from redback.model_library import modules_dict # import model library in function to avoid circular dependency @@ -402,13 +414,13 @@ def thin_shell_supernova_model_bolometric(time, mej, ek, interaction_process=ip. kwargs['vej'] = v_ejecta kwargs['mej'] = mej - lbol = function(time, interaction_process=interaction_process, **kwargs) + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + lbol = function(time, interaction_process=_interaction_process, **kwargs) return lbol @citation_wrapper('redback') -def homologous_expansion_supernova_model(time, redshift, mej, ek, interaction_process=ip.Diffusion, - photosphere=photosphere.TemperatureFloor, sed=sed.Blackbody, **kwargs): +def homologous_expansion_supernova_model(time, redshift, mej, ek, **kwargs): """ Assumes homologous expansion to transform kinetic energy to ejecta velocity @@ -416,25 +428,30 @@ def homologous_expansion_supernova_model(time, redshift, mej, ek, interaction_pr :param redshift: source redshift :param mej: ejecta mass in solar masses :param ek: kinetic energy in ergs + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor + 'base model' from homologous_expansion_models list :param interaction_process: Default is Diffusion. Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. - :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used - e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor - 'base model' from homologous_expansion_models list + :return: flux_density or magnitude depending on output_format kwarg """ + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + _photosphere = kwargs.get("photosphere", photosphere.TemperatureFloor) + _sed = kwargs.get("sed", sed.Blackbody) + frequency = kwargs['frequency'] frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value lbol = homologous_expansion_supernova_model_bolometric(time=time, mej=mej, ek=ek, - interaction_process=interaction_process, **kwargs) - photo = photosphere(time=time, luminosity=lbol, **kwargs) + interaction_process=_interaction_process, **kwargs) + photo = _photosphere(time=time, luminosity=lbol, **kwargs) - sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + sed_1 = _sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, frequency=frequency, luminosity_distance=dl) flux_density = sed_1.flux_density @@ -445,8 +462,7 @@ def homologous_expansion_supernova_model(time, redshift, mej, ek, interaction_pr return flux_density.to(uu.ABmag).value @citation_wrapper('redback') -def thin_shell_supernova_model(time, redshift, mej, ek, interaction_process=ip.Diffusion, - photosphere=photosphere.TemperatureFloor, sed=sed.Blackbody, **kwargs): +def thin_shell_supernova_model(time, redshift, mej, ek, **kwargs): """ Assumes thin shell ejecta to transform kinetic energy into ejecta velocity @@ -454,25 +470,29 @@ def thin_shell_supernova_model(time, redshift, mej, ek, interaction_process=ip.D :param redshift: source redshift :param mej: ejecta mass in solar masses :param ek: kinetic energy in ergs + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor + 'base model' from homologous_expansion_models list :param interaction_process: Default is Diffusion. Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. - :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used - e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor - 'base model' from homologous_expansion_models list :return: flux_density or magnitude depending on output_format kwarg """ + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + _photosphere = kwargs.get("photosphere", photosphere.TemperatureFloor) + _sed = kwargs.get("sed", sed.Blackbody) + frequency = kwargs['frequency'] frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value lbol = thin_shell_supernova_model_bolometric(time=time, mej=mej, ek=ek, - interaction_process=interaction_process, **kwargs) - photo = photosphere(time=time, luminosity=lbol, **kwargs) + interaction_process=_interaction_process, **kwargs) + photo = _photosphere(time=time, luminosity=lbol, **kwargs) - sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + sed_1 = _sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, frequency=frequency, luminosity_distance=dl) flux_density = sed_1.flux_density @@ -572,8 +592,7 @@ def _csm_engine(time, mej, csm_mass, vej, eta, rho, kappa, r0, **kwargs): @citation_wrapper('redback') -def csm_interaction_bolometric(time, mej, csm_mass, vej, eta, rho, kappa, r0, - interaction_process=ip.CSMDiffusion, **kwargs): +def csm_interaction_bolometric(time, mej, csm_mass, vej, eta, rho, kappa, r0, **kwargs): """ :param time: time in days in source frame :param mej: ejecta mass in solar masses @@ -583,23 +602,25 @@ def csm_interaction_bolometric(time, mej, csm_mass, vej, eta, rho, kappa, r0, :param rho: csm density profile amplitude :param kappa: opacity :param r0: radius of csm shell in AU - :param interaction_process: Default is CSMDiffusion. - Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: efficiency: in converting between kinetic energy and luminosity, default 0.5 delta: default 1, nn: default 12, If interaction process is different kwargs must include other keyword arguments that are required. + :param interaction_process: Default is CSMDiffusion. + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :return: bolometric_luminosity """ + _interaction_process = kwargs.get("interaction_process", ip.CSMDiffusion) + csm_output = _csm_engine(time=time, mej=mej, csm_mass=csm_mass, vej=vej, eta=eta, rho=rho, kappa=kappa, r0=r0, **kwargs) lbol = csm_output.lbol r_photosphere = csm_output.r_photosphere mass_csm_threshold = csm_output.mass_csm_threshold - if interaction_process is not None: - interaction_class = interaction_process(time=time, luminosity=lbol, + if _interaction_process is not None: + interaction_class = _interaction_process(time=time, luminosity=lbol, kappa=kappa, r_photosphere=r_photosphere, mass_csm_threshold=mass_csm_threshold, csm_mass=csm_mass, **kwargs) lbol = interaction_class.new_luminosity @@ -607,9 +628,7 @@ def csm_interaction_bolometric(time, mej, csm_mass, vej, eta, rho, kappa, r0, @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2013ApJ...773...76C/abstract') -def csm_interaction(time, redshift, mej, csm_mass, vej, eta, rho, kappa, r0, - interaction_process=ip.CSMDiffusion, photosphere=photosphere.TemperatureFloor, - sed=sed.Blackbody, **kwargs): +def csm_interaction(time, redshift, mej, csm_mass, vej, eta, rho, kappa, r0, **kwargs): """ :param time: time in days in observer frame :param redshift: source redshift @@ -620,26 +639,30 @@ def csm_interaction(time, redshift, mej, csm_mass, vej, eta, rho, kappa, r0, :param rho: csm density profile amplitude :param kappa: opacity :param r0: radius of csm shell in AU + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa_gamma, temperature_floor + 'base model' from homologous_expansion_models list :param interaction_process: Default is CSMDiffusion. Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. - :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used - e.g., for Diffusion and TemperatureFloor: kappa_gamma, temperature_floor - 'base model' from homologous_expansion_models list :return: flux_density or magnitude depending on output_format kwarg """ + _interaction_process = kwargs.get("interaction_process", ip.CSMDiffusion) + _photosphere = kwargs.get("photosphere", photosphere.TemperatureFloor) + _sed = kwargs.get("sed", sed.Blackbody) + frequency = kwargs['frequency'] frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value lbol = csm_interaction_bolometric(time=time, mej=mej, csm_mass=csm_mass, vej=vej, eta=eta, - rho=rho, kappa=kappa, r0=r0, interaction_process=interaction_process, **kwargs) + rho=rho, kappa=kappa, r0=r0, interaction_process=_interaction_process, **kwargs) - photo = photosphere(time=time, luminosity=lbol, vej=vej, **kwargs) + photo = _photosphere(time=time, luminosity=lbol, vej=vej, **kwargs) - sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + sed_1 = _sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, frequency=frequency, luminosity_distance=dl) flux_density = sed_1.flux_density @@ -761,51 +784,56 @@ def type_1c(time, redshift, f_nickel, mej, **kwargs): return flux_density.to(uu.ABmag).value @citation_wrapper('redback') -def general_magnetar_slsn_bolometric(time, l0, tsd, nn, interaction_process=ip.Diffusion, **kwargs): +def general_magnetar_slsn_bolometric(time, l0, tsd, nn, **kwargs): """ :param time: time in days in source frame :param l0: magnetar energy normalisation in ergs :param tsd: magnetar spin down damping timescale in source frame days :param nn: braking index - :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process, e.g., for Diffusion: kappa, kappa_gamma, vej (km/s) - :return: + :param interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. + :return: bolometric_luminosity """ + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + lbol = _magnetar_only(time=time*day_to_s, l0=l0, tsd=tsd*day_to_s, nn=nn) - if interaction_process is not None: - interaction_class = interaction_process(time=time, luminosity=lbol, **kwargs) + if _interaction_process is not None: + interaction_class = _interaction_process(time=time, luminosity=lbol, **kwargs) lbol = interaction_class.new_luminosity return lbol @citation_wrapper('redback') -def general_magnetar_slsn(time, redshift, l0, tsd, nn, interaction_process = ip.Diffusion, - photosphere = photosphere.TemperatureFloor, sed = sed.Blackbody, ** kwargs): +def general_magnetar_slsn(time, redshift, l0, tsd, nn, ** kwargs): """ :param time: time in days in observer frame :param redshift: source redshift :param l0: magnetar energy normalisation in ergs :param tsd: magnetar spin down damping timescale in source frame in days :param nn: braking index + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor :param interaction_process: Default is Diffusion. Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. - :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used - e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor :return: flux_density or magnitude depending on output_format kwarg """ + _interaction_process = kwargs.get("interaction_process", ip.Diffusion) + _photosphere = kwargs.get("photosphere", photosphere.TemperatureFloor) + _sed = kwargs.get("sed", sed.Blackbody) + frequency = kwargs['frequency'] frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value lbol = general_magnetar_slsn_bolometric(time=time, l0=l0, tsd=tsd, nn=nn, - interaction_process = interaction_process, ** kwargs) - photo = photosphere(time=time, luminosity=lbol, **kwargs) + interaction_process = _interaction_process, ** kwargs) + photo = _photosphere(time=time, luminosity=lbol, **kwargs) - sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + sed_1 = _sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, frequency = frequency, luminosity_distance = dl) flux_density = sed_1.flux_density From 489355e5ae9583bb06ed3d0afea7bd256238a85b Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Fri, 11 Mar 2022 15:42:23 +0100 Subject: [PATCH 30/37] Fix blackbody cut off --- redback/sed.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/redback/sed.py b/redback/sed.py index 7c45f37b3..5428c298c 100644 --- a/redback/sed.py +++ b/redback/sed.py @@ -56,7 +56,7 @@ def calculate_flux_density(self): # Mostly from Mosfit/SEDs cutoff_wavelength = self.cutoff_wavelength * angstrom_cgs wavelength = nu_to_lambda(self.frequency) - x_const = planck * speed_of_light * boltzmann_constant + x_const = planck * speed_of_light / boltzmann_constant flux_const = 4 * np.pi * 2*np.pi * planck * speed_of_light**2 * angstrom_cgs mask = wavelength < cutoff_wavelength @@ -75,8 +75,6 @@ def calculate_flux_density(self): norms = self.luminosity[uniq_is] / \ (flux_const / angstrom_cgs * self.r_photosphere[uniq_is]**2 * self.temperature[uniq_is]) - rp2 = self.r_photosphere[uniq_is]**2 - rp2 = rp2.reshape(lu, 1) tp = self.temperature[uniq_is].reshape(lu, 1) tp2 = tp**2 tp3 = tp**3 @@ -89,7 +87,6 @@ def calculate_flux_density(self): (6 * tp3 - np.exp(-nxcs / (cutoff_wavelength * tp)) * (nxcs ** 3 + 3 * nxcs ** 2 * cutoff_wavelength * tp + 6 * (nxcs * cutoff_wavelength**2 * tp2 + cutoff_wavelength**3 *tp3)) / cutoff_wavelength**3) / (nxcs ** 4) ), 1 ) - norms /= f_blue_reds # Apply renormalisation From 6e40d46ea7a7fd895b0883b4d63255c37c9a8d60 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Fri, 11 Mar 2022 15:46:08 +0100 Subject: [PATCH 31/37] tde analytical bolometric prior --- ...roadband_afterglow_private_data_example.py | 20 ++++++++++++------- .../priors/tde_analytical_bolometric.prior | 6 ++++++ redback/transient_models/supernova_models.py | 5 ++--- 3 files changed, 21 insertions(+), 10 deletions(-) create mode 100644 redback/priors/tde_analytical_bolometric.prior diff --git a/examples/broadband_afterglow_private_data_example.py b/examples/broadband_afterglow_private_data_example.py index 853f194db..b71d6c7b4 100644 --- a/examples/broadband_afterglow_private_data_example.py +++ b/examples/broadband_afterglow_private_data_example.py @@ -1,7 +1,7 @@ # What if your data is not public, or if you simulated it yourself. # You can fit it with redback but now by loading the transient object yourself. # This example also shows how to fit afterglow models to broadband data. -# In particular, the afterglow of GRB170817A, the GRB that accompanied the first Binary neutron star merger +# In particular, the afterglow of GRB170817A, the GRB that accompanied the first binary neutron star merger import redback import pandas as pd @@ -13,7 +13,7 @@ frequency = data['frequency'].values flux_density_err = data['flux_err'].values -# we now load the afterglow transient object. We are using flux_density data here so we need to use that data mode +# we now load the afterglow transient object. We are using flux_density data here, so we need to use that data mode data_mode = 'flux_density' # set some other useful things as variables @@ -25,17 +25,20 @@ flux_density=flux_density, flux_density_err=flux_density_err, frequency=frequency) # Now we have loaded the data up, we can plot it. +# We can pass keyword arguments here to not show, save the plot or change the aesthetics. afterglow.plot_data() afterglow.plot_multiband() # now let's actually fit it with data. We will use all the data and a gaussiancore structured jet from afterglowpy. -# Note this is not a fast example, so we will make some sampling sacrifices for speed. +# Note to make the example run quick, we will make some sampling sacrifices for speed and fix several parameters. model = 'gaussiancore' -# use default priors and 'nestle' sampler +# use default priors and 'dynesty' sampler sampler = 'dynesty' priors = redback.priors.get_priors(model=model) + +#fix the redshift priors['redshift'] = redshift # We are gonna fix some of the microphysical parameters for speed @@ -47,11 +50,14 @@ model_kwargs = dict(frequency=frequency, output_format='flux_density') -# returns a supernova result object +# Fitting returns a afterglow result object, which can be used for plotting or doing further analysis. result = redback.fit_model(transient=afterglow, model=model, sampler=sampler, model_kwargs=model_kwargs, - prior=priors, sample='rslice', nlive=50, dlogz=10, resume=True) + prior=priors, sample='rslice', nlive=500, resume=True) # plot corner result.plot_corner() -# plot multiband lightcurve. This will plot a panel for every unique frequency +# plot multiband lightcurve. +# This will plot a panel for every unique frequency, +# along with the fitted lightcurve from a 100 random realisations randomly drawn from the prior. +# We can change other settings by passing in different keyword arguments here. result.plot_multiband_lightcurve(random_models=100) diff --git a/redback/priors/tde_analytical_bolometric.prior b/redback/priors/tde_analytical_bolometric.prior new file mode 100644 index 000000000..a2115e55b --- /dev/null +++ b/redback/priors/tde_analytical_bolometric.prior @@ -0,0 +1,6 @@ +mej = LogUniform(1e-4, 100, 'mej', latex_label = r'$M_{\mathrm{ej}} [M_{\odot}]$') +vej = LogUniform(1e3, 1e5, 'vej', latex_label = r'$v_{\mathrm{ej}} [km/s]$') +kappa = Uniform(0.05, 2, 'kappa', latex_label = r'$\kappa$') +kappa_gamma = Uniform(1e-4, 1e4, 'kappa_gamma', latex_label = r'$\kappa_{\gamma}$') +l0 = LogUniform(1e51, 1e58, "l0", latex_label="$l_0$") +t_0 = LogUniform(1e-4, 1e2, "t_0", latex_label="$t_0$") \ No newline at end of file diff --git a/redback/transient_models/supernova_models.py b/redback/transient_models/supernova_models.py index 86de784ef..851ffd72a 100644 --- a/redback/transient_models/supernova_models.py +++ b/redback/transient_models/supernova_models.py @@ -254,7 +254,7 @@ def slsn_bolometric(time, p0, bp, mass_ns, theta_pb,**kwargs): """ _interaction_process = kwargs.get("interaction_process", ip.Diffusion) return basic_magnetar_powered_bolometric(time=time, p0=p0, bp=bp, mass_ns=mass_ns, - theta_pb=theta_pb, interaction_process=_interaction_process, **kwargs) + theta_pb=theta_pb, **kwargs) @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2017ApJ...850...55N/abstract') @@ -286,8 +286,7 @@ def slsn(time, redshift, p0, bp, mass_ns, theta_pb,**kwargs): frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value - lbol = slsn_bolometric(time=time, p0=p0, bp=bp, mass_ns=mass_ns, theta_pb=theta_pb, - interaction_process=_interaction_process) + lbol = slsn_bolometric(time=time, p0=p0, bp=bp, mass_ns=mass_ns, theta_pb=theta_pb, **kwargs) photo = _photosphere(time=time, luminosity=lbol, **kwargs) sed_1 = _sed(time=time, luminosity=lbol, temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere,frequency=frequency, luminosity_distance=dl, From 6c5b0234844f3286334fad1efef0d11fb1625b47 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Fri, 11 Mar 2022 15:51:40 +0100 Subject: [PATCH 32/37] add pyqt5 --- optional_requirements.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/optional_requirements.txt b/optional_requirements.txt index 26ab1ebb7..eff8f06ae 100644 --- a/optional_requirements.txt +++ b/optional_requirements.txt @@ -3,4 +3,5 @@ sncosmo lalsimulation nestle sherpa -kilonova-heating-rate \ No newline at end of file +kilonova-heating-rate +PyQt5 \ No newline at end of file From 545bba4cb5ab85caaa9815960cac1b613412aa82 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Fri, 11 Mar 2022 16:30:47 +0100 Subject: [PATCH 33/37] change A_1 to a_1 --- redback/transient_models/magnetar_models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/redback/transient_models/magnetar_models.py b/redback/transient_models/magnetar_models.py index aff5c6450..6649ed69d 100644 --- a/redback/transient_models/magnetar_models.py +++ b/redback/transient_models/magnetar_models.py @@ -75,7 +75,7 @@ def evolving_magnetar(time, a_1, alpha_1, mu0, muinf, p0, sinalpha0, tm, II, **k :param kwargs: key word argument for handling plotting :return: luminosity (depending on scaling) as a function of time. """ - pl = one_component_fireball_model(time=time, a_1=A_1, alpha_1=alpha_1) + pl = one_component_fireball_model(time=time, a_1=a_1, alpha_1=alpha_1) magnetar = _evolving_magnetar_only(time=time, mu0=mu0, muinf=muinf, p0=p0, sinalpha0=sinalpha0, tm=tm, II=II) return pl + magnetar @@ -409,7 +409,7 @@ def collapsing_radiative_losses(time, a_1, alpha_1, l0, tau, nn, tcol, kappa, t0 :param l0: initial luminosity parameter :param tau: spin-down damping timescale :param nn: braking index - :param tcol: collapse time in seconds + :param tcol: collapse time in seconds :param kappa: radiative efficiency :param t0: time for radiative losses to start in seconds :param log_e0: log10 E0 to connect curvature effect energy smoothly with transition point energy From 715368bf255ad28eb6cd062b5b70cc4d5fbf2f05 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Fri, 11 Mar 2022 16:31:12 +0100 Subject: [PATCH 34/37] fix selenium version to one compatible with phantomjs --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 9c74eb8e1..82ef688ce 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,7 +2,7 @@ setuptools numpy pandas scipy -selenium +selenium==3.8.0 bilby matplotlib astropy From 57f79baed19b6269a266605212d7ac2de940aad3 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Fri, 11 Mar 2022 16:32:26 +0100 Subject: [PATCH 35/37] Get examples running --- examples/fit_your_own_model_example.py | 2 +- examples/kilonova_example.py | 8 ++++---- examples/magnetar_boosted_example.py | 7 ++++--- examples/magnetar_example.py | 7 ++++--- examples/supernova_example.py | 3 +-- 5 files changed, 14 insertions(+), 13 deletions(-) diff --git a/examples/fit_your_own_model_example.py b/examples/fit_your_own_model_example.py index 09a33d59c..d4a7da8e7 100644 --- a/examples/fit_your_own_model_example.py +++ b/examples/fit_your_own_model_example.py @@ -44,4 +44,4 @@ def my_favourite_model(time, l0, alpha, **kwargs): result = redback.fit_model(model=model, sampler='dynesty', nlive=200, transient=afterglow, prior=priors, sample='rslice', resume=False) -result.plot_lightcurve(random_models=1000, model=my_favourite_model) +result.plot_lightcurve(random_models=100, model=my_favourite_model) diff --git a/examples/kilonova_example.py b/examples/kilonova_example.py index 49d8a6f7f..79a861d0c 100644 --- a/examples/kilonova_example.py +++ b/examples/kilonova_example.py @@ -24,10 +24,10 @@ priors['redshift'] = 1e-2 model_kwargs = dict(frequency=kilonova.filtered_frequencies, output_format='flux_density') -# result = redback.fit_model(transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, -# prior=priors, sample='rslice', nlive=200, resume=True) -result = redback.result.read_in_result("kilonova/one_component_kilonova_model/at2017gfo_result.json") -# result.plot_corner() +result = redback.fit_model(transient=kilonova, model=model, sampler=sampler, model_kwargs=model_kwargs, + prior=priors, sample='rslice', nlive=200, resume=True) +result.plot_corner() # returns a Kilonova result object result.plot_lightcurve(plot_show=False) +# Even though we only fit the 'g' band, we can still plot the fit for different bands. result.plot_multiband_lightcurve(filters=["g", "r", "i", "z", "y", "J"], plot_show=False) diff --git a/examples/magnetar_boosted_example.py b/examples/magnetar_boosted_example.py index ab4331739..63aa7d76c 100644 --- a/examples/magnetar_boosted_example.py +++ b/examples/magnetar_boosted_example.py @@ -12,7 +12,7 @@ # gets the magnitude data for AT2017gfo, the KN associated with GW170817 data = redback.get_data.get_kilonova_data_from_open_transient_catalog_data(transient=kne) # creates a GRBDir with GRB -kilonova = redback.kilonova.Kilonova.from_open_access_catalogue(name=kne, data_mode="flux_density") +kilonova = redback.kilonova.Kilonova.from_open_access_catalogue(name=kne, data_mode="flux_density", active_bands=['g']) # kilonova.flux_density_data = True fig, axes = plt.subplots(3, 2, sharex=True, sharey=True, figsize=(12, 8)) bands = ["g"] @@ -32,5 +32,6 @@ prior=priors, sample='rslice', nlive=200, resume=True) # result.plot_corner() # returns a Kilonova result object -# result.plot_lightcurve(random_models=1000) -result.plot_multiband_lightcurve() +# result.plot_lightcurve(random_models=50) +# Even though we only fit the 'g' band, we can still plot the fit for different bands. +result.plot_multiband_lightcurve(filters=["g", "r", "i", "z", "y", "J"]) diff --git a/examples/magnetar_example.py b/examples/magnetar_example.py index 1d43eb2da..529888cc9 100644 --- a/examples/magnetar_example.py +++ b/examples/magnetar_example.py @@ -1,4 +1,5 @@ import redback +import bilby # We implemented many models implemented, including # afterglow/magnetar varieties/n_dimensional_fireball/shapelets/band function/kilonova/SNe/TDE @@ -21,8 +22,8 @@ # use default priors priors = redback.priors.get_priors(model=model, data_mode='luminosity') -# alternatively can pass in some priors -# priors = {} +# alternatively can create a dictionary of priors in some priors +# priors = bilby.core.prior.PriorDict() # priors['A_1'] = bilby.core.prior.LogUniform(1e-15, 1e15, 'A_1', latex_label = r'$A_{1}$') # priors['alpha_1'] = bilby.core.prior.Uniform(-7, -1, 'alpha_1', latex_label = r'$\alpha_{1}$') # priors['p0'] = bilby.core.prior.Uniform(0.7e-3, 0.1, 'p0', latex_label = r'$P_{0} [s]$') @@ -37,4 +38,4 @@ result = redback.fit_model(model=model, sampler='dynesty', nlive=200, transient=afterglow, prior=priors, sample='rslice', resume=True) -result.plot_lightcurve(random_models=1000) +result.plot_lightcurve(random_models=100) diff --git a/examples/supernova_example.py b/examples/supernova_example.py index 00e6706db..a8a11a197 100644 --- a/examples/supernova_example.py +++ b/examples/supernova_example.py @@ -11,8 +11,7 @@ data = redback.get_data.get_supernova_data_from_open_transient_catalog_data(sne) -supernova = redback.supernova.Supernova.from_open_access_catalogue(name=sne, data_mode='flux_density') -bands = ["g", "i"] +supernova = redback.supernova.Supernova.from_open_access_catalogue(name=sne, data_mode='flux_density', active_bands=["g", "i"]) supernova.plot_multiband(filters=["J", "H", "g", "i"]) # use default priors From 98636dd82501c4b1ac5621dff80c8c9c2e7246cd Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Fri, 11 Mar 2022 17:06:59 +0100 Subject: [PATCH 36/37] get tde example running --- examples/tde_example.py | 2 +- redback/transient_models/tde_models.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/tde_example.py b/examples/tde_example.py index 80cd87732..e4d32befa 100644 --- a/examples/tde_example.py +++ b/examples/tde_example.py @@ -22,7 +22,7 @@ priors = redback.priors.get_priors(model=model) # we know the redshift for PS18kh so we just fix the prior for the redshift to the known value. # We can do this through the metadata that was downloaded alongside the data, or if you just know it. -# priors['redshift'] = tidal_disruption_event.redshift +priors['redshift'] = 0.07 model_kwargs = dict(frequency=tidal_disruption_event.filtered_frequencies, output_format='flux_density') diff --git a/redback/transient_models/tde_models.py b/redback/transient_models/tde_models.py index f039aecdf..4f1f9ae94 100644 --- a/redback/transient_models/tde_models.py +++ b/redback/transient_models/tde_models.py @@ -71,7 +71,7 @@ def tde_analytical(time, redshift, l0, t_0, **kwargs): frequency=frequency, luminosity_distance=dl, cutoff_wavelength=cutoff_wavelength, luminosity=lbol) flux_density = sed_1.flux_density - + flux_density = np.nan_to_num(flux_density) if kwargs['output_format'] == 'flux_density': return flux_density.to(uu.mJy).value elif kwargs['output_format'] == 'magnitude': From a7477ccfa82b735f6e1203ca49abb087bb146ba9 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Fri, 11 Mar 2022 17:11:52 +0100 Subject: [PATCH 37/37] Version 0.2 release --- CHANGELOG.md | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 CHANGELOG.md diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 000000000..82f8154e4 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,9 @@ +# All notable changes will be documented in this file + +## [0.2.0] 2022-03-11 +Version 1.1.5 release of redback + +### Added +- Basic functionality +- Several examples +- Docs \ No newline at end of file