Skip to content

Commit

Permalink
updated documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
@ACCarnall
ACCarnall committed Apr 5, 2018
1 parent d1c9511 commit ad57f6d
Show file tree
Hide file tree
Showing 15 changed files with 477 additions and 184 deletions.
76 changes: 37 additions & 39 deletions bagpipes/fit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,12 +20,6 @@
import model_manager as models
import model_galaxy

"""
# Turn off annoying dividion by zero errors when MultiNest is running
import warnings
warnings.filterwarnings("ignore")
print("Bagpipes: Warning, Python warnings are being ignored.")
"""


class Fit:
Expand Down Expand Up @@ -161,7 +155,6 @@ def process_fit_instructions(self):
self.priors.append(self.fit_instructions[fit_param.split(":")[0]][fit_param.split(":")[1] + "prior"])

else:
print("Bagpipes: Warning, no prior specified on " + fit_param + ", adopting a uniform prior.")
self.priors.append("uniform")

"""
Expand All @@ -179,9 +172,11 @@ def process_fit_instructions(self):
def fit(self, verbose=False, sampling_efficiency="model", n_live=400, const_efficiency_mode=False):
""" Fit the specified model to the input galaxy data. """

print("\nBagpipes: fitting object " + self.Galaxy.ID + "\n")

pmn.run(self.get_lnprob, self.prior_transform, self.ndim, const_efficiency_mode = const_efficiency_mode, importance_nested_sampling = False, verbose = verbose, sampling_efficiency = sampling_efficiency, n_live_points = n_live, outputfiles_basename=models.working_dir + "/pipes/pmn_chains/" + self.run + "/" + self.Galaxy.ID + "-")

a = pmn.Analyzer(n_params = self.ndim, outputfiles_basename=models.working_dir + "/pipes/pmn_chains/" + self.run + "/" + self.Galaxy.ID + "-")
a = pmn.Analyzer(n_params = self.ndim, outputfiles_basename=models.working_dir + "/pipes/pmn_chains/" + self.run + "/" + self.Galaxy.ID + "-", verbose=False)

s = a.get_stats()

Expand Down Expand Up @@ -213,11 +208,10 @@ def fit(self, verbose=False, sampling_efficiency="model", n_live=400, const_effi

self.min_chisq_red = self.min_chisq/float(self.ndof)

if verbose == True:
print("\nBagpipes: Confidence interval:")
for x in range(self.ndim):
print(str(np.round(self.conf_int[x], 4)), np.round(self.posterior_median[x], 4), self.fit_params[x])
print("\n")
print("\nBagpipes: fitting complete, confidence interval")
for x in range(self.ndim):
print(str(np.round(self.conf_int[x], 4)), np.round(self.posterior_median[x], 4), self.fit_params[x])
print("\n")

self.get_model(self.posterior_median)
self.get_post_info()
Expand Down Expand Up @@ -390,7 +384,7 @@ def get_post_info(self):

nsamples = self.posterior["samples"].shape[0]

self.posterior["sfh"] = np.zeros((nsamples, self.Model.sfh.ages.shape[0]))
self.posterior["sfh"] = np.zeros((self.Model.sfh.ages.shape[0], nsamples))
self.posterior["sfr"] = np.zeros(nsamples)
self.posterior["tmw"] = np.zeros(nsamples)
self.posterior["UVJ"] = np.zeros((3, nsamples))
Expand Down Expand Up @@ -421,8 +415,8 @@ def get_post_info(self):
for i in range(nsamples):
self.get_model(self.posterior["samples"][i,:])

self.posterior["sfh"][i,:] = self.Model.sfh.sfr
self.posterior["sfr"][i] = self.posterior["sfh"][i,0]
self.posterior["sfh"][:,i] = self.Model.sfh.sfr
self.posterior["sfr"][i] = self.posterior["sfh"][0,i]
self.posterior["tmw"][i] = np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - (10**-9)*np.sum(self.Model.sfh.sfr*self.Model.sfh.ages*self.Model.sfh.age_widths)/np.sum(self.Model.sfh.sfr*self.Model.sfh.age_widths)
self.posterior["living_stellar_mass"]["total"][i] = self.Model.living_stellar_mass["total"]

Expand Down Expand Up @@ -633,8 +627,6 @@ def plot_corner(self, param_names_tolog=[], truths=None, ranges=None, return_fig
plot_range.append(self.fit_limits[i])

if self.fit_params[i] in param_names_tolog:
#param_names_toplot[-1] = "log_" + param_names_toplot[-1]
#param_truths_toplot[-1] = np.log10(param_truths_toplot[-1])
params_tolog.append(len(param_cols_toplot)-1)
plot_range[-1] = (np.log10(plot_range[-1][0]), np.log10(plot_range[-1][1]))

Expand All @@ -649,9 +641,7 @@ def plot_corner(self, param_names_tolog=[], truths=None, ranges=None, return_fig
if param_names_toplot[i] in reference_param_names:
param_names_toplot[i] = latex_param_names[reference_param_names.index(param_names_toplot[i])]

#ranges = [(0., 0.5), (4., np.interp(self.model_components["redshift"], models.z_array, models.age_at_z)), (-2, 3), (0.2, 0.8), (10.15, 10.65), (0.5, 3)]

fig = corner.corner(self.posterior["samples"][:,param_cols_toplot], labels=param_names_toplot, quantiles=[0.16, 0.5, 0.84], show_titles=True, title_kwargs={"fontsize": 16}, smooth="1.5", smooth1d="0.5", truths=truths, range=ranges)#truths=param_truths_toplot,
fig = corner.corner(self.posterior["samples"][:,param_cols_toplot], labels=param_names_toplot, quantiles=[0.16, 0.5, 0.84], show_titles=True, title_kwargs={"fontsize": 14}, smooth=1.5, smooth1d=0.5, truths=truths, range=ranges)

sfh_ax = fig.add_axes([0.65, 0.59, 0.32, 0.15], zorder=10)
sfr_ax = fig.add_axes([0.82, 0.82, 0.15, 0.15], zorder=10)
Expand All @@ -663,13 +653,13 @@ def plot_corner(self, param_names_tolog=[], truths=None, ranges=None, return_fig
self.plot_sfh_post(sfh_ax, style="step")

# Plot the current star formation rate posterior
sfr_ax.hist(self.posterior["sfr"], bins=15, color="white", normed=True, histtype="step", edgecolor="black")
sfr_ax.hist(self.posterior["sfr"], bins=20, color="white", density=True, histtype="step", edgecolor="black", lw=1.5, range=(np.max([0., np.mean(self.posterior["sfr"]) - 3*np.std(self.posterior["sfr"])]), np.mean(self.posterior["sfr"]) + 3*np.std(self.posterior["sfr"])))
sfr_ax.set_xlabel("$\mathrm{SFR\ /\ M_\odot\ yr^{-1}}$")
sfr_ax.set_xlim(np.max([0., np.mean(self.posterior["sfr"]) - 3*np.std(self.posterior["sfr"])]), np.mean(self.posterior["sfr"]) + 3*np.std(self.posterior["sfr"]))
sfr_ax.set_yticklabels([])

# Plot the mass weighted age posterior
tmw_ax.hist(self.posterior["tmw"], bins=15, color="white", normed=True, histtype="step", edgecolor="black")
tmw_ax.hist(self.posterior["tmw"], bins=20, color="white", density=True, histtype="step", edgecolor="black", lw=1.5, range=(np.max([0., np.mean(self.posterior["tmw"]) - 3*np.std(self.posterior["tmw"])]), np.mean(self.posterior["tmw"]) + 3*np.std(self.posterior["tmw"])))
tmw_ax.set_xlabel("$t(z_\mathrm{form})\ /\ \mathrm{Gyr}$")
tmw_ax.set_xlim(np.mean(self.posterior["tmw"]) - 3*np.std(self.posterior["tmw"]), np.mean(self.posterior["tmw"]) + 3*np.std(self.posterior["tmw"]))
tmw_ax.set_yticklabels([])
Expand All @@ -684,8 +674,8 @@ def plot_corner(self, param_names_tolog=[], truths=None, ranges=None, return_fig
tmw_ax.axvline(np.percentile(self.posterior["tmw"], 84), linestyle="--", color="black")
#mwa_ax.axvline(self.mwa_maxprob, color="#4682b4")

fig.text(0.725, 0.978, "$t(z_\mathrm{form})\ /\ \mathrm{Gyr} =\ " + str(np.round(np.percentile(self.posterior["tmw"], 50), 2)) + "^{+" + str(np.round(np.percentile(self.posterior["tmw"], 84) - np.percentile(self.posterior["tmw"], 50), 2)) + "}_{-" + str(np.round(np.percentile(self.posterior["tmw"], 50) - np.percentile(self.posterior["tmw"], 16), 2)) + "}$", horizontalalignment = "center")
fig.text(0.895, 0.978, "$\mathrm{SFR\ /\ M_\odot\ yr^{-1}}\ =\ " + str(np.round(np.percentile(self.posterior["sfr"], 50), 2)) + "^{+" + str(np.round(np.percentile(self.posterior["sfr"], 84) - np.percentile(self.posterior["sfr"], 50), 2)) + "}_{-" + str(np.round(np.percentile(self.posterior["sfr"], 50) - np.percentile(self.posterior["sfr"], 16), 2)) + "}$", horizontalalignment = "center")
fig.text(0.725, 0.978, "$t(z_\mathrm{form})\ /\ \mathrm{Gyr} =\ " + str(np.round(np.percentile(self.posterior["tmw"], 50), 2)) + "^{+" + str(np.round(np.percentile(self.posterior["tmw"], 84) - np.percentile(self.posterior["tmw"], 50), 2)) + "}_{-" + str(np.round(np.percentile(self.posterior["tmw"], 50) - np.percentile(self.posterior["tmw"], 16), 2)) + "}$", horizontalalignment = "center", size=14)
fig.text(0.895, 0.978, "$\mathrm{SFR\ /\ M_\odot\ yr^{-1}}\ =\ " + str(np.round(np.percentile(self.posterior["sfr"], 50), 2)) + "^{+" + str(np.round(np.percentile(self.posterior["sfr"], 84) - np.percentile(self.posterior["sfr"], 50), 2)) + "}_{-" + str(np.round(np.percentile(self.posterior["sfr"], 50) - np.percentile(self.posterior["sfr"], 16), 2)) + "}$", horizontalalignment = "center", size=14)

fig.savefig(models.working_dir + "/pipes/plots/" + self.run + "/" + self.Galaxy.ID + "_corner.pdf")

Expand All @@ -697,7 +687,15 @@ def plot_corner(self, param_names_tolog=[], truths=None, ranges=None, return_fig



def plot_sfh_post(self, sfh_ax, style="smooth"):
def plot_sfh_post(self, sfh_ax, style="smooth", colorscheme="bw"):

color1 = "black"
color2 = "gray"

if colorscheme == "irnbru":
color1 = "darkorange"
color2 = "navajowhite"

if style == "step":
# Generate and populate sfh arrays which allow the SFH to be plotted with straight lines across bins of SFH
sfh_x = np.zeros(2*self.Model.sfh.ages.shape[0])
Expand All @@ -709,20 +707,20 @@ def plot_sfh_post(self, sfh_ax, style="smooth"):

sfh_x[2*j] = self.Model.sfh.age_lhs[j]

sfh_y[2*j] = np.median(self.posterior["sfh"][:,j])
sfh_y[2*j + 1] = np.median(self.posterior["sfh"][:,j])
sfh_y[2*j] = np.median(self.posterior["sfh"][j,:])
sfh_y[2*j + 1] = np.median(self.posterior["sfh"][j,:])

sfh_y_low[2*j] = np.percentile(self.posterior["sfh"][:,j], 16)
sfh_y_low[2*j + 1] = np.percentile(self.posterior["sfh"][:,j], 16)
sfh_y_low[2*j] = np.percentile(self.posterior["sfh"][j,:], 16)
sfh_y_low[2*j + 1] = np.percentile(self.posterior["sfh"][j,:], 16)

if sfh_y_low[2*j] < 0:
sfh_y_low[2*j] = 0.

if sfh_y_low[2*j+1] < 0:
sfh_y_low[2*j+1] = 0.

sfh_y_high[2*j] = np.percentile(self.posterior["sfh"][:,j], 84)
sfh_y_high[2*j + 1] = np.percentile(self.posterior["sfh"][:,j], 84)
sfh_y_high[2*j] = np.percentile(self.posterior["sfh"][j,:], 84)
sfh_y_high[2*j + 1] = np.percentile(self.posterior["sfh"][j,:], 84)

if j == self.Model.sfh.sfr.shape[0]-1:
sfh_x[-1] = self.Model.sfh.age_lhs[-1] + 2*(self.Model.sfh.ages[-1] - self.Model.sfh.age_lhs[-1])
Expand All @@ -731,16 +729,16 @@ def plot_sfh_post(self, sfh_ax, style="smooth"):
sfh_x[2*j + 1] = self.Model.sfh.age_lhs[j+1]

# Plot the SFH
sfh_ax.fill_between(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - sfh_x*10**-9, sfh_y_low, sfh_y_high, color="gray", alpha=0.75)
sfh_ax.plot(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - sfh_x*10**-9, sfh_y, color="black", zorder=10)
sfh_ax.plot(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - sfh_x*10**-9, sfh_y_high, color="gray")
sfh_ax.plot(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - sfh_x*10**-9, sfh_y_low, color="gray")
sfh_ax.fill_between(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - sfh_x*10**-9, sfh_y_low, sfh_y_high, color=color2, alpha=0.75)
sfh_ax.plot(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - sfh_x*10**-9, sfh_y, color=color1, zorder=10)
sfh_ax.plot(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - sfh_x*10**-9, sfh_y_high, color=color2)
sfh_ax.plot(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - sfh_x*10**-9, sfh_y_low, color=color2)
sfh_ax.set_ylim(0, 1.1*np.max(sfh_y_high))

elif style == "smooth":
sfh_ax.fill_between(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - self.Model.sfh.ages*10**-9, np.percentile(self.posterior["sfh"], 16, axis=0), np.percentile(self.posterior["sfh"], 84, axis=0), color="navajowhite")
sfh_ax.plot(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - self.Model.sfh.ages*10**-9, np.percentile(self.posterior["sfh"], 50, axis=0), color="darkorange")
sfh_ax.set_ylim(0, 1.1*np.max(np.percentile(self.posterior["sfh"], 84, axis=0)))
sfh_ax.fill_between(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - self.Model.sfh.ages*10**-9, np.percentile(self.posterior["sfh"], 16, axis=1), np.percentile(self.posterior["sfh"], 84, axis=1), color=color2)
sfh_ax.plot(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z) - self.Model.sfh.ages*10**-9, np.percentile(self.posterior["sfh"], 50, axis=1), color=color1)
sfh_ax.set_ylim(0, 1.1*np.max(np.percentile(self.posterior["sfh"], 84, axis=1)))

sfh_ax.set_xlim(np.interp(self.model_components["redshift"], models.z_array, models.age_at_z), 0)

Expand Down
20 changes: 10 additions & 10 deletions bagpipes/galaxy.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,9 @@ class Galaxy:
----------
ID : str
string denoting the ID of the object to be loaded. This will be passed to data_load_function.
string denoting the ID of the object to be loaded. This will be passed to load_data.
data_load_function : function
load_data : function
Function which takes ID, filtlist as its two arguments and returns the model spectrum and photometry. Spectrum should come first and
be an array with a column of wavelengths in Angstroms, a column of fluxes in erg/s/cm^2/A and a column of flux errors in the same
units. Photometry should come second and be an array with a column of fluxes in microjanskys and a column of flux errors in the
Expand All @@ -28,15 +28,15 @@ class Galaxy:
The name of the filtlist, must be specified if photometry is to be loaded.
spectrum_exists : bool
If you do not have a spectrum for this object, set this to False. In this case, data_load_function should only return photometry.
If you do not have a spectrum for this object, set this to False. In this case, load_data should only return photometry.
photometry_exists : bool
If you do not have photometry for this object, set this to False. In this case, data_load_function should only return a spectrum.
If you do not have photometry for this object, set this to False. In this case, load_data should only return a spectrum.
"""

def __init__(self, ID, data_load_function, filtlist=None, spectrum_exists=True, photometry_exists=True):
def __init__(self, ID, load_data, filtlist=None, spectrum_exists=True, photometry_exists=True):

out_units="ergscma"
no_of_spectra=1
Expand All @@ -53,25 +53,25 @@ def __init__(self, ID, data_load_function, filtlist=None, spectrum_exists=True,

elif spectrum_exists == True and photometry_exists == False:
if self.no_of_spectra == 1:
self.spectrum = data_load_function(self.ID, self.filtlist)
self.spectrum = load_data(self.ID, self.filtlist)

else:
data = data_load_function(self.ID, self.filtlist)
data = load_data(self.ID, self.filtlist)
self.spectrum = data[0]
self.extra_spectra = []
for i in range(len(data)-1):
self.extra_spectra.append(data[i+1])

elif spectrum_exists == False and photometry_exists == True:
photometry_nowavs = data_load_function(self.ID, self.filtlist)
photometry_nowavs = load_data(self.ID, self.filtlist)
self.no_of_spectra = 0

else:
if self.no_of_spectra == 1:
self.spectrum, photometry_nowavs = data_load_function(self.ID, self.filtlist)
self.spectrum, photometry_nowavs = load_data(self.ID, self.filtlist)

else:
data = data_load_function(self.ID, self.filtlist)
data = load_data(self.ID, self.filtlist)
self.spectrum = data[0]
self.extra_spectra = []
photometry_nowavs = data[-1]
Expand Down
2 changes: 2 additions & 0 deletions bagpipes/model_galaxy.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@
import star_formation_history
import model_manager as models

# Ignore division by zero and overflow warnings
np.seterr(divide='ignore', invalid='ignore', over="ignore")

class Model_Galaxy:
""" Build a model galaxy spectrum.
Expand Down
6 changes: 3 additions & 3 deletions docs/filtlists.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,13 +3,13 @@
Defining filter lists
=====================

Bagpipes uses filter lists to load up filter curve files which are in turn used to generate photometry.
Bagpipes uses filter lists to load up filter curve files, which are in turn used to generate photometry.

To define a filter list (referred to as ``filtlist`` within the code) you'll first have to set up some of the :ref:`directory structure <directory-structure>`. Specifically, within your working directory you'll have to make a ``pipes/`` directory, and within ``pipes/`` a ``filters/`` directory.

Now, within the ``pipes/filters/`` directory you should create a file called ``<filter list name>.filtlist``. For example, if I wanted to set up a PanSTARRS filter list, I could call my file ``PanSTARRS.filtlist``.
Now, within the ``pipes/filters/`` directory create a file called ``<filter list name>.filtlist``. For example, if I wanted to set up a PanSTARRS filter list, I could call my file ``PanSTARRS.filtlist``.

In this file, you'll have to add paths from the ``pipes/filters/`` directory to the locations the filter curves you want to include are stored. In order to find the curves you want I recommend the `SVO filter profile service <http://svo2.cab.inta-csic.es/svo/theory/fps>`_. Bagpipes expects filter curve files to contain a column of wavelengths in Angstroms followed by a column of relative transmission values (normalisation not important).
In this file, you'll have to add paths from the ``pipes/filters/`` directory to the locations the filter curves you want to include are stored. In order to obtain these I recommend the `SVO filter profile service <http://svo2.cab.inta-csic.es/svo/theory/fps>`_. Bagpipes expects filter curve files to contain a column of wavelengths in Angstroms followed by a column of relative transmission values (normalisation not important).

For example, if you downloaded the PS1 grizy filters and put them in a directory called ``PanSTARRS/`` within the ``pipes/filters/`` directory, you'd need the following in ``PanSTARRS.filtlist``:

Expand Down
Loading

0 comments on commit ad57f6d

Please sign in to comment.