Moved some functions to commons and created stub for aperture spectra fitting

import pyspeckit

import logging

import sys

import yaml

import numpy as np

import os

import pickle

import scipy

from astropy import units as u

from astropy.constants import c as speed_of_light

from typing import (List,

                    Union,

                    Tuple)

def setup_logger(name: str,

                 log_level: str = None) -> logging.Logger:

    """

    Configure default logger

    :param name: Logger name

    :param log_level: Logging levels, as defined by logging

    :return: the logger object

    """

    log_level = 'DEBUG' if log_level is None else log_level

    """general logger configurator"""

    logger = logging.getLogger(name)

    logger.setLevel(log_level)

    handler = logging.StreamHandler(sys.stdout)

    handler.setLevel(log_level)

    formatter = logging.Formatter("%(asctime)s - %(name)-7s - %(levelname)s - %(message)s", datefmt="%Y-%m-%d %H:%M:%S")

    handler.setFormatter(formatter)

    logger.addHandler(handler)

    return logger

def validate_parameter(param_to_validate,

                       default):

    """

    Substitute default value if the parameter is set to None

    :param param_to_validate: input to validate

    :param default: default value

    :return: default value if parameter is set to None, else the passed value

    """

    return param_to_validate if param_to_validate is not None else default

def load_config_file(config_file_path: str,

                     override_config: Union[dict, None] = None) -> dict:

    """

    Load the information in the YAML configuration file into a python dictionary

    :param config_file_path: path to the configuration file

    :param override_config: parameters of the input file to override (e.g. for grid creation)

    :return: a dictionary with the parsed information

    """

    _override_config = validate_parameter(override_config, default={})

    with open(config_file_path) as config_file:

        config = yaml.load(config_file, Loader=yaml.FullLoader)

    for key in _override_config:

        if key not in config.keys():

            config[key] = {}

        for subkey in _override_config[key]:

            try:

                config[key][subkey] = _override_config[key][subkey]

            except KeyError:

                config[key] = {subkey: _override_config[key][subkey]}

    return config

def convert_ratio_to_density(ratio):

    return 10 ** ((np.arctanh((ratio - 0.523) / 1.009) + 6.086) ** 0.971)

def gradient_plane_3d(xy, ra_grad, dec_grad, offset):

    return ra_grad * xy[0] - dec_grad * xy[1] + offset

def get_velocities_and_linewidths_from_fit(velocity_indices, linewidth_indices, species, source):

    with open(f'final_results_{source}_{species}.pickle', 'rb') as infile:

        fit_results = pickle.load(infile)

    return (fit_results[velocity_indices[species], :, :],

            fit_results[linewidth_indices[species], :, :])

def convert_dvelocity_to_hz(dvelocity, rest_frequencies):

    return dvelocity / speed_of_light.to(u.km / u.s) * rest_frequencies[0].to(u.Hz)

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

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

    if smooth_data is True:

        smoothed_data = np.ones_like(cube.cube)

        for i in range(cube.cube.shape[0]):

            smoothed_data[i, :, :] = scipy.ndimage.gaussian_filter(cube.cube[i, :, :], sigma=3)

        cube.cube = smoothed_data

    if frequency_limits is not None:

        cube = cube.slice(frequency_limits[0], frequency_limits[1], unit=u.GHz)

    return cube

def fill_guesses(line_peak, doppler_adjusted_frequencies, shift_hz, linewidth_guess_hz, rms,

                 peak_scaling_factors=(1, 1, 0.5, 0.3)):

    guesses = []

    for idx, freq in enumerate(doppler_adjusted_frequencies):

        guesses += [np.max((rms, peak_scaling_factors[idx] * line_peak)),

                    freq.to(u.Hz).value + shift_hz,

                    linewidth_guess_hz]

    return guesses

species: 'N2H+'

#source_names: ['G332.96', 'G28', 'G14']

source_names: ['G14', ]

max_components: 3

rest_frequencies: {

    'CH3OH': [96.739358, 96.74142, 96.744545, 96.755501],

import pyspeckit

import yaml

import os

import matplotlib

import matplotlib.pyplot as plt

import numpy as np

from astropy import units as u

from assets.commons import get_cube, convert_dvelocity_to_hz, fill_guesses

from assets.constants import source_properties, config_mdl, datadir

def fill_tied_pars(tied_pars):

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

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

    return _tied_pars_list

source = 'G28'

species = 'CH3OH'

vlsr = source_properties[source]['vlsr'] * u.km / u.s

rms = source_properties[source][species]['rms']

linewidth_guess = source_properties[source][species]['linewidth_guess_km_s'] * u.km / u.s

frequency_limits = source_properties[source][species]['frequency_limits_ghz']

rest_frequencies = np.array(config_mdl['rest_frequencies'][species]) * u.GHz

freq_to_vel = u.doppler_radio(rest_frequencies[0])

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

doppler_adjusted_frequencies = rest_frequencies + frequency_offset

freq_to_vel = u.doppler_radio(rest_frequencies[0])

linewidth_guess_hz = convert_dvelocity_to_hz(dvelocity=linewidth_guess,

                                             rest_frequencies=rest_frequencies).value / 2.355

delta_freq = rest_frequencies[1:] - rest_frequencies[0]

pars_per_component = 3 * len(rest_frequencies)

tied = []

guesses = []

peak_guess = 0.03

manual_components = 1

for components in range(manual_components):

    idx_shift = components * pars_per_component

    if len(rest_frequencies) > 1:

        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

        }

        tied += fill_tied_pars(tied_pars=tied_pars)

        guesses += fill_guesses(line_peak=peak_guess,

                                doppler_adjusted_frequencies=doppler_adjusted_frequencies,

                                shift_hz=linewidth_guess_hz / 2 * idx_shift,

                                linewidth_guess_hz=linewidth_guess_hz,

                                rms=rms)

matplotlib.use('TkAgg')

source = 'G28'

species = 'CH3OH'

filename = source_properties[source][species]['data']

cube = get_cube(filename=filename,

                datadir=datadir,

                frequency_limits=frequency_limits,

                smooth_data=config_mdl['smooth_data'])

regions = {

    'centre': (256, 256, 20)

}

plotdir = os.path.join('prs', 'figures', 'aperture_spectra_fits')

for region in regions:

    apspec = cube.get_apspec(regions[region])

    apspec.plotter()

    apspec.specfit(fittype='gaussian',

                   guesses=guesses,

                   tied=tied,

                   clear=True,)

    apspec.plotter.savefig(os.path.join(plotdir, f'aperture_spectra_{source}_{species}_{region}.svg'))

import matplotlib.pyplot as plt

import matplotlib

import astropy.units as u

import pyspeckit

import pickle

import yaml

import scipy

from astropy.constants import c as speed_of_light

from scipy.signal import peak_widths

from etl.assets.commons import validate_parameter

from assets.commons import validate_parameter, get_cube, convert_dvelocity_to_hz, fill_guesses

from assets.constants import source_names

from typing import Union, Tuple

matplotlib.use('TkAgg')

    return bic

def convert_dvelocity_to_hz(dvelocity, rest_frequencies):

    return dvelocity / speed_of_light.to(u.km / u.s) * rest_frequencies[0].to(u.Hz)

def create_limits(par_number, channelwidth_hz, rms, max_linewidth: Union[float, None] = None):

    limits = []

    _max_linewidth = validate_parameter(max_linewidth, 6e6)

    return _tied_pars_list

def fill_guesses(line_peak, doppler_adjusted_frequencies, shift_hz, linewidth_guess_hz, rms,

                 peak_scaling_factors=(1, 1, 0.5, 0.3)):

    guesses = []

    for idx, freq in enumerate(doppler_adjusted_frequencies):

        guesses += [np.max((rms, peak_scaling_factors[idx] * line_peak)),

                    freq.to(u.Hz).value + shift_hz,

                    linewidth_guess_hz]

    return guesses

def update_results(fitted_cube,

                   final_pars,

                   indices_to_keep,

rest_frequencies = np.array(config['rest_frequencies'][config['species']]) * u.GHz

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

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

    if smooth_data is True:

        smoothed_data = np.ones_like(cube.cube)

        for i in range(cube.cube.shape[0]):

            smoothed_data[i, :, :] = scipy.ndimage.gaussian_filter(cube.cube[i, :, :], sigma=3)

        cube.cube = smoothed_data

    if frequency_limits is not None:

        cube = cube.slice(frequency_limits[0], frequency_limits[1], unit=u.GHz)

    return cube

def get_mask_and_starting_point(source_properties, data, rms, threshold_sigma=4):

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

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

    if 'starting_point' in source_properties:

        starting_point = source_properties['starting_point']

    else:

        # spec_sum = -np.where(spec_sum < 0, 0, spec_sum)

        # spec_sum[230:240, 295:305] *= -1

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

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

    return mask, starting_point, spec_sum

for source in config['source_names']:

for source in source_names:

    print(f'Processing {source}')

    filename = source_properties[source][config['species']]['data']

    vlsr = source_properties[source]['vlsr'] * u.km / u.s

    cube = get_cube(filename=filename,

                    datadir=datadir,

                    frequency_limits=frequency_limits,

                    smooth_data=config['smooth_data'])

    freq_to_vel = u.doppler_radio(rest_frequencies[0])

                         limits=limits,

                         limited=limited,

                         maskmap=mask,

                         # use_nearest_as_guess=True,

                         use_neighbor_as_guess=True,

                         multicore=10)

            # cube.plot_spectrum(x=starting_point[1], y=starting_point[0])

            # cube.plot_fit(x=starting_point[1], y=starting_point[0])

            # matplotlib.use('Qt5Agg')

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

            starting_point = np.unravel_index(np.argmax(np.nansum(residuals, axis=0)), spec_sum.shape)

            peaks = scipy.signal.find_peaks(residuals[:, starting_point[0], starting_point[1]], height=3 * rms, width=3)

    with open(f'final_results_{source}_{config["species"]}.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()

import pyspeckit

import yaml

import os

import matplotlib

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

with open(os.path.join(datadir, 'source_properties.yml'), 'r') as infile:

    source_properties = yaml.load(infile, Loader=yaml.FullLoader)

with open(os.path.join('mdl', 'config', 'config.yml'), 'r') as infile:

    config_mdl = yaml.load(infile, Loader=yaml.FullLoader)

with open(os.path.join('config', 'config_visfit.yml'), 'r') as infile:

    config = yaml.load(infile, Loader=yaml.FullLoader)

matplotlib.use('TkAgg')

source = 'G28'

species = 'CH3OH'

filename = source_properties[source][species]['data']

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

apspec = cube.get_apspec((256,256,20))