Skip to content

Commit

Permalink
v0.7.0
Browse files Browse the repository at this point in the history 
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!
  • Loading branch information
@joshspeagle
joshspeagle committed Jul 20, 2019
1 parent 09e37d0 commit 40a0bc3
Show file tree
Hide file tree
Showing 11 changed files with 69 additions and 225 deletions.
13 changes: 8 additions & 5 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
3 changes: 1 addition & 2 deletions brutus/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
11 changes: 0 additions & 11 deletions brutus/cluster.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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__ = [""]
10 changes: 0 additions & 10 deletions brutus/dust.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 0 additions & 2 deletions brutus/filters.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
25 changes: 11 additions & 14 deletions brutus/fitting.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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"]

Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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 "
Expand Down Expand Up @@ -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`.")
Expand Down Expand Up @@ -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 "
Expand Down
144 changes: 10 additions & 134 deletions brutus/los.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"]


Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit 40a0bc3

Please sign in to comment.