BLiSS 1.0 usage example¶

from pathlib import Path
import sys

# Use the local package when running the notebook from the repository.
PROJECT_ROOT = Path.cwd().parent if Path.cwd().name == "notebooks" else Path.cwd()
sys.path.insert(0, str(PROJECT_ROOT))

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import bliss


from bliss import (
    load_text_spectrum,
    load_fits_spectrum,
    rebin_bins,
    rebin_snr,
    rebin_resolution,
    find_emission_lines,
    add_most_probable_ion,
    get_all_compatible_lines,
    plot_line_prob,
    create_bliss_results_folder,
)
from bliss.isis_interface import write_isis_line_model_files, run_bliss_for_isis

print("Using BLiSS from:", bliss.__file__)

results_dir = create_bliss_results_folder(PROJECT_ROOT / "results", "bliss")
results_dir

text_path = PROJECT_ROOT / "examples" / "xmm.dat"
energy, values, errors, bin_width = load_text_spectrum(text_path)

print(energy.shape, values.shape, errors.shape)
print("Energy range:", energy.min(), energy.max())

plt.figure(figsize=(8, 4))
plt.errorbar(energy, values, yerr=errors, fmt=".", alpha=0.7)
plt.xlabel("Energy (keV)")
plt.ylabel("Spectrum")
plt.title("Raw text spectrum")
plt.show()

rx_bins, ry_bins, rsy_bins = rebin_bins(energy, values, errors, 1)
rx_snr, ry_snr, rsy_snr = rebin_snr(energy, values, errors, 0.8)

plt.figure(figsize=(8, 4))
plt.errorbar(rx_snr, ry_snr, yerr=rsy_snr, fmt=".")
plt.xlabel("Energy (keV)")
plt.ylabel("Rebinned spectrum")
plt.title("Text spectrum rebinned by S/N")
plt.show()

xmm_lines = find_emission_lines(
    rx_snr,
    ry_snr,
    rsy_snr,
    show_plot=True,
    output_dir=results_dir,
    plot_name="xmm_bliss_fit.png",
)

xmm_lines.head()

xmm_lines.to_csv(results_dir / "xmm_candidate_lines.csv", index=False)
print("Saved:", results_dir / "xmm_candidate_lines.csv")
print("Number of candidates:", len(xmm_lines))

xmm_identified = add_most_probable_ion(xmm_lines, v_doppler_kms=1200)
xmm_identified.to_csv(results_dir / "xmm_identified_lines.csv", index=False)
xmm_identified.head()

all_xmm_identified = get_all_compatible_lines(xmm_lines, v_doppler_kms=200)
all_xmm_identified[0]

fig, ax = plot_line_prob(xmm_identified, show=False)
ax.set_title("BLiSS candidate probabilities and IDs")

fig.savefig(
    results_dir / "xmm_line_probability.png",
    dpi=150,
    bbox_inches="tight"
)

plt.show()

fits_dir = PROJECT_ROOT / "examples" / "fits_spectrum_data"
print(sorted(p.name for p in fits_dir.iterdir()))

plt.figure(figsize=(8, 4))
plt.errorbar(
    fits_spectrum.energy,
    fits_spectrum.values,
    yerr=fits_spectrum.uncertainties,
    fmt=".",
    alpha=0.6,
)
plt.xlabel("Energy (keV)")
plt.ylabel("Background-subtracted spectrum")
plt.title("Spectrum loaded from FITS files")
plt.show()

fx, fy, fsy = fits_spectrum.energy, fits_spectrum.values, fits_spectrum.uncertainties
valid = np.isfinite(fx) & np.isfinite(fy) & np.isfinite(fsy) & (fsy > 0)
fx, fy, fsy = fx[valid], fy[valid], fsy[valid]

# Optional: keep only positive values for the current emission-line workflow.
pos = fy > 0
fx, fy, fsy = fx[pos], fy[pos], fsy[pos]

# Use a modest rebinning for demonstration.
frx, fry, frsy = rebin_resolution(fx, fy, fsy, 0.02)

fits_lines = find_emission_lines(
    frx,
    fry,
    frsy,
    show_plot=True,
    output_dir=results_dir,
    plot_name="fits_bliss_fit.png",
)
fits_lines.to_csv(results_dir / "fits_candidate_lines.csv", index=False)
fits_lines.head()