In [ ]:
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
import h5py
In [34]:
from brutus import filters
from brutus import seds
from brutus import utils as butils
filt = filters.decam[1:] + filters.tmass
# import EEP tracks
nnfile = "brutus/data/DATAFILES/nn_c3k.h5"
mistfile = "brutus/data/DATAFILES/MIST_1.2_EEPtrk.h5"
sedmaker = seds.SEDmaker(filters=filt, nnfile=nnfile, mistfile=mistfile)
# import MIST model grid
gridfile = "brutus/data/DATAFILES/grid_mist_v9.h5"
(models_mist, labels_mist, lmask_mist) = butils.load_models(gridfile, filters=filt)
Filters: ['DECam_g', 'DECam_r', 'DECam_i', 'DECam_z', 'DECam_Y', '2MASS_J', '2MASS_H', '2MASS_Ks'] Constructing MIST library...done! Computing age weights for track (mini, feh, afe) = (300.0, 0.5, 0.0) Initializing FastNN predictor...done! Reading filter 2MASS_Ks
In [38]:
# load offsets
off_file_mist = "brutus/data/DATAFILES/offsets_mist_v8.txt"
print("MIST photometric offsets:")
off_mist = butils.load_offsets(off_file_mist, filters=filt)
MIST photometric offsets:
DECam_g (1.0%) DECam_r (-3.0%) DECam_i (-3.0%) DECam_z (-2.0%) DECam_Y (-3.0%) 2MASS_J (-2.0%) 2MASS_H (1.0%) 2MASS_Ks (-2.0%)
In [39]:
from brutus.fitting import BruteForce
# initialize `BruteForce` class
BF_mist = BruteForce(models_mist, labels_mist, lmask_mist)
In [68]:
# df = pd.read_csv("dust_region_stars_photometry.csv")
df = pd.read_csv("dust_region_stars_photometry.csv")
# replace 0.0 with NaN
df.replace(0.0, np.nan, inplace=True)
df.drop(14)
df
Out[68]:
| ra | dec | l | b | g_maggies | g_err_maggies | r_maggies | r_err_maggies | i_maggies | i_err_maggies | ... | y_maggies | y_err_maggies | j_maggies | j_err_maggies | h_maggies | h_err_maggies | k_maggies | k_err_maggies | parallax_mas | parallax_err_mas | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 158.595688 | -60.640267 | 286.897020 | -2.108413 | 3.262650e-10 | 5.512280e-11 | 1.486350e-09 | 1.103840e-10 | 1.859870e-09 | 9.591350e-11 | ... | 4.244120e-09 | 5.658440e-10 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 1 | 158.593768 | -60.637999 | 286.895067 | -2.106923 | 4.774170e-10 | 3.816080e-11 | 1.518170e-09 | 1.156390e-10 | 3.709260e-09 | 9.456880e-11 | ... | 8.300150e-09 | 3.654590e-10 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 2 | 158.590155 | -60.637249 | 286.893156 | -2.107162 | 4.326790e-10 | 6.233630e-11 | 3.272280e-09 | 1.195550e-10 | 1.479240e-08 | 1.185260e-10 | ... | 7.151570e-08 | 5.488640e-10 | 4.751164e-07 | 3.413267e-08 | 0.000002 | 6.375882e-08 | 0.000003 | 1.270950e-07 | -0.089649 | 0.624679 |
| 3 | 158.591477 | -60.640051 | 286.895124 | -2.109261 | 2.379470e-09 | 2.970760e-11 | 9.653800e-09 | 1.606340e-10 | 1.964540e-08 | 1.326030e-10 | ... | 3.027720e-08 | 4.184930e-10 | NaN | NaN | NaN | NaN | NaN | NaN | 0.549278 | 0.356951 |
| 4 | 158.603323 | -60.652291 | 286.906293 | -2.116940 | 2.010640e-09 | 3.450460e-11 | 5.733520e-09 | 1.626940e-10 | 7.139540e-09 | 1.070160e-10 | ... | 1.172830e-08 | 4.794800e-10 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 185 | 158.577578 | -60.675798 | 286.907169 | -2.143603 | 3.428220e-10 | 3.664220e-11 | 1.289940e-09 | 5.706480e-11 | 3.144120e-09 | 8.704980e-11 | ... | 5.605290e-09 | 4.503020e-10 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 186 | 158.579003 | -60.676241 | 286.907995 | -2.143636 | 1.431000e-08 | 6.734040e-11 | 3.959220e-08 | 1.771800e-10 | 5.490650e-08 | 2.286510e-10 | ... | 7.378240e-08 | 5.629940e-10 | NaN | NaN | NaN | NaN | NaN | NaN | 0.547826 | 0.139770 |
| 187 | 158.576333 | -60.676827 | 286.907158 | -2.144799 | 4.302820e-10 | 3.753040e-11 | 1.594290e-09 | 5.666400e-11 | 3.113620e-09 | 8.523950e-11 | ... | 5.999620e-09 | 3.966240e-10 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 188 | 158.583188 | -60.679296 | 286.911303 | -2.145250 | 2.556870e-09 | 3.030260e-11 | 8.227790e-09 | 7.297320e-11 | 1.425060e-08 | 1.083070e-10 | ... | 2.412240e-08 | 4.076210e-10 | NaN | NaN | NaN | NaN | NaN | NaN | 0.570677 | 0.428726 |
| 189 | 158.580621 | -60.679310 | 286.910222 | -2.145893 | 8.697740e-10 | 2.652580e-11 | 3.727510e-09 | 6.123480e-11 | 7.706750e-09 | 9.265820e-11 | ... | 1.503800e-08 | 3.841620e-10 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
190 rows × 22 columns
In [ ]:
# turn all phtommetry into an numpy array
phot = df[["g_maggies", "r_maggies", "i_maggies", "z_maggies", "y_maggies", "j_maggies", "h_maggies", "k_maggies"]].to_numpy()
err = df[["g_err_maggies", "r_err_maggies", "i_err_maggies", "z_err_maggies", "y_err_maggies", "j_err_maggies", "h_err_maggies", "k_err_maggies"]].to_numpy()
mask_band = ~np.isnan(phot)
coords = df[["l", "b"]].to_numpy()
parallax = df["parallax_mas"].to_numpy()
parallax_err = df["parallax_err_mas"].to_numpy()
In [ ]:
# fit a set of hypothetical objects
BF_mist.fit(
phot, # fluxes (in maggies)
err, # errors (in maggies)
mask_band, # band mask (True/False whether band was observed)
np.arange(len(phot)), # additional information to be written out
"dust_region_brutus_mist", # filename where results are stored (.h5 automatically added)
data_coords=coords, # array of (l, b) coordinates for Galactic prior
parallax=parallax,
parallax_err=parallax_err, # parallax measurements (in mas)
phot_offsets=off_mist, # photometric offsets applied to **data**
# dustfile=dustfile, # 3-D dustmap prior
Ndraws=500, # number of samples to save to disk
Nmc_prior=50, # number of Monte Carlo draws used to incorporate priors
logl_dim_prior=True, # use chi2 distribution instead of Gaussian
save_dar_draws=True, # save (dist, Av, Rv) samples
running_io=True, # write out objects as soon as they finish
verbose=True,
) # print progress
In [5]:
# now open dust_region_brutus_mist.h5 to examine results
hf = h5py.File("dust_region_brutus_mist.h5", "r")
print(list(hf.keys()))
['labels', 'ml_av', 'ml_cov_sar', 'ml_rv', 'ml_scale', 'model_idx', 'obj_Nbands', 'obj_chi2min', 'obj_log_evid', 'obj_log_post', 'samps_dist', 'samps_dred', 'samps_logp', 'samps_red']
In [6]:
# load results
f = h5py.File("dust_region_brutus_mist.h5", "r")
idxs_mist = f["model_idx"][:] # model indices
chi2_mist = f["obj_chi2min"][:] # best-fit chi2
nbands_mist = f["obj_Nbands"][:] # number of bands in fit
dists_mist = f["samps_dist"][:] # distance samples
reds_mist = f["samps_red"][:] # A(V) samples
dreds_mist = f["samps_dred"][:] # R(V) samples
lnps_mist = f["samps_logp"][:] # log-posterior of samples
In [ ]:
from brutus import plotting as bplot
In [ ]:
# pick an object
i = 2
fig = plt.subplots(1, 1, figsize=(6, 6))
fig, ax, parts = bplot.posterior_predictive(
models_mist, # stellar model grid
idxs_mist[i], # model indices
reds_mist[i], # A(V) draws
dreds_mist[i], # R(V) draws
dists_mist[i], # distance draws
data=phot[i],
data_err=err[i], # data
data_mask=mask[i], # band mask
offset=off_mist, # photometric offsets
psig=2.0, # plot 2-sigma errors
labels=filt, # filters
vcolor="blue", # "violin plot" colors for the posteriors
pcolor="black",
fig=fig,
) # photometry colors for the data
In [20]:
# Convert Gaia parallax to distance in kpc
valid_gaia = (parallax > 0) & (~np.isnan(parallax))
gaia_dist = np.full_like(parallax, np.nan)
gaia_dist_err = np.full_like(parallax, np.nan)
gaia_dist[valid_gaia] = 1 / parallax[valid_gaia]
gaia_dist_err[valid_gaia] = (1 / parallax[valid_gaia] ** 2) * parallax_err[valid_gaia]
# Compute percentiles for Gaia distances in the fitted distributions
median_dist = np.median(dists_mist, axis=1)
percentiles_lower = np.full_like(parallax, np.nan)
percentiles_upper = np.full_like(parallax, np.nan)
for i in np.where(valid)[0]:
percentiles_lower[i] = np.median(dists_mist[i, :]) - np.percentile(dists_mist[i, :], 16)
percentiles_upper[i] = np.percentile(dists_mist[i, :], 84) - np.median(dists_mist[i, :])
percentiles = np.vstack([percentiles_lower, percentiles_upper])
i_mag = -2.5 * np.log10(phot[:, 2])
mask = valid_gaia & (~np.isnan(median_dist))
# Plot
plt.figure(figsize=(10, 8))
plt.errorbar(median_dist[mask], gaia_dist[mask], yerr=gaia_dist_err[mask], xerr=percentiles[:, mask], fmt="o", capsize=5)
# plt.colorbar()
plt.xlabel("Fitted Distance (kpc)")
plt.ylabel("Gaia Distance (kpc)")
plt.title("Gaia Distance vs Fitted Distance with Percentiles")
plt.plot([0, 5], [0, 5], "r--", alpha=0.5) # 1:1 line
plt.xlim(0, 5)
plt.ylim(0, 5)
plt.show()
In [28]:
plt.figure(figsize=(10, 6))
# Determine colors based on data availability
has_2mass = np.any(mask_band[:, 5:8], axis=1)
colors = []
for i in range(len(idxs_mist)):
if valid_gaia[i] and has_2mass[i]:
colors.append("purple")
elif valid_gaia[i]:
colors.append("blue")
elif has_2mass[i]:
colors.append("red")
else:
colors.append("black")
for i in range(len(idxs_mist)):
# plt.plot(dists_mist[i, :], reds_mist[i, :], "k.", alpha=0.1)
# plot with average and error bar
plt.errorbar(
np.median(dists_mist[i, :]),
np.median(reds_mist[i, :]),
xerr=[
[np.median(dists_mist[i, :]) - np.percentile(dists_mist[i, :], 16)],
[np.percentile(dists_mist[i, :], 84) - np.median(dists_mist[i, :])],
],
yerr=[
[np.median(reds_mist[i, :]) - np.percentile(reds_mist[i, :], 16)],
[np.percentile(reds_mist[i, :], 84) - np.median(reds_mist[i, :])],
],
fmt="o",
color=colors[i],
alpha=0.7,
)
plt.xlabel("Distance (kpc)")
plt.ylabel("A(V) (mag)")
plt.xlim(0, 5)
# Add legend
import matplotlib.patches as mpatches
black_patch = mpatches.Patch(color="black", label="DECaPS only")
red_patch = mpatches.Patch(color="red", label="D + 2MASS")
blue_patch = mpatches.Patch(color="blue", label="D + Gaia")
purple_patch = mpatches.Patch(color="purple", label="D + Gaia + 2MASS")
plt.legend(handles=[black_patch, red_patch, blue_patch, purple_patch])
plt.show()
