Version 0.1.0

bagpipes/IGM_Inoue2014.py

@@ -1,11 +1,11 @@
-from __future__ import print_function
+from __future__ import print_function, division, absolute_import
 
 import numpy as np 
 import matplotlib.pyplot as plt 
 import sys
 import os
 
-import model_manager as models
+from . import model_manager as models
 
 from astropy.io import fits
 

bagpipes/__init__.py

@@ -3,7 +3,7 @@
 from .model_galaxy import Model_Galaxy
 from .star_formation_history import Star_Formation_History
 from .fit import Fit
-from .compare_fits import compare_fits
-from .catalogue_fit import *
+from .compare_fits import Compare_Fits
+from .catalogue_fit import Catalogue_Fit
 from .make_cloudy_models import *
 from .IGM_Inoue2014 import *

bagpipes/catalogue_fit.py

@@ -1,4 +1,4 @@
-from __future__ import print_function
+from __future__ import print_function, division, absolute_import
 
 import bagpipes as pipes
 import numpy as np
@@ -10,9 +10,9 @@
 from astropy.io import fits
 from subprocess import call
 
-import model_manager as models
-from galaxy import Galaxy
-from fit import Fit
+from . import model_manager as models
+from .galaxy import Galaxy
+from .fit import Fit
 
 
 

bagpipes/compare_fits.py

@@ -1,3 +1,5 @@
+from __future__ import print_function, division, absolute_import
+
 import numpy as np 
 import corner
 import sys
@@ -6,7 +8,7 @@
 import os
 
 
-def compare_fits(fit1, fit2, param_names_tolog=[], truths=None, comp_run="."):
+def Compare_Fits(fit1, fit2, param_names_tolog=[], truths=None, comp_run="."):
 
 	colour1 = "darkorange"
 	colour1_2 = "navajowhite"

bagpipes/fit.py

@@ -1,4 +1,4 @@
-from __future__ import print_function
+from __future__ import print_function, division, absolute_import
 
 import numpy as np
 import os
@@ -17,8 +17,8 @@
 except:
 	print("Bagpipes: MultiNest/PyMultiNest not installed, fitting will not be available.")
 
-import model_manager as models
-import model_galaxy
+from . import model_manager as models
+from . import model_galaxy
 
 
 
@@ -29,13 +29,13 @@ class Fit:
 	Parameters
 	----------
 
-	Galaxy : Bagpipes Galaxy object
+	Galaxy : bagpipes.Galaxy
 		A Galaxy object containing the photomeric and/or spectroscopic data you wish to fit.
 
 	fit_instructions : dict
 		A dictionary containing the details of the model to be fitted to the data.
 
-	run : string (optional)
+	run : string - optional
 		The subfolder into which outputs will be saved, useful e.g. for fitting more than one model configuration to the same data.
 
 	"""
@@ -78,13 +78,13 @@ def __init__(self, Galaxy, fit_instructions, run="."):
 
 		# posterior: Will be used to store posterior samples.
 		self.posterior = {}
-		"""
+
 		# extra_models: a list into which extra model items will be placed if the galaxy object has more thna one spectrum
-		if self.Galaxy.no_of_spectra > 1:
-			self.extra_models = []
-			for i in range(self.Galaxy.no_of_spectra-1):
-				self.extra_models.append(None)
-		"""
+		#if self.Galaxy.no_of_spectra > 1:
+		#	self.extra_models = []
+		#	for i in range(self.Galaxy.no_of_spectra-1):
+		#		self.extra_models.append(None)
+
 		#Set up directories to contain the outputs, if they don't already exist
 		if self.run is not ".":
 
@@ -123,7 +123,7 @@ def __init__(self, Galaxy, fit_instructions, run="."):
 
 
 	def process_fit_instructions(self):
-		"""Sets up the class by generating relevant variables from the input fit_instructions dictionary."""
+		# Sets up the class by generating relevant variables from the input fit_instructions dictionary.
 
 		for key in list(self.fit_instructions):
 			if isinstance(self.fit_instructions[key], tuple):
@@ -170,7 +170,24 @@ 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. """
+		""" Fit the specified model to the input galaxy data using the MultiNest algorithm. 
+
+		Parameters
+		----------
+
+		verbose : bool - optional
+			Set to True to get extra progress updates from MultiNest.
+
+		sampling_efficiency : str - optional
+			MultiNest sampling efficiency parameter (see MultiNest papers), can be changed to "parameter" to speed things up, but may lead to unreliable results.
+
+		n_live : int - optional
+			Number of MultiNest live points. Reducing speeds up the code but may lead to unreliable results.
+
+		const_efficiency_mode : bool - optional
+			Use MultiNest's constant efficiency mode, dramatically speeds things up, but often fails to find the correct solution.
+
+		"""
 
 		print("\nBagpipes: fitting object " + self.Galaxy.ID + "\n")
 
@@ -219,7 +236,7 @@ def fit(self, verbose=False, sampling_efficiency="model", n_live=400, const_effi
 
 
 	def prior_transform(self, cube, ndim, nparam): 
-		""" Prior function for MultiNest algorithm, currently just converts unit cube to uniform prior between set units. """
+		# Prior function for MultiNest algorithm.
 
 		for i in range(self.ndim):
 
@@ -261,7 +278,7 @@ def prior_transform(self, cube, ndim, nparam):
 
 
 	def get_lnprob(self, x, ndim, nparam):
-		""" Returns the log-probability for a given model sfh and parameter vector x. """
+		# Returns the log-probability for a given model sfh and parameter vector x. 
 
 		self.get_model(x)
 
@@ -296,7 +313,7 @@ def get_lnprob(self, x, ndim, nparam):
 
 
 	def get_model(self, param):
-		""" Generates a model object for the a specified set of parameters """
+		# Generates a model object for the a specified set of parameters.
 
 		self.model_components = self.get_model_components(param)
 		"""
@@ -338,7 +355,7 @@ def get_model(self, param):
 
 
 	def get_model_components(self, param):
-		""" Turns a vector of parameters into a model_components dict, if input is already a dict simply returns it. """
+		# Turns a vector of parameters into a model_components dict, if input is already a dict simply returns it. 
 
 		# If param is a model_components dictionary get right on with calculating the chi squared value
 		if isinstance(param, dict):
@@ -366,7 +383,7 @@ def get_model_components(self, param):
 
 
 	def load_posterior(self):
-		""" Load the posterior samples generated by MultiNest. """
+		# Load the posterior samples generated by MultiNest. 
 
 		if "samples" not in list(self.posterior):
 			self.posterior["samples"] = np.loadtxt(models.working_dir + "/pipes/pmn_chains/" + self.run + "/" + self.Galaxy.ID + "-post_equal_weights.dat")[:,:-1]
@@ -376,7 +393,7 @@ def load_posterior(self):
 
 
 	def get_post_info(self):
-		""" Calculates a whole bunch of useful posterior quantities from the MultiNest output. """
+		# Calculates a whole bunch of useful posterior quantities from the MultiNest output. 
 
 		if "sfh" not in list(self.posterior):
 
@@ -451,7 +468,16 @@ def get_post_info(self):
 
 
 	def plot_fit(self, return_fig=False):
-		""" Generate a plot of the input data and fitted posterior spectrum/photometry. """
+		""" Generate a plot of the input data and fitted posterior spectrum/photometry. 
+	
+		Parameters
+		----------
+
+		return_fig : bool - optional
+			If True, returns the figure containing the fit to the user instead of saving it to the pipes/plots/run/ directory.
+
+		"""
+
 		normalisation_factor = 10**18
 
 		self.get_post_info()
@@ -462,9 +488,9 @@ def plot_fit(self, return_fig=False):
 		if self.Galaxy.photometry_exists == True:
 			naxes += 1
 
-		fig, axes = plt.subplots(naxes, figsize=(14, 4.*naxes))
+		fig, axes = plt.subplots(naxes, figsize=(12, 4.*naxes))
 
-		plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.75, hspace=None)
+		plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.25, hspace=0.25)
 
 		if naxes == 1:
 			axes = [axes]
@@ -489,7 +515,7 @@ def plot_fit(self, return_fig=False):
 			ylabel = "$\mathrm{f_{\lambda}}\ \mathrm{/\ erg\ s^{-1}\ \AA^{-1}}$"
 
 		if naxes > 1:
-			fig.text(0.08, 0.55, ylabel, rotation=90)
+			fig.text(0.08, 0.63, ylabel, rotation=90)
 
 		else:
 			ax1.set_ylabel(ylabel, rotation=90)
@@ -507,8 +533,8 @@ def plot_fit(self, return_fig=False):
 
 			#ax1.plot(self.Galaxy.spectrum[:, 0], normalisation_factor*self.Galaxy.spectrum[:, 1], color="red", zorder=10)
 
-			ax1.plot(self.Galaxy.spectrum[:, 0], normalisation_factor*self.Galaxy.spectrum[:, 1]/polynomial, color="dodgerblue", zorder=1)
-			ax1.fill_between(self.Galaxy.spectrum[:, 0], normalisation_factor*(self.Galaxy.spectrum[:, 1]/polynomial - self.Galaxy.spectrum[:, 2]), normalisation_factor*(self.Galaxy.spectrum[:, 1]/polynomial + self.Galaxy.spectrum[:, 2]), color="dodgerblue", zorder=1, alpha=0.75, linewidth=0)
+			ax1.plot(self.Galaxy.spectrum[:, 0], normalisation_factor*self.Galaxy.spectrum[:, 1]/polynomial, color="dodgerblue", zorder=1, lw=2)
+			#ax1.fill_between(self.Galaxy.spectrum[:, 0], normalisation_factor*(self.Galaxy.spectrum[:, 1]/polynomial - self.Galaxy.spectrum[:, 2]), normalisation_factor*(self.Galaxy.spectrum[:, 1]/polynomial + self.Galaxy.spectrum[:, 2]), color="dodgerblue", zorder=1, alpha=0.75, linewidth=0)
 
 			ax1.set_ylim(0, 1.1*normalisation_factor*np.max(self.Galaxy.spectrum[:, 1]/polynomial))
 		"""
@@ -535,12 +561,12 @@ def plot_fit(self, return_fig=False):
 
 			for axis in axes:
 				axis.errorbar(np.log10(self.Galaxy.photometry[:,0][self.Galaxy.photometry[:,1] != 0.]), normalisation_factor*self.Galaxy.photometry[:,1][self.Galaxy.photometry[:,1] != 0.], yerr=normalisation_factor*self.Galaxy.photometry[:,2][self.Galaxy.photometry[:,1] != 0.], lw=1.0, linestyle=" ", capsize=3, capthick=1, zorder=3, color="black")
-				axis.scatter(np.log10(self.Galaxy.photometry[:,0][self.Galaxy.photometry[:,1] != 0.]), normalisation_factor*self.Galaxy.photometry[:,1][self.Galaxy.photometry[:,1] != 0.], color="blue", s=75, zorder=4, linewidth=1, facecolor="blue", edgecolor="black")
+				axis.scatter(np.log10(self.Galaxy.photometry[:,0][self.Galaxy.photometry[:,1] != 0.]), normalisation_factor*self.Galaxy.photometry[:,1][self.Galaxy.photometry[:,1] != 0.], color="blue", s=50, zorder=4, linewidth=1, facecolor="blue", edgecolor="black")
 
 
 		# Add masked regions to plots
-		if os.path.exists(models.working_dir + "/pipes/masks/" + self.Galaxy.ID + "_mask") and self.Galaxy.spectrum_exists:
-			mask = np.loadtxt(models.working_dir + "/pipes/masks/" + self.Galaxy.ID + "_mask")
+		if os.path.exists(models.working_dir + "/pipes/object_masks/" + self.Galaxy.ID + "_mask") and self.Galaxy.spectrum_exists:
+			mask = np.loadtxt(models.working_dir + "/pipes/object_masks/" + self.Galaxy.ID + "_mask")
 
 			for j in range(self.Galaxy.no_of_spectra):
 				if len(mask.shape) == 1:
@@ -569,11 +595,11 @@ def plot_fit(self, return_fig=False):
 
 			for j in range(self.Model.photometry.shape[0]):
 				phot_1sig = self.posterior["photometry"][j,(self.posterior["photometry"][j,:] > np.percentile(self.posterior["photometry"][j,:], 16)) & (self.posterior["photometry"][j,:] < np.percentile(self.posterior["photometry"][j,:], 84))]
-				ax2.scatter(np.log10(np.zeros(phot_1sig.shape[0]) + self.Model.phot_wavs[j]), normalisation_factor*phot_1sig, color="darkorange", zorder=2, alpha=0.05, s=150, rasterized=True)
+				ax2.scatter(np.log10(np.zeros(phot_1sig.shape[0]) + self.Model.phot_wavs[j]), normalisation_factor*phot_1sig, color="darkorange", zorder=2, alpha=0.05, s=100, rasterized=True)
 
 		if self.Galaxy.spectrum_exists == True:
 			ax1.fill_between(self.Model.spectrum[:,0], normalisation_factor*np.percentile(self.posterior["spectrum"], 16, axis=1), normalisation_factor*np.percentile(self.posterior["spectrum"], 84, axis=1), color="sandybrown", zorder=2, alpha=0.75, linewidth=0)
-			ax1.plot(self.Model.spectrum[:,0], normalisation_factor*np.percentile(self.posterior["spectrum"], 50, axis=1), color="sandybrown", zorder=2)
+			ax1.plot(self.Model.spectrum[:,0], normalisation_factor*np.percentile(self.posterior["spectrum"], 50, axis=1), color="sandybrown", zorder=2, lw=1.5)
 
 			ax1.set_ylim(0, np.max([ax1.get_ylim()[1], 1.1*np.max(normalisation_factor*np.percentile(self.posterior["spectrum"], 84, axis=1))]))
 			"""
@@ -595,7 +621,7 @@ def plot_fit(self, return_fig=False):
 
 
 	def plot_poly(self):
-		""" Plot the posterior for the polynomial correction applied to the spectrum. """
+		# Plot the posterior for the polynomial correction applied to the spectrum. 
 
 		self.get_post_info()
 
@@ -609,8 +635,25 @@ def plot_poly(self):
 
 
 
-	def plot_corner(self, param_names_tolog=[], truths=None, ranges=None, return_fig=False):
-		""" Make a corner plot of the fitting parameters. """
+	def plot_corner(self, return_fig=False, truths=None, ranges=None, param_names_tolog=[]):
+		""" Make a corner plot of the fitting parameters. 
+
+		Parameters
+		----------
+
+		return_fig : bool - optional
+			If True, returns the figure containing the fit to the user instead of saving it to the pipes/plots/run/ directory.
+
+		truths : list - optional
+			List of true parameter values in the same order as Fit.fit_params, if supplied adds the true values to the corner plot.
+
+		ranges : list - optional
+			List of tuples containing upper and lower limits for the plot areas for each parameter (same order as Fit.fit_params).
+
+		param_names_tolog : list - optional
+			List of parameter names formatted as in Fit.fit_params. The log_10 of the posteriors for these parameters will be plotted.
+
+		"""
 
 		param_cols_toplot = []
 		param_names_toplot = []

bagpipes/galaxy.py

@@ -1,11 +1,13 @@
+from __future__ import print_function, division, absolute_import
+
 import numpy as np
 import sys
 import os
 
 import matplotlib.pyplot as plt
 from matplotlib import gridspec
 
-import model_manager as models 
+from . import model_manager as models 
 
 class Galaxy:
 
@@ -142,7 +144,7 @@ def plot(self):
         ax2.set_xlabel("$\mathrm{log_{10}}\\Big(\lambda / \mathrm{\AA}\\Big)$", size=18)
 
         if naxes == 2:
-            fig.text(0.06, 0.58, "$\mathrm{f_{\lambda}}\ \mathrm{/\ erg\ s^{-1}\ cm^{-2}\ \AA^{-1}}$", size=18, rotation=90)
+            fig.text(0.06, 0.58, "$\mathrm{f_{\lambda}}\ \mathrm{/\ 10^{-18}\ erg\ s^{-1}\ cm^{-2}\ \AA^{-1}}$", size=18, rotation=90)
 
         else:
             ax1.set_ylabel("$\mathrm{f_{\lambda}}\ \mathrm{/\ 10^{-18}\ erg\ s^{-1}\ cm^{-2}\ \AA^{-1}}$", size=18)
@@ -183,8 +185,8 @@ def plot(self):
 
 
         # Add masked regions to plots
-        if os.path.exists(models.working_dir + "/pipes/masks/" + self.ID + "_mask") and self.spectrum_exists:
-            mask = np.loadtxt(models.working_dir + "/pipes/masks/" + self.ID + "_mask")
+        if os.path.exists(models.working_dir + "/pipes/object_masks/" + self.ID + "_mask") and self.spectrum_exists:
+            mask = np.loadtxt(models.working_dir + "/pipes/object_masks/" + self.ID + "_mask")
 
             for j in range(self.no_of_spectra):
                 if len(mask.shape) == 1:
@@ -193,9 +195,7 @@ def plot(self):
                 if len(mask.shape) == 2:
                     for i in range(mask.shape[0]):
                         axes[j].axvspan(mask[i,0], mask[i,1], color="gray", alpha=0.8, zorder=3)
-
-        plt.savefig("/Users/adam/using_bagpipes/JWST_targets/plots/" + self.ID + "_3dhst_spec.pdf", bbox_inches="tight")
-
+
         plt.show()
         plt.close(fig)
 
@@ -214,11 +214,11 @@ def get_eff_wavs(self):
 
     """ Set the error spectrum to infinity in masked regions. """
     def mask_spectrum(self, spectrum):
-        if not os.path.exists(models.install_dir + "/object_masks/" + self.ID + "_mask"): #" + self.ID + "
+        if not os.path.exists(models.working_dir + "/pipes/object_masks/" + self.ID + "_mask"): #" + self.ID + "
             return spectrum
 
         else:
-            mask = np.loadtxt(models.install_dir + "/object_masks/" + self.ID + "_mask") #" + self.ID + "
+            mask = np.loadtxt(models.working_dir + "/pipes/object_masks/" + self.ID + "_mask") #" + self.ID + "
             if len(mask.shape) == 1:
                 if spectrum[(spectrum[:,0] > mask[0]) & (spectrum[:,0] < mask[1]), 2].shape[0] is not 0:
                     spectrum[(spectrum[:,0] > mask[0]) & (spectrum[:,0] < mask[1]), 2] = 9.9*10**99.