Merge branch 'master' into retrieval

species/core/box.py

@@ -1,14 +1,22 @@
 """
-Box module.
+Module with classes for storing data in `Box` objects.
 """
 
+import numpy as np
+import spectres
+
+import species
+
 
 def create_box(boxtype,
                **kwargs):
     """
+    Function for creating a :class:`species.core.box.Box`.
+
     Returns
     -------
     species.core.box
+        Box with the data and parameters.
     """
 
     if boxtype == 'colormag':
@@ -119,7 +127,7 @@ def create_box(boxtype,
 
 class Box:
     """
-    Text
+    Class for generic methods that can be applied on all `Box` object.
     """
 
     def __init__(self):
@@ -132,6 +140,8 @@ def __init__(self):
 
     def open_box(self):
         """
+        Method for inspecting the content of a `Box`.
+
         Returns
         -------
         NoneType
@@ -146,7 +156,7 @@ def open_box(self):
 
 class ColorMagBox(Box):
     """
-    Text
+    Class for storing color-magnitude data in a :class:`species.core.box.Box`.
     """
 
     def __init__(self):
@@ -169,7 +179,7 @@ def __init__(self):
 
 class ColorColorBox(Box):
     """
-    Text
+    Class for storing color-color data in a :class:`species.core.box.Box`.
     """
 
     def __init__(self):
@@ -191,7 +201,7 @@ def __init__(self):
 
 class IsochroneBox(Box):
     """
-    Text
+    Class for storing isochrone data in a :class:`species.core.box.Box`.
     """
 
     def __init__(self):
@@ -214,7 +224,7 @@ def __init__(self):
 
 class ModelBox(Box):
     """
-    Text
+    Class for storing a model spectrum in a :class:`species.core.box.Box`.
     """
 
     def __init__(self):
@@ -233,10 +243,57 @@ def __init__(self):
         self.quantity = None
         self.contribution = None
 
+    def smooth_spectrum(self,
+                        spec_res: float) -> None:
+        """
+        Method for smoothing the spectrum with a Gaussian kernel to the instrument resolution. The
+        method is best applied on an input spectrum with a logarithmic wavelength sampling (i.e.
+        constant spectral resolution). Alternatively, the wavelength sampling may be linear, but
+        the smoothing is slower in that case.
+
+        Parameters
+        ----------
+        spec_res : float
+            Spectral resolution that is used for the smoothing kernel.
+
+        Returns
+        -------
+        NoneType
+            None
+        """
+
+        self.flux = species.util.read_util.smooth_spectrum(
+            self.wavelength, self.flux, spec_res)
+
+    def resample_spectrum(self,
+                          wavel_resample: np.ndarray) -> None:
+        """
+        Method for resampling the spectrum with ``spectres`` to a new wavelength grid.
+
+        Parameters
+        ----------
+        wavel_resample : np.ndarray
+            Wavelength points (um) to which the spectrum will be resampled.
+
+        Returns
+        -------
+        NoneType
+            None
+        """
+
+        self.flux = spectres.spectres(wavel_resample,
+                                      self.wavelength,
+                                      self.flux,
+                                      spec_errs=None,
+                                      fill=np.nan,
+                                      verbose=True)
+
+        self.wavelength = wavel_resample
+
 
 class ObjectBox(Box):
     """
-    Text
+    Class for storing object data in a :class:`species.core.box.Box`.
     """
 
     def __init__(self):
@@ -257,7 +314,7 @@ def __init__(self):
 
 class PhotometryBox(Box):
     """
-    Text
+    Class for storing photometric data in a :class:`species.core.box.Box`.
     """
 
     def __init__(self):
@@ -279,7 +336,7 @@ def __init__(self):
 
 class ResidualsBox(Box):
     """
-    Text
+    Class for storing best-fit residuals in a :class:`species.core.box.Box`.
     """
 
     def __init__(self):
@@ -297,7 +354,7 @@ def __init__(self):
 
 class SamplesBox(Box):
     """
-    Text
+    Class for storing posterior samples in a :class:`species.core.box.Box`.
     """
 
     def __init__(self):
@@ -319,7 +376,7 @@ def __init__(self):
 
 class SpectrumBox(Box):
     """
-    Text
+    Class for storing spectral data in a :class:`species.core.box.Box`.
     """
 
     def __init__(self):
@@ -342,7 +399,7 @@ def __init__(self):
 
 class SynphotBox(Box):
     """
-    Text
+    Class for storing synthetic photometry in a :class:`species.core.box.Box`.
     """
 
     def __init__(self):

species/data/blackbody.py

@@ -60,11 +60,11 @@ def add_blackbody(input_path: str,
     url = 'https://home.strw.leidenuniv.nl/~stolker/species/blackbody.tgz'
 
     if not os.path.isfile(data_file):
-        print('Downloading blackbody model spectra (46 MB)...', end='', flush=True)
+        print('Downloading blackbody model spectra (56 MB)...', end='', flush=True)
         urllib.request.urlretrieve(url, data_file)
         print(' [DONE]')
 
-    print('Unpacking blackbody model spectra (46 MB)...', end='', flush=True)
+    print('Unpacking blackbody model spectra (56 MB)...', end='', flush=True)
     tar = tarfile.open(data_file)
     tar.extractall(data_folder)
     tar.close()
@@ -90,7 +90,7 @@ def add_blackbody(input_path: str,
                         continue
 
                 print_message = f'Adding blackbody model spectra... {filename}'
-                print(f'\r{print_message:<62}', end='')
+                print(f'\r{print_message:<63}', end='')
 
                 data_wavel, data_flux = np.loadtxt(os.path.join(data_folder, filename), unpack=True)
 
@@ -123,7 +123,7 @@ def add_blackbody(input_path: str,
                     flux.append(flux_resample)  # (W m-2 um-1)
 
     print_message = 'Adding blackbody model spectra... [DONE]'
-    print(f'\r{print_message:<62}')
+    print(f'\r{print_message:<63}')
 
     data_sorted = data_util.sort_data(np.asarray(teff),
                                       None,

species/data/database.py

@@ -1393,27 +1393,34 @@ def get_mcmc_spectra(self,
                          random: int,
                          burnin: Optional[int] = None,
                          wavel_range: Optional[Union[Tuple[float, float], str]] = None,
-                         spec_res: Optional[float] = None) -> Union[List[box.ModelBox],
-                                                                    List[box.SpectrumBox]]:
+                         spec_res: Optional[float] = None,
+                         wavel_resample: Optional[np.ndarray] = None) -> \
+                             Union[List[box.ModelBox], List[box.SpectrumBox]]:
         """
+        Function for drawing random spectra from the sampled posterior distributions.
+
         Parameters
         ----------
         tag : str
             Database tag with the posterior samples.
         random : int
             Number of random samples.
         burnin : int, None
-            Number of burnin steps. No burnin is removed if set to ``None``.
+            Number of burnin steps. No burnin is removed if set to ``None``. Not required when
+            using nested sampling.
         wavel_range : tuple(float, float), str, None
             Wavelength range (um) or filter name. Full spectrum is used if set to ``None``.
         spec_res : float, None
             Spectral resolution that is used for the smoothing with a Gaussian kernel. No smoothing
-            is applied if set to ``None``.
+            is applied if the argument set to ``None``.
+        wavel_resample : np.ndarray, None
+            Wavelength points (um) to which the model spectrum will be resampled. The resampling is
+            applied after the optional smoothing to the resolution of ``spec_res``.
 
         Returns
         -------
         list(species.core.box.ModelBox)
-            Boxes with the randomly sampled spectra.
+            List with ``ModelBox`` objects.
         """
 
         if burnin is None:
@@ -1515,13 +1522,23 @@ def get_mcmc_spectra(self,
 
             if spectrum_type == 'model':
                 if spectrum_name == 'planck':
-                    specbox = readmodel.get_spectrum(model_param, spec_res)
+                    specbox = readmodel.get_spectrum(model_param,
+                                                     spec_res,
+                                                     smooth=True,
+                                                     wavel_resample=wavel_resample)
 
                 elif spectrum_name == 'powerlaw':
+                    if wavel_resample is not None:
+                        warnings.warn('The \'wavel_resample\' parameter is not support by the '
+                                      '\'powerlaw\' model so the argument will be ignored.')
+
                     specbox = read_util.powerlaw_spectrum(wavel_range, model_param)
 
                 else:
-                    specbox = readmodel.get_model(model_param, spec_res=spec_res, smooth=True)
+                    specbox = readmodel.get_model(model_param,
+                                                  spec_res=spec_res,
+                                                  wavel_resample=wavel_resample,
+                                                  smooth=True)
 
             elif spectrum_type == 'calibration':
                 specbox = readcalib.get_spectrum(model_param)
@@ -1536,23 +1553,33 @@ def get_mcmc_spectra(self,
     def get_mcmc_photometry(self,
                             tag: str,
                             filter_name: str,
-                            burnin: Optional[int] = None) -> np.ndarray:
+                            burnin: Optional[int] = None,
+                            phot_type: str = 'magnitude') -> np.ndarray:
         """
+        Function for calculating synthetic magnitudes or fluxes from the posterior samples.
+
         Parameters
         ----------
         tag : str
             Database tag with the posterior samples.
         filter_name : str
-            Filter name for which the photometry will be computed.
+            Filter name for which the synthetic photometry will be computed.
         burnin : int, None
-            Number of burnin steps. No burnin is removed if set to ``None``.
+            Number of burnin steps. No burnin is removed if set to ``None``. Not required when
+            using nested sampling.
+        phot_type : str
+            Photometry type ('magnitude' or 'flux').
 
         Returns
         -------
         np.ndarray
-            Synthetic magnitudes.
+            Synthetic magnitudes or fluxes (W m-2 um-1).
         """
 
+        if phot_type not in ['magnitude', 'flux']:
+            raise ValueError('The argument of \'phot_type\' is not recognized and should be '
+                             'set to \'magnitude\' or \'flux\'.')
+
         if burnin is None:
             burnin = 0
 
@@ -1603,6 +1630,7 @@ def get_mcmc_photometry(self,
 
         for i in tqdm.tqdm(range(samples.shape[0]), desc='Getting MCMC photometry'):
             model_param = {}
+
             for j in range(n_param):
                 model_param[param[j]] = samples[i, j]
 
@@ -1613,16 +1641,27 @@ def get_mcmc_photometry(self,
                 if spectrum_name == 'powerlaw':
                     pl_box = read_util.powerlaw_spectrum(synphot.wavel_range, model_param)
 
-                    app_mag, _ = synphot.spectrum_to_magnitude(pl_box.wavelength, pl_box.flux)
+                    if phot_type == 'magnitude':
+                        app_mag, _ = synphot.spectrum_to_magnitude(pl_box.wavelength, pl_box.flux)
+                        mcmc_phot[i] = app_mag[0]
 
-                    mcmc_phot[i] = app_mag[0]
+                    elif phot_type == 'flux':
+                        mcmc_phot[i], _ = synphot.spectrum_to_flux(pl_box.wavelength, pl_box.flux)
 
                 else:
-                    mcmc_phot[i], _ = readmodel.get_magnitude(model_param)
+                    if phot_type == 'magnitude':
+                        mcmc_phot[i], _ = readmodel.get_magnitude(model_param)
+
+                    elif phot_type == 'flux':
+                        mcmc_phot[i], _ = readmodel.get_flux(model_param)
 
             elif spectrum_type == 'calibration':
-                app_mag, _ = readcalib.get_magnitude(model_param=model_param, distance=None)
-                mcmc_phot[i] = app_mag[0]
+                if phot_type == 'magnitude':
+                    app_mag, _ = readcalib.get_magnitude(model_param=model_param, distance=None)
+                    mcmc_phot[i] = app_mag[0]
+
+                elif phot_type == 'flux':
+                    mcmc_phot[i], _ = readcalib.get_flux(model_param=model_param)
 
         return mcmc_phot
 

species/read/read_calibration.py

@@ -263,8 +263,12 @@ def get_flux(self,
 
         Returns
         -------
-        tuple(float, float)
-            Average flux and uncertainty (W m-2 um-1).
+        Returns
+        -------
+        float
+            Average flux (W m-2 um-1).
+        float
+            Uncertainty (W m-2 um-1).
         """
 
         specbox = self.get_spectrum(model_param=model_param)