Ratios computed through PyMC to avoid artificial peaks in pdfs (b341772e) · Commits · Andrea Giannetti / CoCoA_methanol_archive

etl/assets/commons/init.py

+1 −1

Changes for etl/assets/commons/__init__.py: 1 added line, 1 removed line.
Original line number	Diff line number	Diff line
		@@ -111,7 +111,7 @@ def convert_dvelocity_to_hz(dvelocity, rest_frequencies):


		def convert_dfreq_to_velocity(dfrequency, rest_frequencies):
		return dfrequency / rest_frequencies[0].to(u.Hz) * speed_of_light.to(u.km / u.s)
		return dfrequency * u.Hz / rest_frequencies[0].to(u.Hz) * speed_of_light.to(u.km / u.s)


		def get_cube(filename, datadir, frequency_limits=None, smooth_data=False):

etl/mdl/mdl_extract_ch3oh_spectra.py

+165 −41

Changes for etl/mdl/mdl_extract_ch3oh_spectra.py: 165 added lines, 41 removed lines.
Original line number	Diff line number	Diff line
		@@ -3,13 +3,16 @@ import os
		import matplotlib
		import matplotlib.pyplot as plt
		import numpy as np
		import pymc as pm
		import arviz as az
		from sklearn.neighbors import KernelDensity
		from scipy.integrate import cumulative_trapezoid
		from astropy import units as u
		from astropy.io import fits
		from astropy.wcs import WCS
		from matplotlib.patches import Circle
		from assets.commons import (get_cube,
		from assets.commons import (setup_logger,
		get_cube,
		convert_dvelocity_to_hz,
		convert_dfreq_to_velocity,
		fill_guesses,
		@@ -19,6 +22,107 @@ from matplotlib.collections import PatchCollection
		from assets.constants import source_properties, config_mdl, datadir, continuum_fits
		from typing import Union

		logger = setup_logger(name='AP. SPECTRA', log_level='DEBUG')


		def build_ratio_model(hatI1, hatI2, sigma, lower_limit_ratio):
		"""
		Build a PyMC model for the ratio R = I2 / I1, where:
		- hatI1: observed intensity of line 1
		- sigma1: noise std of line 1
		- L: detection limit (upper limit) for line 2
		- sigma2: noise std of line 2
		Priors for I1 and I2 are LogNormal(mu, sigma_prior).
		"""
		model = pm.Model()
		with model:
		# latent ratio R with a lower bound to forbid R→0
		R = pm.TruncatedNormal("R", mu=hatI2/hatI1, sigma=0.3, lower=lower_limit_ratio, upper=0.95)

		# data container for I1
		I1 = pm.Data("I1", hatI1)

		# predicted I2 = R * I1
		mu2 = R * I1

		# detection mask
		is_det = hatI2 >= 3 * sigma

		# likelihood if detected
		logp_det = pm.Normal.logp(hatI2, mu=mu2, sigma=sigma)

		# censored likelihood if not detected
		z = (3 * sigma - mu2) / sigma
		cdf_z = 0.5 * (1 + pm.math.erf(z / pm.math.sqrt(2)))
		logp_cens = pm.math.log(cdf_z + 1e-15)

		# combine via switch
		total_logp = pm.math.switch(is_det, logp_det, logp_cens)
		pm.Potential("obs_R", total_logp)

		return model


		def sample_ratio_model(model, draws=2000, tune=1000, target_accept=0.9, random_seed=42):
		"""
		Sample from the given PyMC model and return the trace.
		"""
		with model:
		trace = pm.sample(draws=draws, tune=tune, target_accept=target_accept, random_seed=random_seed)
		return trace


		def summarize_ratio_posterior(R_samples):
		"""
		Summarize the ratio posterior with mean, median, and 95% HPD.
		"""
		mean = np.mean(R_samples)
		median = np.median(R_samples)
		hpd = az.hdi(R_samples, hdi_prob=0.95)
		return {"mean": mean, "median": median, "95%_HPD": hpd}


		def compute_reliability_score(component):
		score = 0
		thresholds = {
		'SNR_intensity1': 50,
		'SNR_intensity2': 15,
		'SNR_intensity3': 10,
		'SNR1': 25,
		'SNR2': 10,
		'SNR3': 5,
		'ratio1': [0.1, 0.8],
		'ratio2': [0.1, 0.75],
		'rel_prec_ratio1': 0.7,
		'rel_prec_ratio2': 0.5,
		}
		weights = {
		'SNR_intensity1': 1,
		'SNR_intensity2': 3,
		'SNR_intensity3': 5,
		'SNR1': 1,
		'SNR2': 5,
		'SNR3': 10,
		'ratio1': 1,
		'ratio2': 1,
		'rel_prec_ratio1': 3,
		'rel_prec_ratio2': 3,
		}
		only_thresh = ['ratio1', 'ratio2']
		continuous_thresholds = [thresholds[key] for key in thresholds.keys() if key not in only_thresh]
		continuous_weights = [weights[key] for key in weights.keys() if key not in only_thresh]
		missing_weights = [weights[key] for key in weights.keys() if key in only_thresh]
		max_score = np.sum(np.array(continuous_thresholds) * np.array(continuous_weights)) + np.sum(missing_weights)

		for par in thresholds.keys():
		if type(thresholds[par]) is not list:
		score += weights[par] * np.min((component[par], thresholds[par]))
		else:
		if thresholds[par][0] < component[par] < thresholds[par][1]:
		score += weights[par]

		return score / max_score


		def fill_tied_pars(tied_pars):
		with open(os.path.join('assets', 'templates', 'tied_3mm.template'), 'r') as template:
		@@ -69,8 +173,10 @@ def find_zeros(x, y):


		sources = ['G14', 'G28', 'G332.96']
		sources = ['G28']
		species = 'CH3OH'
		samples = 100
		samples = 1000
		prob_threshold = 0.67
		matplotlib.use('TkAgg')

		plotdir = os.path.join('prs', 'figures', 'aperture_spectra_fits')
		@@ -131,6 +237,7 @@ for source in sources:

		densities = {}
		fwhms = {}
		reliability_score = {}
		for region in regions[source]:
		print(region)
		_region_coordinates = regions[source][region][0:3]
		@@ -140,8 +247,9 @@ for source in sources:
		line_guesses = []
		tied = []

		densities[region] = {'best_fit': [], 'uncertainty': []}
		densities[region] = {'best_fit': [], 'uncertainty': [], 'reliability_score': []}
		fwhms[region] = []
		reliability_score[region] = {}
		for component in range(components):
		idx_shift = component * pars_per_component

		@@ -175,35 +283,61 @@ for source in sources:
		fig = plt.figure(figsize=(15, 5))
		ax_s = fig.add_subplot(122)
		apspec.plotter(title=f'{source} {region}', figure=fig)
		# plt.annotate('a', xy=(0.8, 0.75), xycoords='figure fraction')
		apspec.specfit.plot_fit(show_components=True, annotate=False)
		apspec.plotter.savefig(os.path.join(plotdir, f'aperture_spectra_{source}_{species}_{region}.svg'))

		spec_rms = np.nanstd(apspec.specfit.residuals)

		for vel_component in range(components):
		skipped = False
		if ((source == 'G332.96') or ((source == 'G14') and (region == 'envelope_01'))) and (vel_component > 0):
		pass
		skipped = True
		else:
		idx_shift = vel_component * 4
		line_1 = np.random.default_rng().normal(scale=spec_rms, size=samples) + apspec.specfit.parinfo[
		f'AMPLITUDE{idx_shift}'].value
		line_2 = np.random.default_rng().normal(scale=spec_rms, size=samples) + apspec.specfit.parinfo[
		f'AMPLITUDE{idx_shift + 2}'].value
		line_3 = np.random.default_rng().normal(scale=spec_rms, size=samples) + apspec.specfit.parinfo[
		f'AMPLITUDE{idx_shift + 3}'].value
		ratio_87 = np.where(line_2 <= 0, 0.3 * spec_rms, line_2) / np.where(line_1 <= 0, spec_rms, line_1)
		ratio_88 = np.where(line_3 <= 0, 0.3 * spec_rms, line_3) / np.where(line_1 <= 0, spec_rms, line_1)
		fwhm = convert_dfreq_to_velocity(dfrequency=apspec.specfit.parinfo[f'WIDTH{idx_shift}'].value,
		fwhm_hz = apspec.specfit.parinfo[f'WIDTH{idx_shift}'].value
		fwhm = convert_dfreq_to_velocity(dfrequency=fwhm_hz,
		rest_frequencies=rest_frequencies)
		# if apspec.specfit.parinfo[f'AMPLITUDE{idx_shift+3}'].value > 1.5 * spec_rms:
		# ratio2 = apspec.specfit.parinfo[f'AMPLITUDE{idx_shift+3}'].value / apspec.specfit.parinfo[f'AMPLITUDE{idx_shift}'].value
		# else:
		# ratio2 = np.nan


		nchans = fwhm_hz / np.abs(channelwidth_hz)
		rms_integral = np.sqrt(2 * nchans) * spec_rms * np.abs(channelwidth_hz)
		reliability_score[region][vel_component] = {'pars': {}, 'score': 0}
		if (fwhm < 5 * u.km / u.second):
		int_intensity_86 = 1.0645 * apspec.specfit.parinfo[f'AMPLITUDE{idx_shift}'].value * fwhm_hz
		int_intensity_87 = 1.0645 * apspec.specfit.parinfo[f'AMPLITUDE{idx_shift + 2}'].value * fwhm_hz
		int_intensity_88 = 1.0645 * apspec.specfit.parinfo[f'AMPLITUDE{idx_shift + 3}'].value * fwhm_hz
		reliability_score[region][vel_component]['pars']['SNR_intensity1'] = int_intensity_86 / rms_integral
		reliability_score[region][vel_component]['pars']['SNR_intensity2'] = int_intensity_87 / rms_integral
		reliability_score[region][vel_component]['pars']['SNR_intensity3'] = int_intensity_88 / rms_integral
		reliability_score[region][vel_component]['pars']['SNR1'] = apspec.specfit.parinfo[
		f'AMPLITUDE{idx_shift}'].value / spec_rms
		reliability_score[region][vel_component]['pars']['SNR2'] = apspec.specfit.parinfo[
		f'AMPLITUDE{idx_shift + 2}'].value / spec_rms
		reliability_score[region][vel_component]['pars']['SNR3'] = apspec.specfit.parinfo[
		f'AMPLITUDE{idx_shift + 3}'].value / spec_rms
		model_8786 = build_ratio_model(hatI1=int_intensity_86,
		hatI2=int_intensity_87,
		sigma=rms_integral,
		lower_limit_ratio=0.05)
		sampled_ratio_8786 = sample_ratio_model(model=model_8786)
		ratio_posterior_8786 = sampled_ratio_8786.posterior["R"].values.flatten()
		ratio_posterior_summary_8786 = summarize_ratio_posterior(R_samples=ratio_posterior_8786)

		model_8886 = build_ratio_model(hatI1=int_intensity_86,
		hatI2=int_intensity_88,
		sigma=rms_integral,
		lower_limit_ratio=0.001)
		sampled_ratio_8886 = sample_ratio_model(model=model_8886)
		ratio_posterior_8886 = sampled_ratio_8886.posterior["R"].values.flatten()
		ratio_posterior_summary_8886 = summarize_ratio_posterior(R_samples=ratio_posterior_8886)
		reliability_score[region][vel_component]['pars']['ratio1'] = np.mean(ratio_posterior_8786)
		reliability_score[region][vel_component]['pars']['ratio2'] = np.mean(ratio_posterior_8886)
		reliability_score[region][vel_component]['pars']['rel_prec_ratio1'] = 1 - np.std(
		ratio_posterior_8786) / np.mean(ratio_posterior_8786)
		reliability_score[region][vel_component]['pars']['rel_prec_ratio2'] = 1 - np.std(
		ratio_posterior_8886) / np.mean(ratio_posterior_8886)
		reliability_score[region][vel_component]['score'] = compute_reliability_score(
		component=reliability_score[region][vel_component]['pars'])
		density_pdfs = {}
		for ratio_string, ratio in zip(['87-86', '88-86'], [ratio_87, ratio_88]):
		for ratio_string, ratio in zip(['87-86', '88-86'], [ratio_posterior_8786, ratio_posterior_8886]):
		density = convert_ratio_to_density(ratio, ratio_name=ratio_string)
		if (np.isnan(density).sum() == len(density)) or (np.nanmean(ratio) > 0.95):
		pdf = np.ones(500)
		@@ -211,14 +345,19 @@ for source in sources:
		else:
		density_axis, pdf = compute_density_pdf(density)
		density_pdfs[ratio_string] = {'density': density_axis, 'pdf': pdf}
		final_density_pdf = density_pdfs['87-86']['pdf'] * density_pdfs['88-86']['pdf']
		mix_weight = 0.3
		final_density_pdf = density_pdfs['87-86']['pdf'] * (
		mix_weight * density_pdfs['88-86']['pdf'] + (1 - mix_weight))
		norm = np.trapezoid(final_density_pdf, density_axis)
		final_density_pdf /= norm
		probability_mass = cumulative_trapezoid(final_density_pdf, density_axis)
		if fwhm.value < 5:
		hpd = find_zeros(density_axis, final_density_pdf - 0.67)
		if (fwhm.value < 5) and (skipped is False):
		logger.debug(f'Reliability parameters: {reliability_score[region][vel_component]["pars"]}')
		logger.debug(f'Reliability score: {reliability_score[region][vel_component]["score"]}')
		hpd = find_zeros(density_axis, final_density_pdf - prob_threshold)
		densities[region]['best_fit'].append(float(density_axis[np.argmax(final_density_pdf)]))
		densities[region]['uncertainty'].append(hpd)
		densities[region]['reliability_score'].append(float(reliability_score[region][vel_component]['score']))
		fwhms[region].append(fwhm.value)

		# plt.clf()
		@@ -227,21 +366,6 @@ for source in sources:
		# plt.plot(density_axis, density_pdfs['88-86']['pdf'])
		# plt.show()

		# fig = plt.figure()
		# density_axis, pdf_87 = compute_density_pdf(density)
		# _, pdf_88 = compute_density_pdf(density)
		# kde2 = KernelDensity(bandwidth=0.1, kernel='epanechnikov')
		# kde2.fit(np.log10(density2).reshape(-1, 1))
		# z2 = np.exp(kde2.score_samples(list(density_axis.reshape(-1, 1))))
		# plt.plot(density_axis, z * z2)
		# sns.kdeplot(np.log10(density), figure=fig)
		# sns.kdeplot(np.log10(density2), figure=fig)
		# plt.show()
		# print(np.nanmean(ratio), np.nanmean(ratio2), np.nanmean(density), np.nanstd(density),
		# np.nanmean(density2), np.nanstd(density2), spec_rms)

		# plt.show()

		# plt.clf()
		# hdu_cont = fits.open(os.path.join(datadir, continuum_fits[source]))[0]
		# continuum_data = np.where(hdu_cont.data >= 2.5e-4, hdu_cont.data, np.nan)
		@@ -264,7 +388,6 @@ for source in sources:
		plt.tight_layout()
		plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_{source}_{region}_{species}.svg'))


		fig = plt.figure(figsize=(7, 5))
		plt.clf()
		for region in regions[source]:
		@@ -293,7 +416,8 @@ for source in sources:
		patch = Circle(regions_to_plot[reg][0:2], regions_to_plot[reg][2], edgecolor='black', facecolor='none',
		zorder=10)
		patches.append(patch)
		density_values.append(np.max(densities[reg]['best_fit']))
		most_reliable_density_idx = np.argmax(densities[reg]['reliability_score'])
		density_values.append(densities[reg]['best_fit'][most_reliable_density_idx])
		pc = PatchCollection(patches, cmap=matplotlib.cm.hot, zorder=10)
		pc.set_array(density_values)
		pc.set_facecolor('none')

etl/prs/prs_mom0_dust.py

+10 −7

Changes for etl/prs/prs_mom0_dust.py: 10 added lines, 7 removed lines.
Original line number	Diff line number	Diff line
		@@ -77,14 +77,16 @@ for species in ['N2H+', 'CH3OH']:
		hdu_cont = fits.open(os.path.join(datadir, continuum_fits[source]))[0]
		continuum_data = np.where(hdu_cont.data >= 2.5e-4, hdu_cont.data, np.nan)

		wcs_cont = WCS(hdu_cont.header)
		plt.clf()
		wcs = WCS(hdu_cont.header)
		plt.subplot(projection=wcs)
		wcs = WCS(moment_zero_header)
		ax = plt.subplot(projection=wcs)
		plt.imshow(np.log10(moment_zero.value), origin='lower')
		plt.colorbar(label=r'log(K km s$^{-1}$)')
		plt.contour(continuum_data,
		levels=np.logspace(np.log10(2.5e-4), np.log10(0.9 * np.nanmax(continuum_data)), num=3),
		colors='red',
		levels=np.logspace(np.log10(2.5e-4), np.log10(0.9 * np.nanmax(continuum_data)), num=5),
		transform=ax.get_transform(wcs_cont),
		colors='grey',
		origin='lower')
		plt.ylabel('Dec[J2000]')
		plt.xlabel('RA[J2000]')
		@@ -92,13 +94,14 @@ for species in ['N2H+', 'CH3OH']:
		plt.savefig(os.path.join('prs', 'figures', f'{source}_{species}_mom0_continuum.svg'))

		plt.clf()
		wcs = WCS(hdu_cont.header)
		wcs = WCS(moment_zero_header)
		plt.subplot(projection=wcs)
		plt.imshow(np.log10(moment_zero.value) - np.log10(continuum_data), origin='lower')
		plt.colorbar()
		plt.contour(continuum_data,
		levels=np.logspace(np.log10(2.5e-4), np.log10(0.9 * np.nanmax(continuum_data)), num=3),
		colors='red',
		levels=np.logspace(np.log10(2.5e-4), np.log10(0.9 * np.nanmax(continuum_data)), num=5),
		transform=ax.get_transform(wcs_cont),
		colors='grey',
		origin='lower')
		plt.ylabel('Dec[J2000]')
		plt.xlabel('RA[J2000]')

etl/requirements.txt

+1 −0

Changes for etl/requirements.txt: 1 added line, 0 removed lines.
Original line number	Diff line number	Diff line
		@@ -3,6 +3,7 @@ astroquery==0.4.7
		matplotlib==3.9.2
		numpy==2.1.2
		pandas==2.2.3
		pymc==5.22.0
		pyspeckit==1.0.3
		PyYAML==6.0
		scipy==1.14.1