Tutorial 5: Fitting Individual Stars with BruteForce

Tutorial 5: Fitting Individual Stars with BruteForce#

This tutorial demonstrates how to fit individual stars using brutus’s BruteForce class, which performs fast Bayesian inference over pre-computed model grids.

Topics Covered#

  1. Data preparation (flux units, parallaxes, coordinates)

  2. BruteForce class setup and configuration

  3. Running fits with various options

  4. Result visualization and interpretation

Prerequisites#

This tutorial requires the following brutus data files:

  • grid_mist_v9.h5 - MIST model grid

  • offsets_mist_v9.txt - Photometric offsets

  • Orion_l209.1_b-19.9.h5 - Example Orion field data (included in tutorials/)

  • bayestar2019_v1.h5 (optional) - 3D dust map

If you don’t have these files, run the download cell below.

# Download required data files (if not already available locally)
from tutorial_utils import find_brutus_data_file
from brutus.data import fetch_grids, fetch_offsets

# Download grid and offsets if needed
fetch_grids()       # Downloads MIST model grid
fetch_offsets('mist_v9')  # Downloads photometric offsets

# Check dustmap availability (large file, ~2 GB)
try:
    dustmap_path = find_brutus_data_file('bayestar2019_v1.h5')
    print(f"Dustmap found: {dustmap_path}")
except FileNotFoundError:
    from brutus.data import fetch_dustmaps
    fetch_dustmaps()  # Downloads Bayestar dust map
    print("Dustmap downloaded.")

# Note: The Orion field data is included in the tutorials/ directory.

Dustmap found: /root/.cache/astro-brutus/bayestar2019_v1.h5

# Imports and setup
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
import warnings
warnings.filterwarnings('ignore')

from tutorial_utils import (
    setup_tutorial,
    find_brutus_data_file,
    save_figure as _save_fig,
    print_section,
    load_orion_data,
)

info = setup_tutorial(5, title="Tutorial 05: Fitting Individual Stars with BruteForce")
plots_dir = info['plot_dir']


def save_figure(fig, name):
    """Save figure to this tutorial's plot directory."""
    _save_fig(fig, 5, name)

Tutorial 05: Fitting Individual Stars with BruteForce
=====================================================

Checking data requirements for Tutorial 5
=========================================
  Found: grid_mist_v9.h5
  Found: offsets_mist_v9.txt
  Found: bayestar2019_v1.h5
  Found: Orion_l209.1_b-19.9.h5

  All required files available

1. Data Preparation for BruteForce#

BruteForce requires specific data formats:

  • Photometry in flux units (maggies = 10^(-0.4 * magnitude))

  • Parallaxes in mas (milliarcseconds) from Gaia

  • Galactic coordinates (l, b) for applying priors

Let’s load and examine the Orion field data.

from brutus.utils import inv_magnitude, magnitude

# Load example Orion field data
print("Loading Orion field data...")
data = load_orion_data()

print(f"\nLoaded {len(data['phot'])} sources")
print(f"  Photometry shape: {data['phot'].shape}")
print(f"  Filters: Pan-STARRS (grizy) + 2MASS (JHKs) = 8 bands")
print(f"  Valid parallaxes: {np.sum(np.isfinite(data['parallax']))}")
print(f"  Coordinate range: l=[{data['coords'][:, 0].min():.1f}, {data['coords'][:, 0].max():.1f}], "
      f"b=[{data['coords'][:, 1].min():.1f}, {data['coords'][:, 1].max():.1f}] deg")

Loading Orion field data...

Loaded 207 sources
  Photometry shape: (207, 8)
  Filters: Pan-STARRS (grizy) + 2MASS (JHKs) = 8 bands
  Valid parallaxes: 199
  Coordinate range: l=[208.9, 209.3], b=[-20.1, -19.7] deg

# Add systematic photometric uncertainties in quadrature
# The MIST isochrone corrections have systematic uncertainties of ~0.02 mag
# in the optical and ~0.03 mag in the NIR (Speagle et al. 2025, Table 5).
# These should be added to the measurement errors before fitting.

# Systematic uncertainties per band (mag)
sys_err_mag = np.array([0.02, 0.02, 0.02, 0.02, 0.02,   # PS grizy
                        0.03, 0.03, 0.03])                # 2MASS JHKs

# Convert to flux space and add in quadrature
sys_err_flux = data['phot'] * sys_err_mag[np.newaxis, :] * np.log(10) / 2.5
data['err'] = np.sqrt(data['err']**2 + sys_err_flux**2)

# Show the effect
mag_err_after = 2.5 / np.log(10) * data['err'] / np.where(data['phot'] > 0, data['phot'], 1)
filt_names = ['PS_g', 'PS_r', 'PS_i', 'PS_z', 'PS_y', '2MASS_J', '2MASS_H', '2MASS_Ks']
print("Effective magnitude errors after adding systematics:")
for i, name in enumerate(filt_names):
    valid = data['mask'][:, i]
    if np.sum(valid) > 0:
        print(f"  {name:>10s}: min={np.min(mag_err_after[valid, i]):.4f}, "
              f"median={np.median(mag_err_after[valid, i]):.4f} mag")

Effective magnitude errors after adding systematics:
        PS_g: min=0.0288, median=0.0425 mag
        PS_r: min=0.0286, median=0.0306 mag
        PS_i: min=0.0285, median=0.0288 mag
        PS_z: min=0.0285, median=0.0289 mag
        PS_y: min=0.0286, median=0.0295 mag
     2MASS_J: min=0.0366, median=0.0439 mag
     2MASS_H: min=0.0372, median=0.0461 mag
    2MASS_Ks: min=0.0350, median=0.0431 mag

2. BruteForce Class Setup#

BruteForce performs fast stellar parameter inference by:

  1. Evaluating likelihood over a pre-computed grid of stellar models

  2. Marginalizing analytically over nuisance parameters (distance, extinction)

  3. Incorporating priors (Galactic structure, 3D dust, parallax)

Let’s load the MIST model grid and initialize BruteForce.

from brutus.analysis import BruteForce
from brutus.core import StarGrid
from brutus.data import load_models, filters

# Load MIST model grid
print("Loading MIST model grid...")
grid_file = find_brutus_data_file('grid_mist_v9.h5')

# Define filters (Pan-STARRS + 2MASS)
filt = filters.ps[:-2] + filters.tmass  # Skip PS w and open bands
print(f"\nUsing filters: {filt}")

# Load models
models, labels, label_mask = load_models(grid_file, filters=filt)

print(f"\nLoaded {len(models):,} models")
print(f"  Model shape: {models.shape}")
print(f"  Parameter fields: {list(labels.dtype.names)}")
print(f"  Key parameters:")
print(f"    Mass range: [{labels['mini'].min():.2f}, {labels['mini'].max():.1f}] M_sun")
print(f"    [Fe/H] range: [{labels['feh'].min():.2f}, {labels['feh'].max():.2f}]")
print(f"    log(age) range: [{labels['loga'].min():.2f}, {labels['loga'].max():.2f}]")

# Create StarGrid and initialize BruteForce
grid = StarGrid(models, labels)
bf = BruteForce(grid)

print(f"\nBruteForce initialized with {len(filt)} bands")

Loading MIST model grid...

Using filters: ['PS_g', 'PS_r', 'PS_i', 'PS_z', 'PS_y', '2MASS_J', '2MASS_H', '2MASS_Ks']

Reading entire dataset (49 filters) once...

Extracting 8 requested filters from memory: ['PS_g', 'PS_r', 'PS_i', 'PS_z', 'PS_y', '2MASS_J', '2MASS_H', '2MASS_Ks']

Dropping constant 'afe' label column (single-afe grid, afe=0)...

Loaded 613,530 models
  Model shape: (613530, 8, 3)
  Parameter fields: ['mini', 'feh', 'eep', 'loga', 'logl', 'logt', 'logg', 'agewt']
  Key parameters:
    Mass range: [0.50, 2.0] M_sun
    [Fe/H] range: [-3.00, 0.45]
    log(age) range: [6.46, 10.14]
Loaded StarGrid with 613,530 models, 8 filters, 8 labels
BruteForce initialized with 613,530 models
  Grid parameters (8): mini, feh, eep, loga, logl, logt, logg, agewt

BruteForce initialized with 8 bands

3. Running BruteForce Fits#

Now let’s fit all sources from the Orion field, demonstrating various fitting options.

Key considerations:

  • Photometric offsets calibrate models to match observed photometry

  • 3D dust maps provide extinction priors

  • Parallax constraints from Gaia improve distance estimates

  • Galactic priors based on position help constrain parameters

# Load photometric offsets (calibration corrections)
from brutus.data import load_offsets

print("Loading photometric offsets...")
try:
    offsets_file = find_brutus_data_file('offsets_mist_v9.txt')
    offsets_mist = load_offsets(offsets_file, filters=filt)
    print("\nPhotometric offsets (multiplicative flux):")
    for band, offset in zip(filt, offsets_mist):
        print(f"  {band}: {offset:.4f}")
except FileNotFoundError:
    offsets_mist = np.ones(len(filt))
    print("  Offsets file not found, using unity offsets")

# Select sources with sufficient photometric coverage
good_sources = np.sum(data['mask'], axis=1) >= 4  # At least 4 bands
idx_fit = np.where(good_sources)[0]

print(f"\nSelected {len(idx_fit)} sources for fitting")
print(f"  Criteria: >=4 bands detected")
print(f"  Average bands per source: {np.mean(np.sum(data['mask'][idx_fit], axis=1)):.1f}")

# Try to load 3D dust map
try:
    dustfile = find_brutus_data_file('bayestar2019_v1.h5')
    print("\n3D dust map (Bayestar) available")
except FileNotFoundError:
    dustfile = None
    print("\n3D dust map not available (will use flat prior)")

Loading photometric offsets...

Photometric offsets (multiplicative flux):
  PS_g: 1.0000
  PS_r: 0.9500
  PS_i: 0.9600
  PS_z: 0.9400
  PS_y: 0.9500
  2MASS_J: 0.9800
  2MASS_H: 1.0400
  2MASS_Ks: 1.0300

Selected 207 sources for fitting
  Criteria: >=4 bands detected
  Average bands per source: 7.4

3D dust map (Bayestar) available

PS_g (0.0%)
PS_r (-5.0%)
PS_i (-4.0%)
PS_z (-6.0%)
PS_y (-5.0%)
2MASS_J (-2.0%)
2MASS_H (4.0%)
2MASS_Ks (3.0%)

Note: brutus recommends at least 4 photometric bands for reliable fits. Users can disable the 3D dust map prior (making it flat over A(V)) by not passing the dustfile argument to bf.fit().

# Run BruteForce fits
# Save results to tutorials/ directory for use by later tutorials
tutorials_dir = Path(__file__).parent if '__file__' in dir() else Path('.')
output_file = tutorials_dir / 'Orion_l209.1_b-19.9_mist.h5'

import time

print(f"Running BruteForce fits on {len(idx_fit)} sources...")
print("This may take several minutes depending on your hardware.")

t0 = time.time()

bf.fit(
    data=data['phot'][idx_fit],
    data_err=data['err'][idx_fit],
    data_mask=data['mask'][idx_fit],
    data_labels=idx_fit,
    save_file=str(output_file),
    data_coords=data['coords'][idx_fit],
    parallax=data['parallax'][idx_fit],
    parallax_err=data['parallax_err'][idx_fit],
    phot_offsets=offsets_mist,
    dustfile=dustfile,
    Ndraws=500,
    Nmc_prior=20,
    logl_dim_prior=True,
    save_dar_draws=True,
    running_io=False,
    verbose=True
)

elapsed = time.time() - t0
print(f"\nFitting complete in {elapsed:.1f} seconds ({elapsed/60:.1f} minutes).")
print(f"Results saved to {output_file}")

Running BruteForce fits on 207 sources...
This may take several minutes depending on your hardware.

Fitting object 1/207  

Fitting object 2/207 [chi2/n: 3.3/9] (mean time: 30.330 s/obj, est. remaining: 6247.931 s)    

Fitting object 3/207 [chi2/n: 50.5/9] (mean time: 15.414 s/obj, est. remaining: 3159.811 s)    

Fitting object 4/207 [chi2/n: 44.0/9] (mean time: 10.485 s/obj, est. remaining: 2139.017 s)    

Fitting object 5/207 [chi2/n: 1.8/9] (mean time: 8.014 s/obj, est. remaining: 1626.827 s)    

Fitting object 6/207 [chi2/n: 44.7/9] (mean time: 6.524 s/obj, est. remaining: 1317.779 s)    

Fitting object 7/207 [chi2/n: 60.2/9] (mean time: 5.517 s/obj, est. remaining: 1108.904 s)    

Fitting object 8/207 [chi2/n: 51.4/9] (mean time: 4.825 s/obj, est. remaining: 964.938 s)    

Fitting object 9/207 [chi2/n: 30.8/9] (mean time: 4.293 s/obj, est. remaining: 854.349 s)    

Fitting object 10/207 [chi2/n: 25.2/9] (mean time: 3.879 s/obj, est. remaining: 768.047 s)    

Fitting object 11/207 [chi2/n: 12.6/9] (mean time: 3.608 s/obj, est. remaining: 710.780 s)    

Fitting object 12/207 [chi2/n: 12.0/8] (mean time: 3.354 s/obj, est. remaining: 657.312 s)    

Fitting object 13/207 [chi2/n: 6.0/6] (mean time: 3.110 s/obj, est. remaining: 606.392 s)    

Fitting object 14/207 [chi2/n: 48.3/9] (mean time: 2.908 s/obj, est. remaining: 564.237 s)    

Fitting object 15/207 [chi2/n: 99.2/9] (mean time: 2.735 s/obj, est. remaining: 527.838 s)    

Fitting object 16/207 [chi2/n: 140.2/9] (mean time: 2.585 s/obj, est. remaining: 496.334 s)    

Fitting object 17/207 [chi2/n: 86.1/9] (mean time: 2.460 s/obj, est. remaining: 469.919 s)    

Fitting object 18/207 [chi2/n: 39.3/9] (mean time: 2.348 s/obj, est. remaining: 446.151 s)    

Fitting object 19/207 [chi2/n: 62.8/8] (mean time: 2.297 s/obj, est. remaining: 434.217 s)    

Fitting object 20/207 [chi2/n: 55.6/8] (mean time: 2.208 s/obj, est. remaining: 415.081 s)    

Fitting object 21/207 [chi2/n: 47.3/5] (mean time: 2.130 s/obj, est. remaining: 398.359 s)    

Fitting object 22/207 [chi2/n: 43.3/9] (mean time: 2.064 s/obj, est. remaining: 383.892 s)    

Fitting object 23/207 [chi2/n: 34.6/8] (mean time: 1.994 s/obj, est. remaining: 368.842 s)    

Fitting object 24/207 [chi2/n: 267.1/9] (mean time: 1.928 s/obj, est. remaining: 354.699 s)    

Fitting object 25/207 [chi2/n: 4.6/9] (mean time: 1.866 s/obj, est. remaining: 341.496 s)    

Fitting object 26/207 [chi2/n: 210.0/9] (mean time: 1.811 s/obj, est. remaining: 329.573 s)    

Fitting object 27/207 [chi2/n: 192.6/9] (mean time: 1.759 s/obj, est. remaining: 318.321 s)    

Fitting object 28/207 [chi2/n: 113.3/9] (mean time: 1.712 s/obj, est. remaining: 308.222 s)    

Fitting object 29/207 [chi2/n: 74.9/8] (mean time: 1.669 s/obj, est. remaining: 298.830 s)    

Fitting object 30/207 [chi2/n: 24.8/9] (mean time: 1.629 s/obj, est. remaining: 289.912 s)    

Fitting object 31/207 [chi2/n: 15.9/9] (mean time: 1.614 s/obj, est. remaining: 285.764 s)    

Fitting object 32/207 [chi2/n: 18.9/9] (mean time: 1.580 s/obj, est. remaining: 278.056 s)    

Fitting object 33/207 [chi2/n: 93.8/9] (mean time: 1.545 s/obj, est. remaining: 270.375 s)    

Fitting object 34/207 [chi2/n: 41.6/9] (mean time: 1.515 s/obj, est. remaining: 263.550 s)    

Fitting object 35/207 [chi2/n: 72.4/9] (mean time: 1.487 s/obj, est. remaining: 257.263 s)    

Fitting object 36/207 [chi2/n: 13.1/9] (mean time: 1.459 s/obj, est. remaining: 251.007 s)    

Fitting object 37/207 [chi2/n: 42.7/9] (mean time: 1.432 s/obj, est. remaining: 244.895 s)    

Fitting object 38/207 [chi2/n: 16.5/9] (mean time: 1.417 s/obj, est. remaining: 240.953 s)    

Fitting object 39/207 [chi2/n: 11.2/8] (mean time: 1.400 s/obj, est. remaining: 236.520 s)    

Fitting object 40/207 [chi2/n: 76.0/9] (mean time: 1.375 s/obj, est. remaining: 231.080 s)    

Fitting object 41/207 [chi2/n: 17.8/9] (mean time: 1.363 s/obj, est. remaining: 227.697 s)    

Fitting object 42/207 [chi2/n: 5.2/8] (mean time: 1.351 s/obj, est. remaining: 224.288 s)    

Fitting object 43/207 [chi2/n: 8.6/8] (mean time: 1.341 s/obj, est. remaining: 221.295 s)    

Fitting object 44/207 [chi2/n: 9.0/9] (mean time: 1.332 s/obj, est. remaining: 218.382 s)    

Fitting object 45/207 [chi2/n: 2.9/8] (mean time: 1.318 s/obj, est. remaining: 214.809 s)    

Fitting object 46/207 [chi2/n: 54.3/9] (mean time: 1.299 s/obj, est. remaining: 210.461 s)    

Fitting object 47/207 [chi2/n: 11.4/8] (mean time: 1.290 s/obj, est. remaining: 207.755 s)    

Fitting object 48/207 [chi2/n: 8.8/8] (mean time: 1.283 s/obj, est. remaining: 205.204 s)    

Fitting object 49/207 [chi2/n: 6.3/9] (mean time: 1.282 s/obj, est. remaining: 203.775 s)    

Fitting object 50/207 [chi2/n: 2.2/7] (mean time: 1.273 s/obj, est. remaining: 201.089 s)    

Fitting object 51/207 [chi2/n: 21.0/9] (mean time: 1.272 s/obj, est. remaining: 199.714 s)    

Fitting object 52/207 [chi2/n: 9.9/8] (mean time: 1.265 s/obj, est. remaining: 197.389 s)    

Fitting object 53/207 [chi2/n: 5.2/8] (mean time: 1.260 s/obj, est. remaining: 195.268 s)    

Fitting object 54/207 [chi2/n: 8.0/9] (mean time: 1.255 s/obj, est. remaining: 193.237 s)    

Fitting object 55/207 [chi2/n: 5.2/6] (mean time: 1.243 s/obj, est. remaining: 190.183 s)    

Fitting object 56/207 [chi2/n: 2.5/6] (mean time: 1.243 s/obj, est. remaining: 188.889 s)    

Fitting object 57/207 [chi2/n: 1.1/6] (mean time: 1.234 s/obj, est. remaining: 186.327 s)    

Fitting object 58/207 [chi2/n: 0.5/6] (mean time: 1.225 s/obj, est. remaining: 183.755 s)    

Fitting object 59/207 [chi2/n: 0.2/6] (mean time: 1.217 s/obj, est. remaining: 181.279 s)    

Fitting object 60/207 [chi2/n: 25.0/9] (mean time: 1.206 s/obj, est. remaining: 178.459 s)    

Fitting object 61/207 [chi2/n: 0.3/6] (mean time: 1.198 s/obj, est. remaining: 176.076 s)    

Fitting object 62/207 [chi2/n: 27.5/9] (mean time: 1.187 s/obj, est. remaining: 173.268 s)    

Fitting object 63/207 [chi2/n: 21.8/9] (mean time: 1.176 s/obj, est. remaining: 170.579 s)    

Fitting object 64/207 [chi2/n: 8.3/9] (mean time: 1.168 s/obj, est. remaining: 168.171 s)    

Fitting object 65/207 [chi2/n: 67.4/9] (mean time: 1.158 s/obj, est. remaining: 165.613 s)    

Fitting object 66/207 [chi2/n: 23.0/9] (mean time: 1.152 s/obj, est. remaining: 163.562 s)    

Fitting object 67/207 [chi2/n: 14.4/9] (mean time: 1.152 s/obj, est. remaining: 162.395 s)    

Fitting object 68/207 [chi2/n: 8.1/6] (mean time: 1.144 s/obj, est. remaining: 160.091 s)    

Fitting object 69/207 [chi2/n: 6.8/9] (mean time: 1.151 s/obj, est. remaining: 160.030 s)    

Fitting object 70/207 [chi2/n: 407.0/9] (mean time: 1.142 s/obj, est. remaining: 157.656 s)    

Fitting object 71/207 [chi2/n: 3146.0/9] (mean time: 1.134 s/obj, est. remaining: 155.362 s)    

Fitting object 72/207 [chi2/n: 1537.7/9] (mean time: 1.125 s/obj, est. remaining: 153.062 s)    

Fitting object 73/207 [chi2/n: 1473.8/9] (mean time: 1.117 s/obj, est. remaining: 150.808 s)    

Fitting object 74/207 [chi2/n: 884.2/9] (mean time: 1.109 s/obj, est. remaining: 148.575 s)    

Fitting object 75/207 [chi2/n: 708.9/9] (mean time: 1.101 s/obj, est. remaining: 146.393 s)    

Fitting object 76/207 [chi2/n: 66.8/9] (mean time: 1.093 s/obj, est. remaining: 144.304 s)    

Fitting object 77/207 [chi2/n: 1808.9/9] (mean time: 1.086 s/obj, est. remaining: 142.278 s)    

Fitting object 78/207 [chi2/n: 176.8/9] (mean time: 1.079 s/obj, est. remaining: 140.296 s)    

Fitting object 79/207 [chi2/n: 9.5/9] (mean time: 1.073 s/obj, est. remaining: 138.445 s)    

Fitting object 80/207 [chi2/n: 26.2/9] (mean time: 1.066 s/obj, est. remaining: 136.509 s)    

Fitting object 81/207 [chi2/n: 72.2/9] (mean time: 1.060 s/obj, est. remaining: 134.580 s)    

Fitting object 82/207 [chi2/n: 797.4/9] (mean time: 1.053 s/obj, est. remaining: 132.708 s)    

Fitting object 83/207 [chi2/n: 0.2/6] (mean time: 1.045 s/obj, est. remaining: 130.653 s)    

Fitting object 84/207 [chi2/n: 1553.4/8] (mean time: 1.039 s/obj, est. remaining: 128.837 s)    

Fitting object 85/207 [chi2/n: 20.6/9] (mean time: 1.033 s/obj, est. remaining: 127.082 s)    

Fitting object 86/207 [chi2/n: 40.5/9] (mean time: 1.028 s/obj, est. remaining: 125.460 s)    

Fitting object 87/207 [chi2/n: 148.4/8] (mean time: 1.026 s/obj, est. remaining: 124.136 s)    

Fitting object 88/207 [chi2/n: 1751.2/9] (mean time: 1.021 s/obj, est. remaining: 122.468 s)    

Fitting object 89/207 [chi2/n: 16.5/9] (mean time: 1.022 s/obj, est. remaining: 121.640 s)    

Fitting object 90/207 [chi2/n: 25.8/8] (mean time: 1.015 s/obj, est. remaining: 119.818 s)    

Fitting object 91/207 [chi2/n: 15.0/7] (mean time: 1.014 s/obj, est. remaining: 118.666 s)    

Fitting object 92/207 [chi2/n: 37.5/4] (mean time: 1.007 s/obj, est. remaining: 116.754 s)    

Fitting object 93/207 [chi2/n: 968.9/9] (mean time: 1.002 s/obj, est. remaining: 115.174 s)    

Fitting object 94/207 [chi2/n: 657.6/9] (mean time: 0.997 s/obj, est. remaining: 113.623 s)    

Fitting object 95/207 [chi2/n: 106.8/8] (mean time: 0.992 s/obj, est. remaining: 112.040 s)    

Fitting object 96/207 [chi2/n: 30.5/9] (mean time: 0.987 s/obj, est. remaining: 110.500 s)    

Fitting object 97/207 [chi2/n: 252.8/8] (mean time: 0.982 s/obj, est. remaining: 109.009 s)    

Fitting object 98/207 [chi2/n: 161.7/9] (mean time: 0.978 s/obj, est. remaining: 107.615 s)    

Fitting object 99/207 [chi2/n: 29.2/9] (mean time: 0.975 s/obj, est. remaining: 106.230 s)    

Fitting object 100/207 [chi2/n: 405.7/9] (mean time: 0.971 s/obj, est. remaining: 104.818 s)    

Fitting object 101/207 [chi2/n: 165.0/6] (mean time: 0.966 s/obj, est. remaining: 103.333 s)    

Fitting object 102/207 [chi2/n: 16.9/9] (mean time: 0.962 s/obj, est. remaining: 101.932 s)    

Fitting object 103/207 [chi2/n: 8.4/9] (mean time: 0.958 s/obj, est. remaining: 100.624 s)    

Fitting object 104/207 [chi2/n: 106.1/9] (mean time: 0.955 s/obj, est. remaining: 99.311 s)    

Fitting object 105/207 [chi2/n: 11.4/8] (mean time: 0.951 s/obj, est. remaining: 97.931 s)    

Fitting object 106/207 [chi2/n: 71.5/9] (mean time: 0.947 s/obj, est. remaining: 96.594 s)    

Fitting object 107/207 [chi2/n: 30.0/9] (mean time: 0.943 s/obj, est. remaining: 95.287 s)    

Fitting object 108/207 [chi2/n: 28.2/9] (mean time: 0.941 s/obj, est. remaining: 94.065 s)    

Fitting object 109/207 [chi2/n: 9.1/9] (mean time: 0.940 s/obj, est. remaining: 93.035 s)    

Fitting object 110/207 [chi2/n: 17.1/9] (mean time: 0.942 s/obj, est. remaining: 92.323 s)    

Fitting object 111/207 [chi2/n: 392.7/9] (mean time: 0.938 s/obj, est. remaining: 91.031 s)    

Fitting object 112/207 [chi2/n: 17.2/9] (mean time: 0.942 s/obj, est. remaining: 90.432 s)    

Fitting object 113/207 [chi2/n: 1118.2/9] (mean time: 0.939 s/obj, est. remaining: 89.196 s)    

Fitting object 114/207 [chi2/n: 17.1/7] (mean time: 0.939 s/obj, est. remaining: 88.236 s)    

Fitting object 115/207 [chi2/n: 34.4/6] (mean time: 0.934 s/obj, est. remaining: 86.849 s)    

Fitting object 116/207 [chi2/n: 3018.7/9] (mean time: 0.930 s/obj, est. remaining: 85.526 s)    

Fitting object 117/207 [chi2/n: 0.4/9] (mean time: 0.926 s/obj, est. remaining: 84.229 s)    

Fitting object 118/207 [chi2/n: 82.2/5] (mean time: 0.921 s/obj, est. remaining: 82.911 s)    

Fitting object 119/207 [chi2/n: 13.4/9] (mean time: 0.917 s/obj, est. remaining: 81.643 s)    

Fitting object 120/207 [chi2/n: 24.7/9] (mean time: 0.919 s/obj, est. remaining: 80.875 s)    

Fitting object 121/207 [chi2/n: 33.0/9] (mean time: 0.916 s/obj, est. remaining: 79.720 s)    

Fitting object 122/207 [chi2/n: 29.9/9] (mean time: 0.917 s/obj, est. remaining: 78.823 s)    

Fitting object 123/207 [chi2/n: 107.3/7] (mean time: 0.914 s/obj, est. remaining: 77.655 s)    

Fitting object 124/207 [chi2/n: 27.9/9] (mean time: 0.911 s/obj, est. remaining: 76.526 s)    

Fitting object 125/207 [chi2/n: 1.5/9] (mean time: 0.909 s/obj, est. remaining: 75.421 s)    

Fitting object 126/207 [chi2/n: 1895.5/9] (mean time: 0.906 s/obj, est. remaining: 74.285 s)    

Fitting object 127/207 [chi2/n: 24.5/9] (mean time: 0.903 s/obj, est. remaining: 73.177 s)    

Fitting object 128/207 [chi2/n: 23.3/9] (mean time: 0.901 s/obj, est. remaining: 72.047 s)    

Fitting object 129/207 [chi2/n: 32.3/9] (mean time: 0.898 s/obj, est. remaining: 70.973 s)    

Fitting object 130/207 [chi2/n: 34.5/9] (mean time: 0.896 s/obj, est. remaining: 69.870 s)    

Fitting object 131/207 [chi2/n: 174.1/9] (mean time: 0.894 s/obj, est. remaining: 68.847 s)    

Fitting object 132/207 [chi2/n: 17.0/9] (mean time: 0.892 s/obj, est. remaining: 67.769 s)    

Fitting object 133/207 [chi2/n: 26.5/9] (mean time: 0.890 s/obj, est. remaining: 66.766 s)    

Fitting object 134/207 [chi2/n: 32.9/9] (mean time: 0.889 s/obj, est. remaining: 65.754 s)    

Fitting object 135/207 [chi2/n: 17.9/9] (mean time: 0.886 s/obj, est. remaining: 64.711 s)    

Fitting object 136/207 [chi2/n: 74.7/9] (mean time: 0.884 s/obj, est. remaining: 63.640 s)    

Fitting object 137/207 [chi2/n: 24.0/9] (mean time: 0.882 s/obj, est. remaining: 62.588 s)    

Fitting object 138/207 [chi2/n: 148.2/9] (mean time: 0.879 s/obj, est. remaining: 61.547 s)    

Fitting object 139/207 [chi2/n: 16.4/9] (mean time: 0.877 s/obj, est. remaining: 60.536 s)    

Fitting object 140/207 [chi2/n: 32.3/9] (mean time: 0.875 s/obj, est. remaining: 59.508 s)    

Fitting object 141/207 [chi2/n: 15.9/9] (mean time: 0.874 s/obj, est. remaining: 58.527 s)    

Fitting object 142/207 [chi2/n: 34.5/9] (mean time: 0.873 s/obj, est. remaining: 57.641 s)    

Fitting object 143/207 [chi2/n: 38.0/8] (mean time: 0.875 s/obj, est. remaining: 56.894 s)    

Fitting object 144/207 [chi2/n: 48.7/9] (mean time: 0.873 s/obj, est. remaining: 55.856 s)    

Fitting object 145/207 [chi2/n: 18.8/9] (mean time: 0.871 s/obj, est. remaining: 54.864 s)    

Fitting object 146/207 [chi2/n: 32.9/9] (mean time: 0.869 s/obj, est. remaining: 53.849 s)    

Fitting object 147/207 [chi2/n: 15.9/9] (mean time: 0.866 s/obj, est. remaining: 52.848 s)    

Fitting object 148/207 [chi2/n: 22.0/9] (mean time: 0.864 s/obj, est. remaining: 51.848 s)    

Fitting object 149/207 [chi2/n: 3.9/9] (mean time: 0.865 s/obj, est. remaining: 51.017 s)    

Fitting object 150/207 [chi2/n: 109.9/9] (mean time: 0.862 s/obj, est. remaining: 50.021 s)    

Fitting object 151/207 [chi2/n: 21.6/9] (mean time: 0.861 s/obj, est. remaining: 49.056 s)    

Fitting object 152/207 [chi2/n: 17.1/9] (mean time: 0.858 s/obj, est. remaining: 48.060 s)    

Fitting object 153/207 [chi2/n: 21.8/9] (mean time: 0.856 s/obj, est. remaining: 47.088 s)    

Fitting object 154/207 [chi2/n: 106.0/9] (mean time: 0.854 s/obj, est. remaining: 46.109 s)    

Fitting object 155/207 [chi2/n: 8.9/9] (mean time: 0.856 s/obj, est. remaining: 45.352 s)    

Fitting object 156/207 [chi2/n: 23.1/9] (mean time: 0.854 s/obj, est. remaining: 44.391 s)    

Fitting object 157/207 [chi2/n: 11.1/8] (mean time: 0.855 s/obj, est. remaining: 43.596 s)    

Fitting object 158/207 [chi2/n: 2.7/8] (mean time: 0.856 s/obj, est. remaining: 42.796 s)    

Fitting object 159/207 [chi2/n: 11.6/8] (mean time: 0.856 s/obj, est. remaining: 41.965 s)    

Fitting object 160/207 [chi2/n: 21.2/9] (mean time: 0.854 s/obj, est. remaining: 41.014 s)    

Fitting object 161/207 [chi2/n: 15.9/9] (mean time: 0.858 s/obj, est. remaining: 40.333 s)    

Fitting object 162/207 [chi2/n: 16.5/9] (mean time: 0.858 s/obj, est. remaining: 39.475 s)    

Fitting object 163/207 [chi2/n: 11.8/8] (mean time: 0.856 s/obj, est. remaining: 38.502 s)    

Fitting object 164/207 [chi2/n: 0.9/7] (mean time: 0.856 s/obj, est. remaining: 37.664 s)    

Fitting object 165/207 [chi2/n: 9.1/8] (mean time: 0.857 s/obj, est. remaining: 36.849 s)    

Fitting object 166/207 [chi2/n: 3.0/6] (mean time: 0.857 s/obj, est. remaining: 36.006 s)    

Fitting object 167/207 [chi2/n: 28.8/8] (mean time: 0.859 s/obj, est. remaining: 35.226 s)    

Fitting object 168/207 [chi2/n: 0.7/6] (mean time: 0.858 s/obj, est. remaining: 34.334 s)    

Fitting object 169/207 [chi2/n: 1.6/6] (mean time: 0.859 s/obj, est. remaining: 33.504 s)    

Fitting object 170/207 [chi2/n: 0.2/6] (mean time: 0.859 s/obj, est. remaining: 32.632 s)    

Fitting object 171/207 [chi2/n: 1.3/6] (mean time: 0.861 s/obj, est. remaining: 31.874 s)    

Fitting object 172/207 [chi2/n: 1.5/6] (mean time: 0.861 s/obj, est. remaining: 30.994 s)    

Fitting object 173/207 [chi2/n: 0.7/6] (mean time: 0.864 s/obj, est. remaining: 30.225 s)    

Fitting object 174/207 [chi2/n: 1.8/6] (mean time: 0.863 s/obj, est. remaining: 29.354 s)    

Fitting object 175/207 [chi2/n: 1.1/6] (mean time: 0.863 s/obj, est. remaining: 28.469 s)    

Fitting object 176/207 [chi2/n: 1.9/6] (mean time: 0.862 s/obj, est. remaining: 27.590 s)    

Fitting object 177/207 [chi2/n: 0.8/6] (mean time: 0.862 s/obj, est. remaining: 26.737 s)    

Fitting object 178/207 [chi2/n: 0.7/6] (mean time: 0.862 s/obj, est. remaining: 25.863 s)    

Fitting object 179/207 [chi2/n: 10.4/9] (mean time: 0.860 s/obj, est. remaining: 24.943 s)    

Fitting object 180/207 [chi2/n: 10.4/9] (mean time: 0.858 s/obj, est. remaining: 24.025 s)    

Fitting object 181/207 [chi2/n: 9.4/9] (mean time: 0.856 s/obj, est. remaining: 23.116 s)    

Fitting object 182/207 [chi2/n: 16.4/9] (mean time: 0.855 s/obj, est. remaining: 22.219 s)    

Fitting object 183/207 [chi2/n: 7.2/9] (mean time: 0.853 s/obj, est. remaining: 21.315 s)    

Fitting object 184/207 [chi2/n: 59.0/9] (mean time: 0.851 s/obj, est. remaining: 20.419 s)    

Fitting object 185/207 [chi2/n: 9.6/9] (mean time: 0.849 s/obj, est. remaining: 19.530 s)    

Fitting object 186/207 [chi2/n: 22.0/9] (mean time: 0.848 s/obj, est. remaining: 18.646 s)    

Fitting object 187/207 [chi2/n: 25.4/9] (mean time: 0.851 s/obj, est. remaining: 17.878 s)    

Fitting object 188/207 [chi2/n: 16.5/6] (mean time: 0.849 s/obj, est. remaining: 16.986 s)    

Fitting object 189/207 [chi2/n: 39.7/9] (mean time: 0.847 s/obj, est. remaining: 16.096 s)    

Fitting object 190/207 [chi2/n: 9.3/9] (mean time: 0.847 s/obj, est. remaining: 15.243 s)    

Fitting object 191/207 [chi2/n: 3.9/9] (mean time: 0.848 s/obj, est. remaining: 14.420 s)    

Fitting object 192/207 [chi2/n: 11.7/9] (mean time: 0.846 s/obj, est. remaining: 13.544 s)    

Fitting object 193/207 [chi2/n: 13.3/9] (mean time: 0.846 s/obj, est. remaining: 12.696 s)    

Fitting object 194/207 [chi2/n: 16.8/9] (mean time: 0.845 s/obj, est. remaining: 11.826 s)    

Fitting object 195/207 [chi2/n: 8.4/9] (mean time: 0.851 s/obj, est. remaining: 11.059 s)    

Fitting object 196/207 [chi2/n: 19.5/6] (mean time: 0.849 s/obj, est. remaining: 10.191 s)    

Fitting object 197/207 [chi2/n: 2.5/6] (mean time: 0.848 s/obj, est. remaining: 9.332 s)    

Fitting object 198/207 [chi2/n: 10.6/9] (mean time: 0.846 s/obj, est. remaining: 8.464 s)    

Fitting object 199/207 [chi2/n: 30.8/8] (mean time: 0.844 s/obj, est. remaining: 7.600 s)    

Fitting object 200/207 [chi2/n: 13.4/9] (mean time: 0.845 s/obj, est. remaining: 6.762 s)    

Fitting object 201/207 [chi2/n: 4.3/6] (mean time: 0.843 s/obj, est. remaining: 5.902 s)    

Fitting object 202/207 [chi2/n: 76.3/9] (mean time: 0.841 s/obj, est. remaining: 5.048 s)    

Fitting object 203/207 [chi2/n: 20.9/9] (mean time: 0.840 s/obj, est. remaining: 4.199 s)    

Fitting object 204/207 [chi2/n: 13.2/9] (mean time: 0.842 s/obj, est. remaining: 3.366 s)    

Fitting object 205/207 [chi2/n: 114.8/8] (mean time: 0.840 s/obj, est. remaining: 2.519 s)    

Fitting object 206/207 [chi2/n: 3.3/6] (mean time: 0.839 s/obj, est. remaining: 1.678 s)    

Fitting object 207/207 [chi2/n: 15.6/6] (mean time: 0.837 s/obj, est. remaining: 0.837 s)    

Fitting object 207/207 [chi2/n: 19.3/9] (mean time: 0.841 s/obj, total: 174.002 s)    

Fitting complete in 185.4 seconds (3.1 minutes).
Results saved to Orion_l209.1_b-19.9_mist.h5

Hide code cell source

 
# Load and analyze results
import h5py

with h5py.File(output_file, 'r') as f:
    chi2 = f['obj_chi2min'][:]
    nbands = f['obj_Nbands'][:]
    model_idx = f['model_idx'][:]
    ml_scale = f['ml_scale'][:]
    ml_av = f['ml_av'][:]
    ml_rv = f['ml_rv'][:]
    log_post = f['obj_log_post'][:]
    dists = f['samps_dist'][:]
    reds = f['samps_red'][:]
    dreds = f['samps_dred'][:]

print("Loaded fitting results:")
print(f"  chi2 values: {chi2.shape}")
print(f"  Model indices: {model_idx.shape}")
print(f"  ML scale/AV/RV: {ml_scale.shape}")
print(f"  Distance samples: {dists.shape}")

# Compute summary statistics
chi2_per_band = chi2 / nbands  # shown below as a diagnostic only
# Fit quality is judged by the goodness-of-fit p-value (dof = Nbands - 3
# free params: scale, A(V), R(V)), NOT chi2/Nbands -- the latter is not
# variance-stabilizing across different band counts (see CLAUDE.md).
from scipy.stats import chi2 as chi2_dist
dof = np.maximum(nbands - 3, 1)
pvals = chi2_dist.sf(chi2, dof)
mean_dists = np.mean(dists, axis=1)
mean_av = np.mean(reds, axis=1)

# Create fit quality visualization
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

# Panel 1: Chi2 distribution
ax = axes[0, 0]
ax.hist(chi2_per_band, bins=30, alpha=0.7, color='blue', edgecolor='darkblue')
ax.axvline(1.0, color='red', ls='--', lw=2, label='Expected (chi2/n = 1)')
ax.axvline(np.median(chi2_per_band), color='green', ls='--', lw=2, label=f'Median = {np.median(chi2_per_band):.2f}')
ax.set_xlabel('chi2/band')
ax.set_ylabel('Number of Sources')
ax.set_title('Goodness of Fit')
ax.set_xlim(0, np.percentile(chi2_per_band, 99))
ax.legend()
ax.grid(True, alpha=0.3)

# Panel 2: Distance distribution
ax = axes[0, 1]
ax.hist(mean_dists, bins=30, alpha=0.7, color='green', edgecolor='darkgreen')
ax.axvline(np.median(mean_dists), color='red', ls='--', lw=2, label=f'Median = {np.median(mean_dists):.2f} kpc')
ax.set_xlabel('Distance (kpc)')
ax.set_ylabel('Number of Sources')
ax.set_title('Distance Distribution')
ax.set_xlim(0, np.percentile(mean_dists, 99))
ax.legend()
ax.grid(True, alpha=0.3)

# Panel 3: Extinction distribution
ax = axes[0, 2]
ax.hist(mean_av, bins=30, alpha=0.7, color='brown', edgecolor='darkred')
ax.axvline(np.median(mean_av), color='blue', ls='--', lw=2, label=f'Median = {np.median(mean_av):.2f} mag')
ax.set_xlabel('A(V) (mag)')
ax.set_ylabel('Number of Sources')
ax.set_title('Extinction Distribution')
ax.set_xlim(0, np.percentile(mean_av, 99))
ax.legend()
ax.grid(True, alpha=0.3)

# Panel 4: Distance vs extinction
ax = axes[1, 0]
scatter = ax.scatter(mean_dists, mean_av, alpha=0.5, s=20,
                    c=chi2_per_band, cmap='viridis', vmin=0, vmax=3)
ax.set_xlabel('Distance (kpc)')
ax.set_ylabel('A(V) (mag)')
ax.set_title('Distance-Extinction Relation')
ax.set_xlim(0, np.percentile(mean_dists, 99))
ax.set_ylim(0, np.percentile(mean_av, 99))
plt.colorbar(scatter, ax=ax, label='chi2/band')
ax.grid(True, alpha=0.3)

# Panel 5: Stellar parameters (using best-fit model per source)
ax = axes[1, 1]

best_draw = np.argmax(log_post, axis=1)
best_model_per_source = model_idx[np.arange(len(chi2)), best_draw]
unique_masses = []
for idx in best_model_per_source:
    if idx < len(labels):
        unique_masses.append(labels[idx]['mini'])

if unique_masses:
    ax.hist(unique_masses, bins=np.logspace(-0.5, 1, 30),
           alpha=0.7, color='purple', edgecolor='darkviolet')
    ax.set_xlabel('Initial Mass (M_sun)')
    ax.set_ylabel('Number of Sources')
    ax.set_title('Mass Distribution')
    ax.set_xscale('log')
    ax.grid(True, alpha=0.3)

# Panel 6: Fit statistics
ax = axes[1, 2]
ax.axis('off')

# Calculate statistics
good_fits = pvals > 1e-3
stats_text = f"""
Fit Statistics:

Total sources: {len(chi2)}
Good fits (p > 1e-3): {good_fits.sum()} ({100*good_fits.sum()/len(chi2):.1f}%)

Distance:
  Mean: {np.mean(mean_dists):.2f} kpc
  Median: {np.median(mean_dists):.2f} kpc
  Std: {np.std(mean_dists):.2f} kpc

Extinction:
  Mean: {np.mean(mean_av):.2f} mag
  Median: {np.median(mean_av):.2f} mag
  Std: {np.std(mean_av):.2f} mag

chi2/band:
  Mean: {np.mean(chi2_per_band):.2f}
  Median: {np.median(chi2_per_band):.2f}

R(V):
  Mean: {np.mean(dreds):.2f}
  Std: {np.std(dreds):.2f}
"""

ax.text(0.05, 0.95, stats_text, transform=ax.transAxes,
       fontsize=11, va='top', family='monospace')

plt.suptitle('BruteForce Fitting Results', fontsize=16, fontweight='bold')
save_figure(fig, 'fitting_results')
plt.show()

print("\nResult analysis complete")

Loaded fitting results:
  chi2 values: (207,)
  Model indices: (207, 500)
  ML scale/AV/RV: (207, 500)
  Distance samples: (207, 500)

  Saved: /home/user/brutus/tutorials/plots/tutorial_05/fitting_results.png
../_images/3ac40491c726a4083bc971af27b554662a00e10e38bc30dc02047defb2c47182.png 
Result analysis complete

4. Visualizing Individual Results#

Let’s examine individual star fits in detail to understand the posterior distributions and parameter correlations.

# Select a typical example star that passes the goodness-of-fit p-value cut
good_fit = pvals > 1e-3  # p-value threshold for a good fit

# Use the first source that passes the p-value cut
if good_fit[0]:
    star_idx = 0
else:
    star_idx = np.where(good_fit)[0][0]

data_idx = idx_fit[star_idx]

print(f"Selected example star (index {star_idx}):")
print(f"  p-value = {pvals[star_idx]:.3g}")
print(f"  N_bands = {nbands[star_idx]}")
if np.isfinite(data['parallax'][data_idx]):
    plx_snr = data['parallax'][data_idx] / data['parallax_err'][data_idx]
    print(f"  Parallax = {data['parallax'][data_idx]:.3f} +/- {data['parallax_err'][data_idx]:.3f} mas (SNR = {plx_snr:.1f})")
else:
    print(f"  Parallax = not available")
print(f"  Distance = {mean_dists[star_idx]:.2f} +/- {np.std(dists[star_idx]):.3f} kpc")
print(f"  A(V) = {mean_av[star_idx]:.2f} mag")

Selected example star (index 0):
  p-value = 0.767
  N_bands = 9
  Parallax = 3.495 +/- 0.048 mas (SNR = 72.9)
  Distance = 0.28 +/- 0.005 kpc
  A(V) = 0.26 mag

4.1 Visualizing a Single Star#

brutus provides built-in plotting functions for visualizing individual fitting results. We demonstrate three key plots for the selected example star.

from brutus.plotting import dist_vs_red, posterior_predictive, cornerplot

print(f"Visualizing star {star_idx}")
print(f"  p-value = {pvals[star_idx]:.3g}")
print(f"  Distance = {mean_dists[star_idx]:.2f} +/- {np.std(dists[star_idx]):.3f} kpc")
print(f"  A(V) = {mean_av[star_idx]:.2f} mag")

Visualizing star 0
  p-value = 0.767
  Distance = 0.28 +/- 0.005 kpc
  A(V) = 0.26 mag

Hide code cell source

 
# Distance vs reddening posterior
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
plt.sca(ax)

dist_vs_red(
    (dists[star_idx], reds[star_idx], dreds[star_idx]),
    parallax=data['parallax'][data_idx],
    parallax_err=data['parallax_err'][data_idx],
    coord=(data['coords'][data_idx, 0], data['coords'][data_idx, 1]),
)

plt.title('Distance vs Reddening Posterior')
save_figure(fig, 'dist_vs_red')
plt.show()

  Saved: /home/user/brutus/tutorials/plots/tutorial_05/dist_vs_red.png
../_images/3b708e98004cd1ad681e0ec457f80a4fde9c471498a0f2df2f730475e4a232ae.png

Hide code cell source

 
# Posterior predictive SED
fig, ax, parts = posterior_predictive(
    models,
    idxs=model_idx[star_idx],
    reds=reds[star_idx],
    dreds=dreds[star_idx],
    dists=dists[star_idx],
    data=data['phot'][data_idx],
    data_err=data['err'][data_idx],
    data_mask=data['mask'][data_idx],
    offset=offsets_mist,
    labels=filt,
)

plt.title('Posterior Predictive SED')
save_figure(fig, 'posterior_predictive')
plt.show()

  Saved: /home/user/brutus/tutorials/plots/tutorial_05/posterior_predictive.png
../_images/a57f1131f9bbf2c39bfaa37e966d24fe9ee700b9068d7be7654ebc38ca9634a2.png

Hide code cell source

 
# Corner plot of stellar parameters
fig, axes = cornerplot(
    idxs=model_idx[star_idx],
    data=(dists[star_idx], reds[star_idx], dreds[star_idx]),
    params=labels,
    coord=(data['coords'][data_idx, 0], data['coords'][data_idx, 1]),
    parallax=data['parallax'][data_idx],
    parallax_err=data['parallax_err'][data_idx],
    show_titles=True,
)

# Add padding to prevent axis label overlap with tick numbers
fig.subplots_adjust(bottom=0.15, left=0.15, hspace=0.1, wspace=0.1)

save_figure(fig, 'cornerplot')
plt.show()

  Saved: /home/user/brutus/tutorials/plots/tutorial_05/cornerplot.png
../_images/b7f5aed136d2db887252207a6a45e3e69fffd5e5ed2c847f8dc43ca40af66cf3.png

5. Impact of the 3D Dust Map Prior#

The Bayestar 3D dust map provides distance-dependent extinction constraints that can significantly affect the inferred stellar parameters. When a dust map is available, BruteForce evaluates the extinction prior at each candidate distance, preferring (distance, A_V) combinations consistent with the cumulative dust profile along the line of sight.

To quantify this effect, we re-fit the same sources without the dust map prior and compare the resulting distance and extinction posteriors. Sources in highly extincted regions (like Orion) should show the largest differences.

Hide code cell source

 
# Re-fit the same sources WITHOUT the dust map prior
if dustfile is not None:
    output_file_nodust = tutorials_dir / 'Orion_l209.1_b-19.9_mist_nodust.h5'

    print(f"Running comparison fit WITHOUT dust map prior on {len(idx_fit)} sources...")
    print("This may take several minutes depending on your hardware.")

    t0 = time.time()

    bf.fit(
        data=data['phot'][idx_fit],
        data_err=data['err'][idx_fit],
        data_mask=data['mask'][idx_fit],
        data_labels=idx_fit,
        save_file=str(output_file_nodust),
        data_coords=data['coords'][idx_fit],
        parallax=data['parallax'][idx_fit],
        parallax_err=data['parallax_err'][idx_fit],
        phot_offsets=offsets_mist,
        dustfile=None,  # <-- No dust map
        Ndraws=500,
        Nmc_prior=20,
        logl_dim_prior=True,
        save_dar_draws=True,
        running_io=False,
        verbose=True
    )

    elapsed = time.time() - t0
    print(f"\nComparison fit complete in {elapsed:.1f} seconds ({elapsed/60:.1f} minutes).")
    print(f"Saved to {output_file_nodust}")

    # --- Load both result sets ---
    with h5py.File(output_file, 'r') as f:
        reds_dust = f['samps_red'][:]

    with h5py.File(output_file_nodust, 'r') as f:
        reds_nodust = f['samps_red'][:]

    mean_av_dust = np.mean(reds_dust, axis=1)
    mean_av_nodust = np.mean(reds_nodust, axis=1)

    # --- Single-panel A(V) comparison ---
    fig, ax = plt.subplots(figsize=(6, 6))
    ax.scatter(mean_av_nodust, mean_av_dust, s=30, alpha=0.7,
               c='firebrick', edgecolor='k', lw=0.5)
    ax.plot([0, 6], [0, 6], 'k--', lw=1, alpha=0.5, label='1:1')
    ax.set_xlabel('A(V) (mag) -- No Dust Prior')
    ax.set_ylabel('A(V) (mag) -- With Dust Prior')
    ax.set_title('Effect of 3D Dust Map Prior on Extinction')
    ax.set_xlim(0, 6)
    ax.set_ylim(0, 6)
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_aspect('equal')

    save_figure(fig, 'dust_prior_comparison')
    plt.show()

    print("\nDust prior comparison complete")
else:
    print("Dust map not available -- skipping comparison.")
    print("The dust prior comparison requires bayestar2019_v1.h5")

Running comparison fit WITHOUT dust map prior on 207 sources...
This may take several minutes depending on your hardware.

Fitting object 1/207  

Fitting object 2/207 [chi2/n: 3.3/9] (mean time: 0.452 s/obj, est. remaining: 93.164 s)    

Fitting object 3/207 [chi2/n: 50.5/9] (mean time: 0.461 s/obj, est. remaining: 94.453 s)    

Fitting object 4/207 [chi2/n: 44.0/9] (mean time: 0.471 s/obj, est. remaining: 96.110 s)    

Fitting object 5/207 [chi2/n: 1.8/9] (mean time: 0.469 s/obj, est. remaining: 95.142 s)    

Fitting object 6/207 [chi2/n: 44.7/9] (mean time: 0.472 s/obj, est. remaining: 95.392 s)    

Fitting object 7/207 [chi2/n: 60.2/9] (mean time: 0.469 s/obj, est. remaining: 94.256 s)    

Fitting object 8/207 [chi2/n: 51.4/9] (mean time: 0.487 s/obj, est. remaining: 97.419 s)    

Fitting object 9/207 [chi2/n: 30.8/9] (mean time: 0.493 s/obj, est. remaining: 98.080 s)    

Fitting object 10/207 [chi2/n: 25.2/9] (mean time: 0.492 s/obj, est. remaining: 97.362 s)    

Fitting object 11/207 [chi2/n: 12.6/9] (mean time: 0.533 s/obj, est. remaining: 104.907 s)    

Fitting object 12/207 [chi2/n: 12.0/8] (mean time: 0.556 s/obj, est. remaining: 108.968 s)    

Fitting object 13/207 [chi2/n: 6.0/6] (mean time: 0.544 s/obj, est. remaining: 106.053 s)    

Fitting object 14/207 [chi2/n: 48.3/9] (mean time: 0.539 s/obj, est. remaining: 104.521 s)    

Fitting object 15/207 [chi2/n: 99.2/9] (mean time: 0.536 s/obj, est. remaining: 103.411 s)    

Fitting object 16/207 [chi2/n: 140.2/9] (mean time: 0.537 s/obj, est. remaining: 103.043 s)    

Fitting object 17/207 [chi2/n: 86.1/9] (mean time: 0.548 s/obj, est. remaining: 104.710 s)    

Fitting object 18/207 [chi2/n: 39.3/9] (mean time: 0.552 s/obj, est. remaining: 104.842 s)    

Fitting object 19/207 [chi2/n: 62.8/8] (mean time: 0.555 s/obj, est. remaining: 104.955 s)    

Fitting object 20/207 [chi2/n: 55.6/8] (mean time: 0.558 s/obj, est. remaining: 104.944 s)    

Fitting object 21/207 [chi2/n: 47.3/5] (mean time: 0.565 s/obj, est. remaining: 105.577 s)    

Fitting object 22/207 [chi2/n: 43.3/9] (mean time: 0.566 s/obj, est. remaining: 105.202 s)    

Fitting object 23/207 [chi2/n: 34.6/8] (mean time: 0.563 s/obj, est. remaining: 104.117 s)    

Fitting object 24/207 [chi2/n: 267.1/9] (mean time: 0.564 s/obj, est. remaining: 103.766 s)    

Fitting object 25/207 [chi2/n: 4.6/9] (mean time: 0.563 s/obj, est. remaining: 103.087 s)    

Fitting object 26/207 [chi2/n: 210.0/9] (mean time: 0.565 s/obj, est. remaining: 102.882 s)    

Fitting object 27/207 [chi2/n: 192.6/9] (mean time: 0.565 s/obj, est. remaining: 102.177 s)    

Fitting object 28/207 [chi2/n: 113.3/9] (mean time: 0.565 s/obj, est. remaining: 101.622 s)    

Fitting object 29/207 [chi2/n: 74.9/8] (mean time: 0.563 s/obj, est. remaining: 100.714 s)    

Fitting object 30/207 [chi2/n: 24.8/9] (mean time: 0.561 s/obj, est. remaining: 99.934 s)    

Fitting object 31/207 [chi2/n: 15.9/9] (mean time: 0.573 s/obj, est. remaining: 101.373 s)    

Fitting object 32/207 [chi2/n: 18.9/9] (mean time: 0.574 s/obj, est. remaining: 101.030 s)    

Fitting object 33/207 [chi2/n: 93.8/9] (mean time: 0.575 s/obj, est. remaining: 100.554 s)    

Fitting object 34/207 [chi2/n: 41.6/9] (mean time: 0.574 s/obj, est. remaining: 99.831 s)    

Fitting object 35/207 [chi2/n: 72.4/9] (mean time: 0.572 s/obj, est. remaining: 98.955 s)    

Fitting object 36/207 [chi2/n: 13.1/9] (mean time: 0.575 s/obj, est. remaining: 98.901 s)    

Fitting object 37/207 [chi2/n: 42.7/9] (mean time: 0.575 s/obj, est. remaining: 98.260 s)    

Fitting object 38/207 [chi2/n: 16.5/9] (mean time: 0.589 s/obj, est. remaining: 100.057 s)    

Fitting object 39/207 [chi2/n: 11.2/8] (mean time: 0.594 s/obj, est. remaining: 100.408 s)    

Fitting object 40/207 [chi2/n: 76.0/9] (mean time: 0.594 s/obj, est. remaining: 99.738 s)    

Fitting object 41/207 [chi2/n: 17.8/9] (mean time: 0.607 s/obj, est. remaining: 101.420 s)    

Fitting object 42/207 [chi2/n: 5.2/8] (mean time: 0.617 s/obj, est. remaining: 102.412 s)    

Fitting object 43/207 [chi2/n: 8.6/8] (mean time: 0.633 s/obj, est. remaining: 104.419 s)    

Fitting object 44/207 [chi2/n: 9.0/9] (mean time: 0.639 s/obj, est. remaining: 104.767 s)    

Fitting object 45/207 [chi2/n: 2.9/8] (mean time: 0.643 s/obj, est. remaining: 104.810 s)    

Fitting object 46/207 [chi2/n: 54.3/9] (mean time: 0.642 s/obj, est. remaining: 104.057 s)    

Fitting object 47/207 [chi2/n: 11.4/8] (mean time: 0.649 s/obj, est. remaining: 104.558 s)    

Fitting object 48/207 [chi2/n: 8.8/8] (mean time: 0.657 s/obj, est. remaining: 105.079 s)    

Fitting object 49/207 [chi2/n: 6.3/9] (mean time: 0.667 s/obj, est. remaining: 105.981 s)    

Fitting object 50/207 [chi2/n: 2.2/7] (mean time: 0.670 s/obj, est. remaining: 105.881 s)    

Fitting object 51/207 [chi2/n: 21.0/9] (mean time: 0.683 s/obj, est. remaining: 107.179 s)    

Fitting object 52/207 [chi2/n: 9.9/8] (mean time: 0.689 s/obj, est. remaining: 107.450 s)    

Fitting object 53/207 [chi2/n: 5.2/8] (mean time: 0.696 s/obj, est. remaining: 107.840 s)    

Fitting object 54/207 [chi2/n: 8.0/9] (mean time: 0.703 s/obj, est. remaining: 108.255 s)    

Fitting object 55/207 [chi2/n: 5.2/6] (mean time: 0.702 s/obj, est. remaining: 107.406 s)    

Fitting object 56/207 [chi2/n: 2.5/6] (mean time: 0.714 s/obj, est. remaining: 108.459 s)    

Fitting object 57/207 [chi2/n: 1.1/6] (mean time: 0.713 s/obj, est. remaining: 107.628 s)    

Fitting object 58/207 [chi2/n: 0.5/6] (mean time: 0.714 s/obj, est. remaining: 107.066 s)    

Fitting object 59/207 [chi2/n: 0.2/6] (mean time: 0.718 s/obj, est. remaining: 106.971 s)    

Fitting object 60/207 [chi2/n: 25.0/9] (mean time: 0.714 s/obj, est. remaining: 105.669 s)    

Fitting object 61/207 [chi2/n: 0.3/6] (mean time: 0.714 s/obj, est. remaining: 104.940 s)    

Fitting object 62/207 [chi2/n: 27.5/9] (mean time: 0.711 s/obj, est. remaining: 103.757 s)    

Fitting object 63/207 [chi2/n: 21.8/9] (mean time: 0.707 s/obj, est. remaining: 102.472 s)    

Fitting object 64/207 [chi2/n: 8.3/9] (mean time: 0.703 s/obj, est. remaining: 101.270 s)    

Fitting object 65/207 [chi2/n: 67.4/9] (mean time: 0.700 s/obj, est. remaining: 100.097 s)    

Fitting object 66/207 [chi2/n: 23.0/9] (mean time: 0.701 s/obj, est. remaining: 99.504 s)    

Fitting object 67/207 [chi2/n: 14.4/9] (mean time: 0.705 s/obj, est. remaining: 99.339 s)    

Fitting object 68/207 [chi2/n: 8.1/6] (mean time: 0.703 s/obj, est. remaining: 98.405 s)    

Fitting object 69/207 [chi2/n: 6.8/9] (mean time: 0.722 s/obj, est. remaining: 100.339 s)    

Fitting object 70/207 [chi2/n: 407.0/9] (mean time: 0.721 s/obj, est. remaining: 99.435 s)    

Fitting object 71/207 [chi2/n: 3146.0/9] (mean time: 0.719 s/obj, est. remaining: 98.494 s)    

Fitting object 72/207 [chi2/n: 1537.7/9] (mean time: 0.718 s/obj, est. remaining: 97.629 s)    

Fitting object 73/207 [chi2/n: 1473.8/9] (mean time: 0.716 s/obj, est. remaining: 96.626 s)    

Fitting object 74/207 [chi2/n: 884.2/9] (mean time: 0.713 s/obj, est. remaining: 95.586 s)    

Fitting object 75/207 [chi2/n: 708.9/9] (mean time: 0.711 s/obj, est. remaining: 94.575 s)    

Fitting object 76/207 [chi2/n: 66.8/9] (mean time: 0.709 s/obj, est. remaining: 93.562 s)    

Fitting object 77/207 [chi2/n: 1808.9/9] (mean time: 0.707 s/obj, est. remaining: 92.630 s)    

Fitting object 78/207 [chi2/n: 176.8/9] (mean time: 0.705 s/obj, est. remaining: 91.587 s)    

Fitting object 79/207 [chi2/n: 9.5/9] (mean time: 0.704 s/obj, est. remaining: 90.862 s)    

Fitting object 80/207 [chi2/n: 26.2/9] (mean time: 0.702 s/obj, est. remaining: 89.892 s)    

Fitting object 81/207 [chi2/n: 72.2/9] (mean time: 0.701 s/obj, est. remaining: 89.032 s)    

Fitting object 82/207 [chi2/n: 797.4/9] (mean time: 0.699 s/obj, est. remaining: 88.126 s)    

Fitting object 83/207 [chi2/n: 0.2/6] (mean time: 0.696 s/obj, est. remaining: 87.051 s)    

Fitting object 84/207 [chi2/n: 1553.4/8] (mean time: 0.695 s/obj, est. remaining: 86.136 s)    

Fitting object 85/207 [chi2/n: 20.6/9] (mean time: 0.693 s/obj, est. remaining: 85.271 s)    

Fitting object 86/207 [chi2/n: 40.5/9] (mean time: 0.694 s/obj, est. remaining: 84.629 s)    

Fitting object 87/207 [chi2/n: 148.4/8] (mean time: 0.696 s/obj, est. remaining: 84.196 s)    

Fitting object 88/207 [chi2/n: 1751.2/9] (mean time: 0.694 s/obj, est. remaining: 83.310 s)    

Fitting object 89/207 [chi2/n: 16.5/9] (mean time: 0.700 s/obj, est. remaining: 83.319 s)    

Fitting object 90/207 [chi2/n: 25.8/8] (mean time: 0.698 s/obj, est. remaining: 82.375 s)    

Fitting object 91/207 [chi2/n: 15.0/7] (mean time: 0.703 s/obj, est. remaining: 82.215 s)    

Fitting object 92/207 [chi2/n: 37.5/4] (mean time: 0.700 s/obj, est. remaining: 81.203 s)    

Fitting object 93/207 [chi2/n: 968.9/9] (mean time: 0.699 s/obj, est. remaining: 80.352 s)    

Fitting object 94/207 [chi2/n: 657.6/9] (mean time: 0.697 s/obj, est. remaining: 79.441 s)    

Fitting object 95/207 [chi2/n: 106.8/8] (mean time: 0.695 s/obj, est. remaining: 78.549 s)    

Fitting object 96/207 [chi2/n: 30.5/9] (mean time: 0.694 s/obj, est. remaining: 77.714 s)    

Fitting object 97/207 [chi2/n: 252.8/8] (mean time: 0.692 s/obj, est. remaining: 76.849 s)    

Fitting object 98/207 [chi2/n: 161.7/9] (mean time: 0.691 s/obj, est. remaining: 76.062 s)    

Fitting object 99/207 [chi2/n: 29.2/9] (mean time: 0.690 s/obj, est. remaining: 75.221 s)    

Fitting object 100/207 [chi2/n: 405.7/9] (mean time: 0.689 s/obj, est. remaining: 74.409 s)    

Fitting object 101/207 [chi2/n: 165.0/6] (mean time: 0.686 s/obj, est. remaining: 73.447 s)    

Fitting object 102/207 [chi2/n: 16.9/9] (mean time: 0.685 s/obj, est. remaining: 72.610 s)    

Fitting object 103/207 [chi2/n: 8.4/9] (mean time: 0.685 s/obj, est. remaining: 71.899 s)    

Fitting object 104/207 [chi2/n: 106.1/9] (mean time: 0.684 s/obj, est. remaining: 71.090 s)    

Fitting object 105/207 [chi2/n: 11.4/8] (mean time: 0.682 s/obj, est. remaining: 70.294 s)    

Fitting object 106/207 [chi2/n: 71.5/9] (mean time: 0.682 s/obj, est. remaining: 69.567 s)    

Fitting object 107/207 [chi2/n: 30.0/9] (mean time: 0.680 s/obj, est. remaining: 68.723 s)    

Fitting object 108/207 [chi2/n: 28.2/9] (mean time: 0.679 s/obj, est. remaining: 67.936 s)    

Fitting object 109/207 [chi2/n: 9.1/9] (mean time: 0.678 s/obj, est. remaining: 67.139 s)    

Fitting object 110/207 [chi2/n: 17.1/9] (mean time: 0.681 s/obj, est. remaining: 66.770 s)    

Fitting object 111/207 [chi2/n: 392.7/9] (mean time: 0.679 s/obj, est. remaining: 65.908 s)    

Fitting object 112/207 [chi2/n: 17.2/9] (mean time: 0.683 s/obj, est. remaining: 65.549 s)    

Fitting object 113/207 [chi2/n: 1118.2/9] (mean time: 0.682 s/obj, est. remaining: 64.744 s)    

Fitting object 114/207 [chi2/n: 17.1/7] (mean time: 0.684 s/obj, est. remaining: 64.254 s)    

Fitting object 115/207 [chi2/n: 34.4/6] (mean time: 0.681 s/obj, est. remaining: 63.339 s)    

Fitting object 116/207 [chi2/n: 3018.7/9] (mean time: 0.679 s/obj, est. remaining: 62.498 s)    

Fitting object 117/207 [chi2/n: 0.4/9] (mean time: 0.678 s/obj, est. remaining: 61.672 s)    

Fitting object 118/207 [chi2/n: 82.2/5] (mean time: 0.676 s/obj, est. remaining: 60.807 s)    

Fitting object 119/207 [chi2/n: 13.4/9] (mean time: 0.674 s/obj, est. remaining: 60.020 s)    

Fitting object 120/207 [chi2/n: 24.7/9] (mean time: 0.677 s/obj, est. remaining: 59.540 s)    

Fitting object 121/207 [chi2/n: 33.0/9] (mean time: 0.675 s/obj, est. remaining: 58.742 s)    

Fitting object 122/207 [chi2/n: 29.9/9] (mean time: 0.676 s/obj, est. remaining: 58.124 s)    

Fitting object 123/207 [chi2/n: 107.3/7] (mean time: 0.674 s/obj, est. remaining: 57.286 s)    

Fitting object 124/207 [chi2/n: 27.9/9] (mean time: 0.673 s/obj, est. remaining: 56.538 s)    

Fitting object 125/207 [chi2/n: 1.5/9] (mean time: 0.672 s/obj, est. remaining: 55.776 s)    

Fitting object 126/207 [chi2/n: 1895.5/9] (mean time: 0.671 s/obj, est. remaining: 54.989 s)    

Fitting object 127/207 [chi2/n: 24.5/9] (mean time: 0.669 s/obj, est. remaining: 54.194 s)    

Fitting object 128/207 [chi2/n: 23.3/9] (mean time: 0.668 s/obj, est. remaining: 53.418 s)    

Fitting object 129/207 [chi2/n: 32.3/9] (mean time: 0.666 s/obj, est. remaining: 52.620 s)    

Fitting object 130/207 [chi2/n: 34.5/9] (mean time: 0.665 s/obj, est. remaining: 51.834 s)    

Fitting object 131/207 [chi2/n: 174.1/9] (mean time: 0.664 s/obj, est. remaining: 51.091 s)    

Fitting object 132/207 [chi2/n: 17.0/9] (mean time: 0.663 s/obj, est. remaining: 50.363 s)    

Fitting object 133/207 [chi2/n: 26.5/9] (mean time: 0.661 s/obj, est. remaining: 49.594 s)    

Fitting object 134/207 [chi2/n: 32.9/9] (mean time: 0.660 s/obj, est. remaining: 48.863 s)    

Fitting object 135/207 [chi2/n: 17.9/9] (mean time: 0.659 s/obj, est. remaining: 48.106 s)    

Fitting object 136/207 [chi2/n: 74.7/9] (mean time: 0.658 s/obj, est. remaining: 47.371 s)    

Fitting object 137/207 [chi2/n: 24.0/9] (mean time: 0.657 s/obj, est. remaining: 46.624 s)    

Fitting object 138/207 [chi2/n: 148.2/9] (mean time: 0.656 s/obj, est. remaining: 45.885 s)    

Fitting object 139/207 [chi2/n: 16.4/9] (mean time: 0.654 s/obj, est. remaining: 45.148 s)    

Fitting object 140/207 [chi2/n: 32.3/9] (mean time: 0.653 s/obj, est. remaining: 44.407 s)    

Fitting object 141/207 [chi2/n: 15.9/9] (mean time: 0.652 s/obj, est. remaining: 43.695 s)    

Fitting object 142/207 [chi2/n: 34.5/9] (mean time: 0.653 s/obj, est. remaining: 43.092 s)    

Fitting object 143/207 [chi2/n: 38.0/8] (mean time: 0.655 s/obj, est. remaining: 42.550 s)    

Fitting object 144/207 [chi2/n: 48.7/9] (mean time: 0.654 s/obj, est. remaining: 41.847 s)    

Fitting object 145/207 [chi2/n: 18.8/9] (mean time: 0.653 s/obj, est. remaining: 41.129 s)    

Fitting object 146/207 [chi2/n: 32.9/9] (mean time: 0.652 s/obj, est. remaining: 40.406 s)    

Fitting object 147/207 [chi2/n: 15.9/9] (mean time: 0.651 s/obj, est. remaining: 39.702 s)    

Fitting object 148/207 [chi2/n: 22.0/9] (mean time: 0.650 s/obj, est. remaining: 39.009 s)    

Fitting object 149/207 [chi2/n: 3.9/9] (mean time: 0.654 s/obj, est. remaining: 38.561 s)    

Fitting object 150/207 [chi2/n: 109.9/9] (mean time: 0.653 s/obj, est. remaining: 37.882 s)    

Fitting object 151/207 [chi2/n: 21.6/9] (mean time: 0.653 s/obj, est. remaining: 37.213 s)    

Fitting object 152/207 [chi2/n: 17.1/9] (mean time: 0.652 s/obj, est. remaining: 36.523 s)    

Fitting object 153/207 [chi2/n: 21.8/9] (mean time: 0.651 s/obj, est. remaining: 35.827 s)    

Fitting object 154/207 [chi2/n: 106.0/9] (mean time: 0.651 s/obj, est. remaining: 35.129 s)    

Fitting object 155/207 [chi2/n: 8.9/9] (mean time: 0.654 s/obj, est. remaining: 34.648 s)    

Fitting object 156/207 [chi2/n: 23.1/9] (mean time: 0.653 s/obj, est. remaining: 33.942 s)    

Fitting object 157/207 [chi2/n: 11.1/8] (mean time: 0.654 s/obj, est. remaining: 33.331 s)    

Fitting object 158/207 [chi2/n: 2.7/8] (mean time: 0.655 s/obj, est. remaining: 32.731 s)    

Fitting object 159/207 [chi2/n: 11.6/8] (mean time: 0.657 s/obj, est. remaining: 32.187 s)    

Fitting object 160/207 [chi2/n: 21.2/9] (mean time: 0.656 s/obj, est. remaining: 31.500 s)    

Fitting object 161/207 [chi2/n: 15.9/9] (mean time: 0.662 s/obj, est. remaining: 31.092 s)    

Fitting object 162/207 [chi2/n: 16.5/9] (mean time: 0.663 s/obj, est. remaining: 30.483 s)    

Fitting object 163/207 [chi2/n: 11.8/8] (mean time: 0.662 s/obj, est. remaining: 29.773 s)    

Fitting object 164/207 [chi2/n: 0.9/7] (mean time: 0.663 s/obj, est. remaining: 29.150 s)    

Fitting object 165/207 [chi2/n: 9.1/8] (mean time: 0.664 s/obj, est. remaining: 28.543 s)    

Fitting object 166/207 [chi2/n: 3.0/6] (mean time: 0.664 s/obj, est. remaining: 27.892 s)    

Fitting object 167/207 [chi2/n: 28.8/8] (mean time: 0.666 s/obj, est. remaining: 27.292 s)    

Fitting object 168/207 [chi2/n: 0.7/6] (mean time: 0.665 s/obj, est. remaining: 26.609 s)    

Fitting object 169/207 [chi2/n: 1.6/6] (mean time: 0.666 s/obj, est. remaining: 25.991 s)    

Fitting object 170/207 [chi2/n: 0.2/6] (mean time: 0.667 s/obj, est. remaining: 25.345 s)    

Fitting object 171/207 [chi2/n: 1.3/6] (mean time: 0.670 s/obj, est. remaining: 24.799 s)    

Fitting object 172/207 [chi2/n: 1.5/6] (mean time: 0.671 s/obj, est. remaining: 24.148 s)    

Fitting object 173/207 [chi2/n: 0.7/6] (mean time: 0.674 s/obj, est. remaining: 23.587 s)    

Fitting object 174/207 [chi2/n: 1.8/6] (mean time: 0.674 s/obj, est. remaining: 22.900 s)    

Fitting object 175/207 [chi2/n: 1.1/6] (mean time: 0.673 s/obj, est. remaining: 22.221 s)    

Fitting object 176/207 [chi2/n: 1.9/6] (mean time: 0.673 s/obj, est. remaining: 21.545 s)    

Fitting object 177/207 [chi2/n: 0.8/6] (mean time: 0.674 s/obj, est. remaining: 20.879 s)    

Fitting object 178/207 [chi2/n: 0.7/6] (mean time: 0.674 s/obj, est. remaining: 20.228 s)    

Fitting object 179/207 [chi2/n: 10.4/9] (mean time: 0.673 s/obj, est. remaining: 19.527 s)    

Fitting object 180/207 [chi2/n: 10.4/9] (mean time: 0.673 s/obj, est. remaining: 18.835 s)    

Fitting object 181/207 [chi2/n: 9.4/9] (mean time: 0.672 s/obj, est. remaining: 18.139 s)    

Fitting object 182/207 [chi2/n: 16.4/9] (mean time: 0.671 s/obj, est. remaining: 17.445 s)    

Fitting object 183/207 [chi2/n: 7.2/9] (mean time: 0.670 s/obj, est. remaining: 16.753 s)    

Fitting object 184/207 [chi2/n: 59.0/9] (mean time: 0.669 s/obj, est. remaining: 16.061 s)    

Fitting object 185/207 [chi2/n: 9.6/9] (mean time: 0.669 s/obj, est. remaining: 15.384 s)    

Fitting object 186/207 [chi2/n: 22.0/9] (mean time: 0.668 s/obj, est. remaining: 14.697 s)    

Fitting object 187/207 [chi2/n: 25.4/9] (mean time: 0.672 s/obj, est. remaining: 14.119 s)    

Fitting object 188/207 [chi2/n: 16.5/6] (mean time: 0.671 s/obj, est. remaining: 13.420 s)    

Fitting object 189/207 [chi2/n: 39.7/9] (mean time: 0.671 s/obj, est. remaining: 12.740 s)    

Fitting object 190/207 [chi2/n: 9.3/9] (mean time: 0.672 s/obj, est. remaining: 12.091 s)    

Fitting object 191/207 [chi2/n: 3.9/9] (mean time: 0.675 s/obj, est. remaining: 11.471 s)    

Fitting object 192/207 [chi2/n: 11.7/9] (mean time: 0.674 s/obj, est. remaining: 10.785 s)    

Fitting object 193/207 [chi2/n: 13.3/9] (mean time: 0.675 s/obj, est. remaining: 10.127 s)    

Fitting object 194/207 [chi2/n: 16.8/9] (mean time: 0.675 s/obj, est. remaining: 9.443 s)    

Fitting object 195/207 [chi2/n: 8.4/9] (mean time: 0.678 s/obj, est. remaining: 8.811 s)    

Fitting object 196/207 [chi2/n: 19.5/6] (mean time: 0.677 s/obj, est. remaining: 8.126 s)    

Fitting object 197/207 [chi2/n: 2.5/6] (mean time: 0.678 s/obj, est. remaining: 7.456 s)    

Fitting object 198/207 [chi2/n: 10.6/9] (mean time: 0.678 s/obj, est. remaining: 6.783 s)    

Fitting object 199/207 [chi2/n: 30.8/8] (mean time: 0.678 s/obj, est. remaining: 6.101 s)    

Fitting object 200/207 [chi2/n: 13.4/9] (mean time: 0.680 s/obj, est. remaining: 5.442 s)    

Fitting object 201/207 [chi2/n: 4.3/6] (mean time: 0.679 s/obj, est. remaining: 4.753 s)    

Fitting object 202/207 [chi2/n: 76.3/9] (mean time: 0.678 s/obj, est. remaining: 4.070 s)    

Fitting object 203/207 [chi2/n: 20.9/9] (mean time: 0.678 s/obj, est. remaining: 3.388 s)    

Fitting object 204/207 [chi2/n: 13.2/9] (mean time: 0.681 s/obj, est. remaining: 2.722 s)    

Fitting object 205/207 [chi2/n: 114.8/8] (mean time: 0.680 s/obj, est. remaining: 2.039 s)    

Fitting object 206/207 [chi2/n: 3.3/6] (mean time: 0.681 s/obj, est. remaining: 1.361 s)    

Fitting object 207/207 [chi2/n: 15.6/6] (mean time: 0.680 s/obj, est. remaining: 0.680 s)    

Fitting object 207/207 [chi2/n: 19.3/9] (mean time: 0.684 s/obj, total: 141.581 s)    

Comparison fit complete in 142.3 seconds (2.4 minutes).
Saved to Orion_l209.1_b-19.9_mist_nodust.h5

  Saved: /home/user/brutus/tutorials/plots/tutorial_05/dust_prior_comparison.png
../_images/8998ba5bba5410a3d482a5d03e0e789abd7f12d09acc7bfec2fcfbf1153af559.png 
Dust prior comparison complete

6. Advanced BruteForce Options#

BruteForce.fit() accepts many parameters beyond the basics shown earlier. Understanding these options lets you tune performance, control output size, and adjust the statistical model to match your science case.

The parameters fall into several categories:

  • Likelihood control: logl_dim_prior switches between a standard Gaussian log-likelihood and a chi-square dimensionality-corrected version that penalizes over-fitting (or under-fitting) relative to the number of observed bands.

  • Prior adjustments: apply_agewt adds an age-based weight from the MIST evolutionary tracks (accounting for time spent in each evolutionary phase). apply_grad applies a gradient-based correction for non-uniform grid spacing. av_gauss imposes a Gaussian prior on A(V) during the per-model likelihood optimization. Note that it is applied in addition to (not instead of) any dust-map prior from dustfile, so passing an informative av_gauss together with a dust map counts the A(V) constraint twice.

  • Output thresholds: wt_thresh and cdf_thresh control which models are kept in the posterior sampling. Raising these thresholds reduces output file size at the cost of potentially trimming the tails of the posterior.

The cell below prints a reference table of these advanced parameters. It does not run a fit (which would be too slow for an interactive tutorial).

# Advanced BruteForce.fit() parameter reference
# This cell documents parameters -- it does NOT run a fit.

advanced_params = [
    # (parameter, default, description)
    ("logl_dim_prior", "True",
     "Use chi-square dimensionality prior on the likelihood.\n"
     "True  -> penalizes models based on chi2 relative to DOF (recommended).\n"
     "False -> standard Gaussian log-likelihood (no DOF correction)."),

    ("apply_agewt", "True",
     "Apply age-based weighting from MIST evolutionary tracks.\n"
     "Accounts for the different amounts of time a star spends\n"
     "in each evolutionary phase (e.g., main sequence vs. giant branch)."),

    ("apply_grad", "True",
     "Apply gradient-based correction for non-uniform grid spacing.\n"
     "Corrects for the fact that model grid points are not evenly\n"
     "distributed in parameter space (mass, age, metallicity)."),

    ("av_gauss", "None",
     "Gaussian prior on A(V) as a (mean, std) tuple.\n"
     "Example: av_gauss=(0.5, 0.2) centers extinction at A(V)=0.5 mag.\n"
     "If None, an essentially flat (0.0, 1e6) prior is used. Applied in\n"
     "ADDITION to any dust-map prior from dustfile (not exclusive), so an\n"
     "informative av_gauss plus dustfile counts the constraint twice."),

    ("wt_thresh", "1e-3",
     "Weight threshold for model selection: keeps models with\n"
     "weight > wt_thresh * max(weight). Lower values keep more\n"
     "models (broader posterior tails) but increase output size."),

    ("cdf_thresh", "2e-3",
     "CDF threshold for model selection. Models contributing less\n"
     "than this fraction of cumulative probability are discarded.\n"
     "Used as fallback if wt_thresh is None."),

    ("Ndraws", "250",
     "Number of posterior samples saved per object.\n"
     "More draws give smoother posteriors but larger output files.\n"
     "Typical range: 100 (quick) to 500 (publication quality)."),

    ("Nmc_prior", "50",
     "Number of Monte Carlo samples for prior integration.\n"
     "Higher values improve accuracy of the distance/extinction\n"
     "marginalization at the cost of computation time."),

    ("rv_gauss", "(3.32, 0.18)",
     "Gaussian prior on R(V) as a (mean, std) tuple.\n"
     "Default based on Schlafly et al. (2016).\n"
     "Set to None to use a flat prior within rvlim bounds."),

    ("running_io", "True",
     "If True, writes each object to disk immediately after fitting.\n"
     "Safer for long runs (recoverable on crash) but slower on\n"
     "network filesystems. Set False to batch-write at the end."),
]

# Print formatted table
print("Advanced BruteForce.fit() Parameters")
print("=" * 72)

for param, default, desc in advanced_params:
    print(f"\n  {param}  (default: {default})")
    print("  " + "-" * 68)
    for line in desc.split("\n"):
        print(f"    {line}")

print("\n" + "=" * 72)
print("\nExample usage with advanced options:")
print("""
    bf.fit(
        data=flux, data_err=flux_err, data_mask=mask,
        data_labels=obj_ids, save_file='results.h5',
        # --- advanced options ---
        logl_dim_prior=True,      # recommended for robust fitting
        apply_agewt=True,         # weight by evolutionary timescales
        apply_grad=True,          # correct for grid spacing
        av_gauss=(0.3, 0.1),     # informative A(V) prior
        rv_gauss=(3.32, 0.18),   # Schlafly+2016 R(V) prior
        wt_thresh=1e-3,           # keep top 0.1% of models
        cdf_thresh=2e-3,          # CDF cutoff
        Ndraws=250,               # posterior samples per source
        Nmc_prior=50,             # MC prior samples
        running_io=True,          # write results incrementally
    )
""")

Advanced BruteForce.fit() Parameters
========================================================================

  logl_dim_prior  (default: True)
  --------------------------------------------------------------------
    Use chi-square dimensionality prior on the likelihood.
    True  -> penalizes models based on chi2 relative to DOF (recommended).
    False -> standard Gaussian log-likelihood (no DOF correction).

  apply_agewt  (default: True)
  --------------------------------------------------------------------
    Apply age-based weighting from MIST evolutionary tracks.
    Accounts for the different amounts of time a star spends
    in each evolutionary phase (e.g., main sequence vs. giant branch).

  apply_grad  (default: True)
  --------------------------------------------------------------------
    Apply gradient-based correction for non-uniform grid spacing.
    Corrects for the fact that model grid points are not evenly
    distributed in parameter space (mass, age, metallicity).

  av_gauss  (default: None)
  --------------------------------------------------------------------
    Gaussian prior on A(V) as a (mean, std) tuple.
    Example: av_gauss=(0.5, 0.2) centers extinction at A(V)=0.5 mag.
    If None, an essentially flat (0.0, 1e6) prior is used. Applied in
    ADDITION to any dust-map prior from dustfile (not exclusive), so an
    informative av_gauss plus dustfile counts the constraint twice.

  wt_thresh  (default: 1e-3)
  --------------------------------------------------------------------
    Weight threshold for model selection: keeps models with
    weight > wt_thresh * max(weight). Lower values keep more
    models (broader posterior tails) but increase output size.

  cdf_thresh  (default: 2e-3)
  --------------------------------------------------------------------
    CDF threshold for model selection. Models contributing less
    than this fraction of cumulative probability are discarded.
    Used as fallback if wt_thresh is None.

  Ndraws  (default: 250)
  --------------------------------------------------------------------
    Number of posterior samples saved per object.
    More draws give smoother posteriors but larger output files.
    Typical range: 100 (quick) to 500 (publication quality).

  Nmc_prior  (default: 50)
  --------------------------------------------------------------------
    Number of Monte Carlo samples for prior integration.
    Higher values improve accuracy of the distance/extinction
    marginalization at the cost of computation time.

  rv_gauss  (default: (3.32, 0.18))
  --------------------------------------------------------------------
    Gaussian prior on R(V) as a (mean, std) tuple.
    Default based on Schlafly et al. (2016).
    Set to None to use a flat prior within rvlim bounds.

  running_io  (default: True)
  --------------------------------------------------------------------
    If True, writes each object to disk immediately after fitting.
    Safer for long runs (recoverable on crash) but slower on
    network filesystems. Set False to batch-write at the end.

========================================================================

Example usage with advanced options:

    bf.fit(
        data=flux, data_err=flux_err, data_mask=mask,
        data_labels=obj_ids, save_file='results.h5',
        # --- advanced options ---
        logl_dim_prior=True,      # recommended for robust fitting
        apply_agewt=True,         # weight by evolutionary timescales
        apply_grad=True,          # correct for grid spacing
        av_gauss=(0.3, 0.1),     # informative A(V) prior
        rv_gauss=(3.32, 0.18),   # Schlafly+2016 R(V) prior
        wt_thresh=1e-3,           # keep top 0.1% of models
        cdf_thresh=2e-3,          # CDF cutoff
        Ndraws=250,               # posterior samples per source
        Nmc_prior=50,             # MC prior samples
        running_io=True,          # write results incrementally
    )

7. Summary and Key Takeaways#

This tutorial has demonstrated how to fit individual stars using BruteForce:

Key Steps#

  1. Data Preparation

    • Convert magnitudes to flux (maggies)

    • Ensure proper error arrays and masks

    • Include parallax and coordinates when available

  2. Model Setup

    • Load pre-computed model grids (MIST)

    • Apply photometric offsets for calibration

    • Initialize BruteForce with models

  3. Running Fits

    • Configure priors (Galactic, dust, parallax)

    • Set sampling parameters (Ndraws, Nmc_prior)

    • Save results to HDF5 for analysis

  4. Result Analysis

    • Assess fit quality with goodness-of-fit p-values (dof = Nbands - 3)

    • Analyze parameter distributions and correlations

    • Visualize individual SEDs and posteriors

Best Practices#

  • Data Quality: Require SNR > 5 and >= 4 bands per source

  • Photometric Offsets: Always apply appropriate calibrations

  • Priors: Use 3D dust maps and parallax when available

  • Validation: Use goodness-of-fit p-values (p > 1e-3 indicates a good fit); chi2/Nbands is not a reliable quality metric across different band counts

  • Uncertainties: Use posterior samples for proper error estimation

Next Steps#

  • Tutorial 6: Cluster Analysis and Population Fitting

  • Tutorial 7: 3D Dust Mapping

  • Tutorial 8: Photometric Calibration

print("Tutorial 5 Complete!")
print("="*60)

tutorials_dir = Path(__file__).parent if '__file__' in dir() else Path('.')

print("\nOutput files:")
for h5_file in sorted(tutorials_dir.glob('Orion_l209.1_b-19.9_mist*.h5')):
    print(f"  - {h5_file.name}")

print("\nGenerated plots:")
for plot_file in sorted(plots_dir.glob('*.png')):
    print(f"  - {plot_file.name}")

Tutorial 5 Complete!
============================================================

Output files:
  - Orion_l209.1_b-19.9_mist.h5
  - Orion_l209.1_b-19.9_mist_nodust.h5

Generated plots:
  - cornerplot.png
  - dist_vs_red.png
  - dust_prior_comparison.png
  - fitting_results.png
  - posterior_predictive.png