v0.7.0

Fixes or introduces a number of features that needed to be resolved before
trying to move to the initial production run. These include the following
changes:

- Resolves #31 by introducing a new `title_quantiles` argument. Also modifies
  the default plotting options so now it just shows bins.
- Removed the binned LOS log-likelihood, since this has now been superseeded
  by the sample-based version. Also added in an option to treat the foreground
  as an additive component, rather than just a separate number.
- Modifications to some of the internals in the `seds` module to make things
  easier to pull out. Should resolve #28.
- Cleaned up imports and similar things so the code should now be
  more self-contained.

Next up is working on the clusters module so I can fit some isochrones!

Author: joshspeagle

README.md

@@ -25,23 +25,26 @@ stellar models, it is configured for two by default:
 and [Bayestar](https://arxiv.org/pdf/1401.1508.pdf).
 
 #### Zero-points
-Zero-point offsets in several bands have been derived using Gaia data
-and can be included during runtime.
+Zero-point offsets in several bands have been estimated using Gaia data
+and can be included during runtime. 
+**These are currently not thoroughly vetted, so use at your own risk.**
 
 #### Dust Map
 `brutus` is able to incorporate a 3-D dust prior. The current prior is
-based on the "Bayestar17" dust map from
-[Green et al. (2018)](https://arxiv.org/abs/1801.03555).
+based on the "Bayestar19" dust map from
+[Green et al. (2019)](https://arxiv.org/abs/1905.02734).
 
 #### Generating SEDs
 `brutus` contains built-in SED generation utilities based on the MIST
 stellar models, modeled off of
 [`minesweeper`](https://github.com/pacargile/MINESweeper).
 These are optimized for either generating photometry from stellar mass
-tracks or for a single-age stellar isochrone, and are based on
+tracks or for a single-age stellar isochrone based on
 artificial neural networks trained on bolometric correction tables.
+
 An empirical correction table to the models derived using several clusters is
 also provided, which improves the models down to ~0.5 solar masses.
+**These are currently not thoroughly vetted, so use at your own risk.**
 
 Please contact Phil Cargile (pcargile@cfa.harvard.edu) and Josh Speagle
 (jspeagle@cfa.harvard.edu) for more information on the provided data files.

brutus/__init__.py

@@ -2,9 +2,8 @@
 # -*- coding: utf-8 -*-
 
 from __future__ import (division, print_function)
-from six.moves import range
 
 from .fitting import *
 from .utils import *
 
-__version__ = "0.6.9"
+__version__ = "0.7.0"

brutus/cluster.py

@@ -7,16 +7,5 @@
 """
 
 from __future__ import (print_function, division)
-import six
-from six.moves import range
-
-import sys
-import os
-import warnings
-import math
-import numpy as np
-import warnings
-
-from .utils import get_seds, add_mag
 
 __all__ = [""]

brutus/dust.py

@@ -8,15 +8,8 @@
 """
 
 from __future__ import (print_function, division)
-import six
-from six.moves import range
 
-import sys
-import os
-import warnings
-import math
 import numpy as np
-import warnings
 import h5py
 
 import astropy.coordinates as coordinates
@@ -287,9 +280,6 @@ def query(self, coords):
 
         """
 
-        # Get number of coordinates requested.
-        n_coords_ret = coords.shape[0]
-
         # Extract the correct angular pixel(s).
         try:
             pix_idx = self._find_data_idx(coords.l.deg, coords.b.deg)

brutus/filters.py

@@ -7,8 +7,6 @@
 """
 
 from __future__ import (print_function, division)
-import six
-from six.moves import range
 
 # Define the set of filters available for the provided MIST models.
 # The Bayestar models are only defined for PanSTARRS `grizy` and 2MASS.

brutus/fitting.py

@@ -7,14 +7,10 @@
 """
 
 from __future__ import (print_function, division)
-from six.moves import range
 
 import sys
-import os
 import warnings
-import math
 import numpy as np
-import warnings
 import h5py
 import time
 from scipy.special import xlogy, gammaln
@@ -24,8 +20,10 @@
 except ImportError:
     from scipy.misc import logsumexp
 
-from .pdf import *
-from .utils import *
+from .pdf import imf_lnprior, ps1_MrLF_lnprior
+from .pdf import parallax_lnprior, scale_parallax_lnprior
+from .pdf import gal_lnprior, dust_lnprior
+from .utils import _function_wrapper, _inverse3, magnitude, get_seds
 
 __all__ = ["loglike", "_optimize_fit", "BruteForce", "_lnpost"]
 
@@ -608,9 +606,9 @@ def __init__(self, models, models_labels, labels_mask, pool=None):
 
     def fit(self, data, data_err, data_mask, data_labels, save_file,
             phot_offsets=None, parallax=None, parallax_err=None,
-            Nmc_prior=100, avlim=(0., 6.), av_gauss=None,
+            Nmc_prior=50, avlim=(0., 6.), av_gauss=None,
             rvlim=(1., 8.), rv_gauss=(3.32, 0.18), binary_frac=0.5,
-            lnprior=None, wt_thresh=1e-3, cdf_thresh=2e-4, Ndraws=2000,
+            lnprior=None, wt_thresh=1e-3, cdf_thresh=2e-4, Ndraws=500,
             apply_agewt=True, apply_grad=True,
             lndistprior=None, lndustprior=None, dustfile='bayestar2017_v1.h5',
             apply_dlabels=True, data_coords=None, logl_dim_prior=True,
@@ -837,9 +835,9 @@ def fit(self, data, data_err, data_mask, data_labels, save_file,
             try:
                 # Try reading in parallel-friendly way if possible.
                 try:
-                    ft = h5py.File(dustfile, 'r', libver='latest', swmr=True)
+                    h5py.File(dustfile, 'r', libver='latest', swmr=True)
                 except:
-                    ft = h5py.File(dustfile, 'r')
+                    h5py.File(dustfile, 'r')
                     pass
             except:
                 raise ValueError("The default dust prior is being used but "
@@ -1152,13 +1150,12 @@ def _fit(self, data, data_err, data_mask,
 
         """
 
-        Ndata, Nmodels = len(data), self.NMODEL
+        Ndata = len(data)
         models = np.array(self.models)
         if wt_thresh is None and cdf_thresh is None:
             wt_thresh = -np.inf  # default to no clipping/thresholding
         if rstate is None:
             rstate = np.random
-        mvn = rstate.multivariate_normal
         if parallax is not None and parallax_err is None:
             raise ValueError("Must provide both `parallax` and "
                              "`parallax_err`.")
@@ -1189,9 +1186,9 @@ def _fit(self, data, data_err, data_mask,
             try:
                 # Try reading in parallel-friendly way if possible.
                 try:
-                    ft = h5py.File(dustfile, 'r', libver='latest', swmr=True)
+                    h5py.File(dustfile, 'r', libver='latest', swmr=True)
                 except:
-                    ft = h5py.File(dustfile, 'r')
+                    h5py.File(dustfile, 'r')
                     pass
             except:
                 raise ValueError("The default dust prior is being used but "

brutus/los.py

@@ -7,24 +7,17 @@
 """
 
 from __future__ import (print_function, division)
-import six
-from six.moves import range
 
-import sys
-import os
 import warnings
-import math
 import numpy as np
-import warnings
 from scipy.stats import truncnorm
 
 try:
     from scipy.special import logsumexp
 except ImportError:
     from scipy.misc import logsumexp
 
-__all__ = ["LOS_clouds_priortransform", "LOS_clouds_loglike_bin",
-           "LOS_clouds_loglike_samples",
+__all__ = ["LOS_clouds_priortransform", "LOS_clouds_loglike_samples",
            "kernel_tophat", "kernel_gauss", "kernel_lorentz"]
 
 
@@ -125,134 +118,9 @@ def LOS_clouds_priortransform(u, rlims=(0., 6.), dlims=(4., 19.),
     return x
 
 
-def LOS_clouds_loglike_bin(theta, cdfs, xedges, yedges, interpolate=True):
-    """
-    Compute the log-likelihood for the cumulative reddening along the
-    line of sight (LOS) parameterized by `theta`, given input binned
-    stellar posteriors/bounds generated by `.pdf.bin_pdfs_distred`. Assumes
-    a uniform outlier model in distance and reddening across our binned
-    posteriors.
-
-    Parameters
-    ----------
-    theta : `~numpy.ndarray` of shape `(Nparams,)`
-        A collection of parameters that characterizes the cumulative
-        reddening along the LOS. Contains the fraction of outliers `P_b`
-        followed by the foreground reddening `fred` and a series of
-        `(dist, red)` pairs for each "cloud" along the LOS.
-
-    cdfs : `~numpy.ndarray` of shape `(Nobj, Nxbin, Nybin)`
-        Binned versions of the CDFs.
-
-    xedges : `~numpy.ndarray` of shape `(Nxbin+1,)`
-        The edges defining the bins in distance.
-
-    yedges : `~numpy.ndarray` of shape `(Nybin+1,)`
-        The edges defining the bins in reddening.
-
-    interpolate : bool, optional
-        Whether to linearly interpolate between bins (defined by their central
-        positions). Default is `True`.
-
-    Returns
-    -------
-    loglike : float
-        The computed log-likelihood.
-
-    """
-
-    # Grab parameters.
-    pb = theta[0]
-    reds, dists = np.atleast_1d(theta[1::2]), np.atleast_1d(theta[2::2])
-
-    # Convert to bin coordinates.
-    dx, dy = xedges[1] - xedges[0], yedges[1] - yedges[0]
-    Nxedges, Nyedges = len(xedges), len(yedges)
-    Nxbin, Nybin = Nxedges - 1, Nyedges - 1
-    x_ctrs = np.arange(0.5, Nxbin, 1.)
-    y_ctrs = np.arange(0.5, Nybin, 1.)
-    x = np.concatenate(([0], (dists - xedges[0]) / dx, [Nxbin]))
-    y = (reds - yedges[0]) / dy
-
-    # Find (x,y) neighbors in bin space.
-    x1 = np.array(np.floor(x - 0.5), dtype='int')
-    x2 = np.array(np.ceil(x - 0.5), dtype='int')
-    y1 = np.array(np.floor(y - 0.5), dtype='int')
-    y2 = np.array(np.ceil(y - 0.5), dtype='int')
-
-    # Threshold values to bin edges.
-    x1[np.where(x1 < 0)] = 0
-    x1[np.where(x1 >= Nxbin)] = Nxbin - 1
-    x2[np.where(x2 < 0)] = 0
-    x2[np.where(x2 >= Nxbin)] = Nxbin - 1
-    y1[np.where(y1 < 0)] = 0
-    y1[np.where(y1 >= Nybin)] = Nybin - 1
-    y2[np.where(y2 < 0)] = 0
-    y2[np.where(y2 >= Nybin)] = Nybin - 1
-
-    # Threshold (x,y) to edges (and shift to deal with centers).
-    x[np.where(x < 0.5)] = 0.5
-    x[np.where(x > Nxbin - 0.5)] = Nxbin - 0.5
-    y[np.where(y < 0.5)] = 0.5
-    y[np.where(y > Nybin - 0.5)] = Nybin - 0.5
-
-    # Define "left" and "right" versions for xs.
-    x1l, x1r = x1[:-1], x1[1:]
-    x2l, x2r = x2[:-1], x2[1:]
-    xl, xr = x[:-1], x[1:]
-
-    # Compute integral along LOS using the provided CDFs.
-    if interpolate:
-        # Interpolate between bins (left side).
-        # Define values q(x_i, y_i).
-        q11, q12, q21, q22 = (cdfs[:, x1l, y1], cdfs[:, x1l, y2],
-                              cdfs[:, x2l, y1], cdfs[:, x2l, y2])
-        # Compute areas.
-        dx1, dx2 = (xl - x_ctrs[x1l]), (x_ctrs[x2l] - xl)
-        dy1, dy2 = (y - y_ctrs[y1]), (y_ctrs[y2] - y)
-        # Interpolate in x.
-        qp1, qp2 = (q11 * dx2 + q21 * dx1), (q12 * dx2 + q22 * dx1)
-        xsel = np.where(~((dx1 > 0.) & (dx2 > 0.)))  # deal with edges
-        qp1[:, xsel], qp2[:, xsel] = q11[:, xsel], q12[:, xsel]  # replace
-        # Interpolate in y.
-        cdf_left = qp1 * dy2 + qp2 * dy1
-        ysel = np.where(~((dy1 > 0.) & (dy2 > 0.)))  # deal with edges
-        cdf_left[ysel] = qp1[ysel]  # replace edges
-
-        # Interpolate between the bins (right side).
-        # Define values q(x_i, y_i).
-        q11, q12, q21, q22 = (cdfs[:, x1r, y1], cdfs[:, x1r, y2],
-                              cdfs[:, x2r, y1], cdfs[:, x2r, y2])
-        # Compute areas.
-        dx1, dx2 = (xr - x_ctrs[x1r]), (x_ctrs[x2r] - xr)
-        dy1, dy2 = (y - y_ctrs[y1]), (y_ctrs[y2] - y)
-        # Interpolate in x.
-        qp1, qp2 = (q11 * dx2 + q21 * dx1), (q12 * dx2 + q22 * dx1)
-        xsel = np.where(~((dx1 > 0.) & (dx2 > 0.)))  # deal with edges
-        qp1[:, xsel], qp2[:, xsel] = q11[:, xsel], q12[:, xsel]  # replace
-        # Interpolate in y.
-        cdf_right = qp1 * dy2 + qp2 * dy1
-        ysel = np.where(~((dy1 > 0.) & (dy2 > 0.)))  # deal with edges
-        cdf_right[ysel] = qp1[ysel]  # replace edges
-    else:
-        # Just use the values from the bin we're currently in.
-        cdf_left, cdf_right = cdfs[:, x1l, y1], cdfs[:, x1r, y1]
-
-    # Compute likelihood.
-    likes = np.sum(cdf_right - cdf_left, axis=1)
-
-    # Add in outlier mixture model. Assume uniform in (x, y) with `pb` weight.
-    likes = pb * (1. / (Nybin * Nxbin)) + (1. - pb) * likes
-
-    # Compute total log-likelihood.
-    loglike = np.sum(np.log(likes))
-
-    return loglike
-
-
 def LOS_clouds_loglike_samples(theta, dsamps, rsamps, kernel='gauss',
                                rlims=(0., 6.), template_reds=None,
-                               Ndraws=25):
+                               Ndraws=25, additive_foreground=False):
     """
     Compute the log-likelihood for the cumulative reddening along the
     line of sight (LOS) parameterized by `theta`, given a set of input
@@ -293,6 +161,10 @@ def LOS_clouds_loglike_samples(theta, dsamps, rsamps, kernel='gauss',
     Ndraws : int, optional
         The number of draws to use for each star. Default is `25`.
 
+    additive_foreground : bool, optional
+        Whether the foreground is treated as just another value or added
+        to all background values. Default is `False`.
+
     Returns
     -------
     loglike : float
@@ -336,6 +208,10 @@ def LOS_clouds_loglike_samples(theta, dsamps, rsamps, kernel='gauss',
     if template_reds is not None:
         reds[1:] *= template_reds[None, :, None]  # reds[1:] are rescalings
 
+    # Adjust reddenings after the foreground if needed.
+    if additive_foreground:
+        reds[1:] += reds[0]  # add foreground to background
+
     # Define kernel parameters (mean, sigma) per LOS chunk.
     kparams = np.array([(r, rsmooth) for r in reds])
     kparams[0][1] = rsmooth0