First commit

,,,

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

,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

import os

import warnings

import numpy as np

import matplotlib.pyplot as plt

import matplotlib

import astropy.units as u

import pyspeckit

import pickle

from astropy.constants import c as speed_of_light

matplotlib.use('Qt5Agg')

datadir = os.path.join('assets', 'data')

filename = 'G332.96_CH3OH_3mm_feather.fits'

# filename = 'test_subcube.fits'

cube = pyspeckit.Cube(os.path.join(datadir, filename), format='fits')

rest_frequencies = np.array([96.739358, 96.74142, 96.744545, 96.755501]) * u.GHz

freq_to_vel = u.doppler_radio(rest_frequencies[0])

vlsr = -67.6 * u.km / u.s

frequency_offset = vlsr.to(u.Hz, equivalencies=freq_to_vel) - rest_frequencies[0]

doppler_adjusted_frequencies = rest_frequencies + frequency_offset

cube = cube.slice(96.755, 96.775, unit=u.GHz)

delta_freq = rest_frequencies.diff()

rms = 0.0034

# Test double component

# cube.cube[:, 0, 0] = np.roll(cube.cube[:, 0, 0], -15) + cube.cube[:, 0, 0]

spec_sum = np.nansum(cube.cube, axis=0)

starting_point = np.unravel_index(np.argmax(spec_sum), spec_sum.shape)

threshold_spec_sum = 3 * np.sqrt(cube.shape[0]) * rms

fit_mask = np.where(spec_sum > threshold_spec_sum, True, False)

tied_pars = {

    'PAR_VELO': 1,

    'DELTA_FREQ_1': delta_freq[0].to(u.Hz).value,

    'DELTA_FREQ_2': delta_freq[1].to(u.Hz).value,

    'DELTA_FREQ_3': delta_freq[2].to(u.Hz).value,

    'PAR_SIGMA_LINE': 2

}

with open(os.path.join('assets', 'templates', 'tied_3mm.template'), 'r') as template:

    tied = template.read().replace('\n', '').format(**tied_pars).split(',')

limits = [(rms, 0) for idx in range(len(tied))]

limited = [(True, False) for idx in range(len(tied))]

guesses = [0.029, doppler_adjusted_frequencies[0].to(u.Hz).value, 50000,

           0.029, doppler_adjusted_frequencies[1].to(u.Hz).value, 50000,

           0.01, doppler_adjusted_frequencies[2].to(u.Hz).value, 50000,

           0.005, doppler_adjusted_frequencies[3].to(u.Hz).value, 50000]

cube.fiteach(guesses=guesses,

             start_from_point=(starting_point[1], starting_point[0]),

             tied=tied,

             limits=limits,

             limited=limited,

             maskmap=fit_mask,

             use_nearest_as_guess=True)

# cube.plot_spectrum(6, 6)

# cube.plot_fit(6, 6)

residuals = cube.cube - cube.get_modelcube(update=True)

number_of_parameters = tied.count('')

bic = - np.nansum(-residuals ** 2, axis=0) / (rms ** 2) - cube.shape[0] * np.log(np.sqrt(2 * np.pi) * rms) + np.log(

    cube.shape[0]) * number_of_parameters

threshold = 2 * (cube.shape[0]) - cube.shape[0] * np.log(np.sqrt(2 * np.pi) * rms) + np.log(

    cube.shape[0]) * number_of_parameters

fit_results = cube.parcube.copy()

for idx_frequency in range(number_of_parameters // 6):

    idx_shift = idx_frequency * 12

    velocity_map = (cube.parcube[1 + idx_shift, :, :] * u.Hz).to(u.km / u.s, equivalencies=freq_to_vel).value

    linewidth_map = 2.355 * ((cube.parcube[2 + idx_shift, :, :] * u.Hz) / rest_frequencies[0].to(u.Hz)).value * speed_of_light.to(

        u.km / u.s).value

    ratio_map_87_86 = (cube.parcube[6 + idx_shift, :, :] / cube.parcube[idx_shift, :, :])

    fit_results[1 + idx_shift, :, :] = velocity_map

    fit_results[2 + idx_shift, :, :] = linewidth_map

# plt.clf()

# plt.imshow(velocity_map.value)

# plt.colorbar()

# plt.show()

# plt.clf()

# plt.imshow(linewidth_map.value)

# plt.colorbar()

# plt.show()

# plt.clf()

# plt.imshow(ratio_map_87_86)

# plt.colorbar()

# plt.show()

max_components = 3

final_pars = np.zeros((number_of_parameters * max_components, cube.parcube.shape[1], cube.parcube.shape[2])) * np.nan

indices_to_keep = np.array((0, 1, 2, 3, 6, 9))

_indices_to_keep = np.array((0, 1, 2, 3, 6, 9))

final_pars[0:number_of_parameters, :, :] = np.where(fit_mask, fit_results[indices_to_keep, :, :], np.nan)

for components in range(max_components - 1):

    refit_mask = np.where(bic > threshold, True, False)

    if np.sum(refit_mask) > 0:

        idx_shift = (components + 1) * 12

        tied_pars = {

            'PAR_VELO': 1 + idx_shift,

            'DELTA_FREQ_1': delta_freq[0].to(u.Hz).value,

            'DELTA_FREQ_2': delta_freq[1].to(u.Hz).value,

            'DELTA_FREQ_3': delta_freq[2].to(u.Hz).value,

            'PAR_SIGMA_LINE': 2 + idx_shift

        }

        with open(os.path.join('assets', 'templates', 'tied_3mm.template'), 'r') as template:

            tied += template.read().replace('\n', '').format(**tied_pars).split(',')

        guesses += [0.029, doppler_adjusted_frequencies[0].to(u.Hz).value + 1e5 * idx_shift, 50000,

                   0.029, doppler_adjusted_frequencies[1].to(u.Hz).value + 1e5 * idx_shift, 50000,

                   0.01, doppler_adjusted_frequencies[2].to(u.Hz).value + 1e5 * idx_shift, 50000,

                   0.005, doppler_adjusted_frequencies[3].to(u.Hz).value + 1e5 * idx_shift, 50000]

        limits = [(rms, 0) for idx in range(len(tied))]

        limited = [(True, False) for idx in range(len(tied))]

        starting_point = np.unravel_index(np.argmax(spec_sum[refit_mask]), spec_sum.shape)

        cube.fiteach(guesses=guesses,

                     start_from_point=(starting_point[1], starting_point[0]),

                     tied=tied,

                     limits=limits,

                     limited=limited,

                     maskmap=refit_mask,

                     use_nearest_as_guess=True)

        residuals = cube.cube - cube.get_modelcube(update=True)

        number_of_parameters = tied.count('')

        new_bic = - np.nansum(-residuals ** 2, axis=0) / (rms ** 2) - cube.shape[0] * np.log(

            np.sqrt(2 * np.pi) * rms) + np.log(cube.shape[0]) * number_of_parameters

        fit_results = cube.parcube.copy()

        for idx_to_change in range(components + 2):

            velocity_map = (cube.parcube[1 + (12 * idx_to_change), :, :] * u.Hz).to(u.km / u.s, equivalencies=freq_to_vel).value

            linewidth_map = 2.355 * ((cube.parcube[2 + (12 * idx_to_change), :, :] * u.Hz) / rest_frequencies[0].to(u.Hz)).value * speed_of_light.to(

                    u.km / u.s).value

            ratio_map_87_86 = (cube.parcube[6 + (12 * idx_to_change), :, :] / cube.parcube[(12 * idx_to_change), :, :])

            fit_results[1 + (12 * idx_to_change), :, :] = velocity_map

            fit_results[2 + (12 * idx_to_change), :, :] = linewidth_map

        _indices_to_keep = np.concatenate((_indices_to_keep, (indices_to_keep + idx_shift)))

        final_pars[0:6 + 6 * (components + 1), :, :] = np.where((bic / new_bic > 1) & (refit_mask),

                                                          fit_results[

                                                          _indices_to_keep, :, :],

                                                          final_pars[0:6 + 6 * (components + 1), :, :])

        print(bic/new_bic)

        bic = new_bic

with open('final_results.pickle', 'wb') as handle:

    pickle.dump(final_pars, handle, protocol=pickle.HIGHEST_PROTOCOL)

# spec = cube.get_spectrum(5,5)

# spec.specfit(guesses = guesses, tied = tied, annotate = False)

# spec.plotter()

# spec.specfit.plot_fit()

# spec.plotter.refresh()

numpy

matplotlib

astropy

PyYAML==6.0

scipy

astroquery

pandas

pyspeckit

spectral-cube

 No newline at end of file