Commit 303b5de3 authored by Andrea Giannetti's avatar Andrea Giannetti
Browse files

Added first version of the script to get clump average density

parent d3e8292a

etl/assets/commons/__init__.py

+28 −0
Changes for etl/assets/commons/__init__.py: 28 added lines, 0 removed lines.
Original line number Diff line number Diff line
@@ -66,6 +66,34 @@ def load_config_file(config_file_path: str, 
    return config





def filter_dict_by_key(d, key_root):

    return {k: v for k, v in d.items() if k.startswith(key_root)}





def find_zeros(x, y):

    """

    Finds the zeros of a discretely sampled function using linear interpolation.



    Parameters:

    x (array-like): The x-values of the sampled function.

    y (array-like): The y-values of the sampled function.



    Returns:

    np.ndarray: The estimated x-values where the function crosses zero.

    """

    x = np.asarray(x)

    y = np.asarray(y)

    zeros = []



    for i in range(len(y) - 1):

        if y[i] * y[i + 1] < 0:  # Sign change indicates a zero crossing

            # Linear interpolation formula: x_zero = x1 - y1 * (x2 - x1) / (y2 - y1)

            x_zero = x[i] - y[i] * (x[i + 1] - x[i]) / (y[i + 1] - y[i])

            zeros.append(float(x_zero))



    return zeros





def convert_ratio_to_density(ratio_value, ratio_name=None):

    """

    Convert the line ratio of methanol to a guess for the number density of H2

etl/assets/templates/tied_3mm.template

+2 −2
Changes for etl/assets/templates/tied_3mm.template: 2 added lines, 2 removed lines.
Original line number Diff line number Diff line 
,,,

,p[{PAR_VELO}]+{DELTA_FREQ_1},p[{PAR_SIGMA_LINE}],

,,p[{PAR_SIGMA_LINE}],

,p[{PAR_VELO}]+{DELTA_FREQ_1},,

,p[{PAR_VELO}]+{DELTA_FREQ_2},p[{PAR_SIGMA_LINE}],

,p[{PAR_VELO}]+{DELTA_FREQ_3},p[{PAR_SIGMA_LINE}]
 
 No newline at end of file

etl/config/ch3oh_regions.yml

+1 −1
Changes for etl/config/ch3oh_regions.yml: 1 added line, 1 removed line.
Original line number Diff line number Diff line
@@ -14,7 +14,7 @@ G28: { 
  'envelope_05': [317, 259, 25, 1],

  'envelope_06': [200, 194, 25, 1],

  'envelope_07': [167, 336, 25, 3],

  'envelope_core_01': [150, 200, 10, 1],

  'core_envelope_01': [150, 200, 10, 1],

}

G14: {

    'core_cont_01': [220, 260, 8, 1],

etl/config/guesses.yml

+1 −1
Changes for etl/config/guesses.yml: 1 added line, 1 removed line.
Original line number Diff line number Diff line
@@ -74,7 +74,7 @@ G28: { 
    [ 6.0e+5, -6.0e+5, -2.0e+5 ],

    [ 2.0e+5, 2.0e+5, 2.0e+5 ]

  ],

  envelope_core_01: [

  core_envelope_01: [

    [ [ 0.025, 0.035, 0.003, 0.001 ], ],

    [ 0, ],

    [ 2.5e+5, ]

etl/mdl/mdl_extract_ch3oh_spectra.py

+44 −63
Changes for etl/mdl/mdl_extract_ch3oh_spectra.py: 44 added lines, 63 removed lines.
Original line number Diff line number Diff line
@@ -5,8 +5,9 @@ import matplotlib.pyplot as plt 
import numpy as np

import pymc as pm

import arviz as az

import pandas as pd

from sklearn.neighbors import KernelDensity

from scipy.integrate import cumulative_trapezoid

from astropy.coordinates import SkyCoord

from astropy import units as u

from astropy.io import fits

from astropy.wcs import WCS
@@ -17,11 +18,15 @@ from assets.commons import (setup_logger, 
                            convert_dfreq_to_velocity,

                            fill_guesses,

                            validate_parameter,

                            convert_ratio_to_density)

                            convert_ratio_to_density,

                            filter_dict_by_key,

                            find_zeros)

from matplotlib.collections import PatchCollection

from assets.constants import source_properties, config_mdl, datadir, continuum_fits

from typing import Union



plt.rcParams['font.size'] = 11



logger = setup_logger(name='AP. SPECTRA', log_level='DEBUG')
@@ -136,34 +141,6 @@ def create_limits(par_number, channelwidth_hz, rms, max_linewidth: Union[float, 
    return limits, limited





def filter_dict_by_key(d, key_root):

    return {k: v for k, v in d.items() if k.startswith(key_root)}





def find_zeros(x, y):

    """

    Finds the zeros of a discretely sampled function using linear interpolation.



    Parameters:

    x (array-like): The x-values of the sampled function.

    y (array-like): The y-values of the sampled function.



    Returns:

    np.ndarray: The estimated x-values where the function crosses zero.

    """

    x = np.asarray(x)

    y = np.asarray(y)

    zeros = []



    for i in range(len(y) - 1):

        if y[i] * y[i + 1] < 0:  # Sign change indicates a zero crossing

            # Linear interpolation formula: x_zero = x1 - y1 * (x2 - x1) / (y2 - y1)

            x_zero = x[i] - y[i] * (x[i + 1] - x[i]) / (y[i + 1] - y[i])

            zeros.append(float(x_zero))



    return zeros





def compute_density_pdf(density, density_axis=None):

    _density_axis = validate_parameter(density_axis, np.linspace(2.5, 7.5, 500))

    kde = KernelDensity(bandwidth=0.1, kernel='epanechnikov')
@@ -229,10 +206,11 @@ def extract_fit_parameters(parinfo, idx_shift: int) -> dict: 
        idx_shift (int): Index offset corresponding to the velocity component.



    Returns:

        dict: Contains WIDTH and AMPLITUDE values for transitions 86, 87, and 88.

        dict: Contains WIDTH, SHIFT, and AMPLITUDE values for transitions 86, 87, and 88.

    """

    try:

        fwhm_hz = parinfo[f'WIDTH{idx_shift}'].value

        frequency = parinfo[f'SHIFT{idx_shift}'].value

        amp_86 = parinfo[f'AMPLITUDE{idx_shift}'].value

        amp_87 = parinfo[f'AMPLITUDE{idx_shift + 2}'].value

        amp_88 = parinfo[f'AMPLITUDE{idx_shift + 3}'].value
@@ -241,6 +219,7 @@ def extract_fit_parameters(parinfo, idx_shift: int) -> dict: 


    return {

        'fwhm_hz': fwhm_hz,

        'frequency': frequency,

        'amplitudes': {

            '86': amp_86,

            '87': amp_87,
@@ -434,11 +413,11 @@ def prepare_fit_guesses_and_ties(source, region, component_count, 




def run_and_plot_specfit(apspec, tied, line_guesses, limits, limited,

                         source, region, species, plotdir):

                         source, region, species, plotdir, name_prefix=None):

    """

    Execute the Gaussian specfit, plot the result, and return the residual RMS.

    """

    fig = plt.figure()

    _name_prefix = validate_parameter(name_prefix, default="aperture_spectra")

    apspec.specfit(

        fittype="gaussian",

        guesses=line_guesses,
@@ -456,13 +435,12 @@ def run_and_plot_specfit(apspec, tied, line_guesses, limits, limited, 
    # plot fit components

    apspec.plotter(title=f"{source} {region}")

    apspec.specfit.plot_fit(show_components=True, annotate=False)

    fig.savefig(os.path.join(

        plotdir, f"aperture_spectra_{source}_{species}_{region}.svg"

    ))

    plt.close(fig)

    plt.savefig(os.path.join(

        plotdir, f"{_name_prefix}_{source}_{species}_{region}.svg"

    ), bbox_inches='tight')



    # return the RMS of the residuals

    return float(np.nanstd(apspec.specfit.residuals)), fig

    return float(np.nanstd(apspec.specfit.residuals))





def overlay_region_on_map(apspec,
@@ -504,12 +482,13 @@ def overlay_region_on_map(apspec, 
    ax_map.set_xlabel("RA [J2000]")

    ax_map.set_ylabel("Dec [J2000]")

    fig_map.tight_layout()

    fig_map.savefig(plotfile)

    fig_map.savefig(plotfile, bbox_inches='tight')

    plt.close(fig_map)



def process_source(source, species, guesses):

    # 1. Load maps and prepare line models

    moment_map, wcs = load_moment_map(source, species)

    freq_to_vel = u.doppler_radio(config_mdl["rest_frequencies"][species][0] * u.GHz)

    freqs_dop, lw_guess_hz, max_lw_hz, rms, freq_limits = prepare_line_model_params(

        source, species, config_mdl, source_properties

    )
@@ -529,7 +508,7 @@ def process_source(source, species, guesses): 
    # 4. Loop over regions

    for region, props in regions[source].items():

        coords, components = props[:3], props[3]

        densities[region] = {'best_fit': [], 'uncertainty': [], 'reliability_score': [], 'fwhm': []}

        densities[region] = {'best_fit': [], 'uncertainty': [], 'reliability_score': [], 'fwhm': [], 'vlsr': []}

        fwhms[region], reliability_score[region] = [], {}



        # 4a. Build & run the spec‐fit
@@ -538,7 +517,7 @@ def process_source(source, species, guesses): 
            source, region, components, guesses, delta_freq, rms, pars_per_component, freqs_dop,

            channelwidth_hz, coords, max_lw_hz

        )

        spec_rms, _ = run_and_plot_specfit(

        spec_rms = run_and_plot_specfit(

            apspec=apspec,

            tied=tied,

            line_guesses=line_guesses,
@@ -555,15 +534,16 @@ def process_source(source, species, guesses): 
            idx = vel_comp * 4



            # extract fit params & convert to physical units

            fitp = extract_fit_parameters(apspec.specfit.parinfo, idx)

            fwhm = convert_dfreq_to_velocity(fitp['fwhm_hz'], freqs_dop)

            fitpars = extract_fit_parameters(apspec.specfit.parinfo, idx)

            fwhm = convert_dfreq_to_velocity(fitpars['fwhm_hz'], freqs_dop)

            velocity = (fitpars['frequency']* u.Hz).to(u.km / u.s, equivalencies=freq_to_vel).value

            if should_skip_component(source, region, vel_comp, fwhm):

                continue



            # compute SNR, ratios, reliability

            rms_int = compute_rms_integral(fitp['fwhm_hz'], channelwidth_hz, spec_rms)

            ints   = compute_integrated_intensities(fitp['amplitudes'], fitp['fwhm_hz'])

            snr_p  = compute_snr_params(integrated_intensities=ints, amplitudes=fitp['amplitudes'], spec_rms=spec_rms, rms_integral=rms_int)

            rms_int = compute_rms_integral(fitpars['fwhm_hz'], channelwidth_hz, spec_rms)

            ints   = compute_integrated_intensities(fitpars['amplitudes'], fitpars['fwhm_hz'])

            snr_p  = compute_snr_params(integrated_intensities=ints, amplitudes=fitpars['amplitudes'], spec_rms=spec_rms, rms_integral=rms_int)



            r87, s87 = compute_ratio_posteriors(ints['86'], ints['87'], rms_int, 0.05)

            r88, s88 = compute_ratio_posteriors(ints['86'], ints['88'], rms_int, 0.001)
@@ -596,6 +576,7 @@ def process_source(source, species, guesses): 
            densities[region]['uncertainty'].append(hpd)

            densities[region]['reliability_score'].append(float(score))

            densities[region]['fwhm'].append(float(fwhm.value))

            densities[region]['vlsr'].append(float(velocity))

            fwhms[region].append(fwhm.value)



        # 4c. Overlay region on map
@@ -611,9 +592,7 @@ def process_source(source, species, guesses): 
                            regions=regions,

                            densities=densities,

                            moment_map=moment_map,

                            wcs=wcs,

                            continuum_data=continuum_data,

                            wcs_cont=wcs_cont)

                            wcs=wcs)

    with open(os.path.join('mdl', 'inferred_densities', f'densities_{source}.yml'), 'w') as outfile:

        yaml.dump(densities, outfile, default_flow_style=False)
@@ -626,47 +605,49 @@ def plot_density_vs_fwhm(densities, fwhms, source, species, regions): 
            color = 'blue' if region.startswith('core') else 'cyan'

            plt.scatter(densities[region]['best_fit'], fwhms[region], label=region, edgecolors=color, figure=fig)

    plt.legend()

    plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_density_FWHM_{source}_{species}.svg'))

    plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_density_FWHM_{source}_{species}.svg'), bbox_inches='tight')

    plt.close(fig)





def overlay_regions_density(source, regions, densities, moment_map, wcs, continuum_data, wcs_cont):

def overlay_regions_density(source, regions, densities, moment_map, wcs):

    region_types = ['envelope', 'core']



    for rtype in region_types:

    plt.clf()

    ax = plt.subplot(projection=wcs)

    plt.imshow(np.log10(moment_map.data), origin='lower')

        plt.colorbar(label=r'log(K km s$^{-1}$)')

    plt.colorbar(label=r'log(K km s$^{-1}$)', pad=0.10, shrink=0.68)

    ax.set_autoscale_on(False)

    cores_data = pd.read_csv(os.path.join(datadir, 'ashes_cores.csv'))

    core_data_filtered = cores_data[(cores_data['Clump'] == source) & (cores_data['Class'] > 0)]



    for rtype in region_types:

        linestyle = '-' if rtype == 'core' else '--'

        regions_to_plot = filter_dict_by_key(regions[source], rtype)

        patches = []

        density_values = []

        for reg in regions_to_plot:

            patch = Circle(regions_to_plot[reg][0:2], regions_to_plot[reg][2], edgecolor='black', facecolor='none',

                           zorder=10)

                           zorder=1)

            patches.append(patch)

            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 = PatchCollection(patches, cmap=matplotlib.cm.hot, zorder=2, linestyles=linestyle)

        pc.set_array(density_values)

        pc.set_facecolor('none')

        pc.set_clim(4.8, 5.8)

        ax.add_collection(pc)

        # plt.contour(continuum_data,

        #             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')

    coord = SkyCoord(core_data_filtered['RA'], core_data_filtered['DEC'], frame='icrs')

    ax.scatter(coord.ra.deg, coord.dec.deg, transform=ax.get_transform('world'), s=100, marker='*',

               edgecolor='lightcyan', facecolor='none', zorder=1, linewidths=1)

    plt.ylabel('Dec')

    plt.xlabel('RA')

        plt.colorbar(pc, label=r'log($n($H$_2)$ [cm$^{-3}$])')

    plt.colorbar(pc, label=r'log($n($H$_2)$ [cm$^{-3}$])', shrink=0.68)

    plt.tight_layout()

        plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_{source}_{rtype}_{species}.svg'))

    plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_{source}_{species}.svg'), bbox_inches='tight')

    plt.close()



sources = ['G14', 'G28', 'G332.96']

# sources = ['G28']

# sources = ['G332.96']

species = 'CH3OH'

samples = 1000

prob_threshold = 0.67