Commit fabc232d authored by Andrea Giannetti's avatar Andrea Giannetti
Browse files

Added documentation; fixed bugs; updated tests.

parent 054bcbe6
Loading
Loading
Loading
Loading
+92 −21
Original line number Diff line number Diff line
@@ -3,7 +3,8 @@ import yaml
import numpy as np
import os
import urllib.request
from typing import List
from typing import (List,
                    Union)
from astropy import units as u
from assets.constants import (radmc_grid_map,
                              radmc_coord_system_map,
@@ -11,6 +12,11 @@ from assets.constants import (radmc_grid_map,


def load_config_file(config_file_path: str) -> dict:
    """
    Load the information in the YAML configuration file into a python dictionary
    :param config_file_path: path to the configuration file
    :return: a dictionary with the parsed information
    """
    with open(config_file_path) as config_file:
        config = yaml.load(config_file, Loader=yaml.FullLoader)
    return config
@@ -20,16 +26,39 @@ def compute_power_law_radial_profile(
        central_value: float,
        power_law_index: float,
        distance_matrix: np.array,
        value_at_reference: Union[float, None] = None,
        distance_reference: float = 1.0,
        fill_reference_pixel: bool = True
) -> np.array:
        fill_reference_pixel: bool = True) -> np.array:
    """
    Compute a power law distribution over the specified grid
    :param central_value: value in the center of the grid
    :param power_law_index: index of the power law used to scale the profile
    :param distance_matrix: the matrix of distances from the reference point
    :param value_at_reference: value of the profile at the reference distance; defaults to central value if not provided
    :param distance_reference: value of the reference distance
    :param fill_reference_pixel: whether to fill the reference point with central_value and set this to the maximum in
        the grid
    :return: the distance matrix
    """
    _value_at_reference = value_at_reference if value_at_reference is not None else central_value
    _distance_matrix = np.where(distance_matrix == 0, 1,
                                distance_matrix) if fill_reference_pixel is True else distance_matrix
    return central_value * (_distance_matrix / distance_reference) ** power_law_index
    profile = _value_at_reference * (_distance_matrix / distance_reference) ** power_law_index
    # If the routine fills the 0-distance point (the centre), it fixes the profile making the central value the maximum
    if fill_reference_pixel is True:
        profile = np.where(profile > central_value, central_value, profile)
    return profile


def get_grid_properties(grid_config: dict,
                        keyword: str) -> list:
    """
    Compute the grid axes properties using the grid configuration dictionary. If less than 3 axes are provided, fill
        the missing ones with the first axis.
    :param grid_config: the dictionary containing the grid configuration metadata
    :param keyword: the keyword to read
    :return: a list of the requested keyword values over the three spatial axes
    """
    grid_property_list = [grid_config['dim1'][keyword]]
    for axis_idx in range(2):
        try:
@@ -41,28 +70,38 @@ def get_grid_properties(grid_config: dict,

def compute_cartesian_coordinate_grid(indices: np.array,
                                      physical_px_size: np.array) -> np.array:
    """
    Compute the physical coordinate grid, for a regular grid
    :param indices: the indices of the grid pixels
    :param physical_px_size: the array of physical pixel sizes or the pixel size, as astropy.Quantities
    :return: the numpy.array with the physical grid coordinates
    """
    try:
        _physical_px_size_cm = [px_size.to(u.cm).value for px_size in physical_px_size]
    return indices * _physical_px_size_cm
    except TypeError:
        _physical_px_size_cm = physical_px_size.to(u.cm).value
    return (indices.T * _physical_px_size_cm).T


def extract_grid_metadata(config: dict) -> dict:
    """
    Enrich grid metadata from the information in the configuration file
    :param config: configuration dictionary
    :return: a dictionary with the metadata
    """
    grid_config = config['grid']

    grid_metadata = grid_config.copy()
    grid_metadata['grid_type'] = radmc_grid_map[grid_config['grid_type']]
    grid_metadata['coordinate_system'] = radmc_coord_system_map[grid_config['coordinate_system']]
    grid_metadata['grid_shape'] = get_grid_properties(grid_config=grid_config,
                                                      keyword='steps')
    grid_metadata['refpix'] = get_grid_properties(grid_config=grid_config,
                                                  keyword='refpix')
    grid_metadata['grid_size'] = get_grid_properties(grid_config=grid_config,
                                                     keyword='size')
    grid_metadata['grid_size_units'] = get_grid_properties(grid_config=grid_config,
                                                           keyword='unit')
    grid_properties_keywords = ['shape', 'refpix', 'size', 'size_units']
    for key in grid_properties_keywords:
        grid_metadata[f'grid_{key}'] = get_grid_properties(grid_config=grid_config,
                                                          keyword=key)

    grid_metadata['physical_px_size'] = get_physical_px_size(grid_metadata)

    indices = np.indices(grid_metadata['grid_shape']).T - grid_metadata['refpix']
    indices = get_centered_indices(grid_metadata)

    grid_metadata['coordinate_grid'] = compute_cartesian_coordinate_grid(indices=indices,
                                                                         physical_px_size=grid_metadata[
@@ -74,18 +113,40 @@ def extract_grid_metadata(config: dict) -> dict:
    return grid_metadata


def get_grid_edges(grid_metadata):
def get_centered_indices(grid_metadata: dict) -> np.array:
    """
    Recompute the indices using as reference the reference pixel in the configuration file
    :param grid_metadata: the dictionary with the grid metadata, to obtain the grid shape and the reference pixel
    :return: a numpy array of indices, wrt. the reference
    """
    return (np.indices(grid_metadata['grid_shape']).T - grid_metadata['grid_refpix']).T


def get_grid_edges(grid_metadata: dict) -> np.array:
    """
    Compute the physical coordinates at the grid boundaries
    :param grid_metadata: the dictionary with the grid metadata, to obtain the grid shape, the reference pixel, and the
        physical pixel size
    :return: a numpy array with the coordinates of the edges, per axis
    """
    grid_edges = []
    for axis_idx in range(len(grid_metadata['grid_shape'])):
        indices = np.arange(grid_metadata['grid_shape'][axis_idx] + 1) - grid_metadata['grid_refpix'][axis_idx] - 0.5
        # valid for regular grid
        grid_edges.append(
            (np.arange(grid_metadata['grid_shape'][axis_idx] + 1) - grid_metadata['refpix'][axis_idx] - 0.5) *
            grid_metadata['physical_px_size'][axis_idx])
        grid_edges.append(compute_cartesian_coordinate_grid(indices=indices,
                                                            physical_px_size=grid_metadata['physical_px_size'][axis_idx]))
    # Transposed for consistent flattening
    return np.array(grid_edges).T


def get_distance_matrix(grid_metadata, indices):
def get_distance_matrix(grid_metadata: dict,
                        indices: np.array) -> np.array:
    """
    Compute physical distance from the reference pixel, in cm
    :param grid_metadata: the dictionary with the grid metadata, to obtain the grid shape, and the physical pixel size
    :param indices: numpy array of the grid pixel indices
    :return: a numpy array with the euclidean distance from the reference pixel
    """
    distance_matrix = np.zeros(grid_metadata['grid_shape'])
    for axis_idx in range(len(grid_metadata['grid_shape'])):
        distance_matrix += (indices[axis_idx, ...] * grid_metadata['physical_px_size'][axis_idx].to(u.cm).value) ** 2
@@ -93,6 +154,12 @@ def get_distance_matrix(grid_metadata, indices):


def get_physical_px_size(grid_metadata: dict) -> List[u.Quantity]:
    """
    Compute the physical pixel size
    :param grid_metadata: the dictionary with the grid metadata, to obtain the grid shape, the grid size, and the grid
        size units
    :return: a list of physical pixel sizes, per axis
    """
    physical_px_size = []
    for axis_idx in range(len(grid_metadata['grid_shape'])):
        physical_px_size.append(
@@ -102,6 +169,11 @@ def get_physical_px_size(grid_metadata: dict) -> List[u.Quantity]:


def get_moldata(species_names: list):
    """
    Downloads the molecular data from the Leiden molecular database; check whether the molecule in mapped in
        leiden_url_mapping
    :param species_names: the names of the species for which to download data
    """
    for species in species_names:
        data = urllib.request.urlopen(leiden_url_mapping[species]).read().decode()
        with open(os.path.join('mdl', 'radmc_input_files', f'molecule_{species}.inp'), 'w') as outfile:
@@ -111,4 +183,3 @@ def get_moldata(species_names: list):
def convert_frequency_to_wavelength(frequency: astropy.units.Quantity,
                                    output_units: astropy.units.Unit):
    return frequency.to(output_units, equivalencies=u.spectral())
+1 −1
Original line number Diff line number Diff line
@@ -4,7 +4,7 @@ grid:
    central_density: 1e6
    density_unit: "cm^-3"
    density_profile: "constant"
    dim1: {"size":0.1, "unit": "pc", "steps": 5, "refpix": 2}
    dim1: {"size":0.1, "size_units": "pc", "shape": 5, "refpix": 2}
    source_distance: 1000
    source_distance_unit: "pc"

+64 −23
Original line number Diff line number Diff line
@@ -7,7 +7,7 @@ from assets.commons import (compute_power_law_radial_profile,
                            get_grid_edges,
                            extract_grid_metadata,
                            get_physical_px_size,
                            compute_cartesian_coordinate_grid, get_distance_matrix)
                            compute_cartesian_coordinate_grid, get_distance_matrix, get_centered_indices)
from unittest import TestCase


@@ -20,6 +20,9 @@ def create_test_config(config_dict: dict,


class Test(TestCase):
    def setUp(self):
        self.config_filename = 'config.yml'

    def test_compute_power_law_radial_profile_flat(self):
        indices = np.indices([3, 3]) - 1
        distance_matrix = np.sqrt(indices[0, :, :] ** 2 + indices[1, :, :] ** 2)
@@ -42,43 +45,43 @@ class Test(TestCase):
        self.assertTrue(np.array_equal(radial_profile, 1.0 / distance_matrix))

    def test_get_grid_properties_uniform(self):
        config_filename = os.path.join('test_files', 'config.yml')
        _config_filename = os.path.join('test_files', self.config_filename)
        config_dict = {
            'dim1': {"size": 0.1, "unit": "pc", "steps": 5, "refpix": 2},
            'dim1': {"size": 0.1, "unit": "pc", "shape": 5, "refpix": 2},
        }
        keywords = ['size', 'unit', 'steps', 'refpix']
        keywords = ['size', 'unit', 'shape', 'refpix']
        expected_result = {
            'size': [0.1] * 3,
            'unit': ['pc'] * 3,
            'steps': [5] * 3,
            'shape': [5] * 3,
            'refpix': [2] * 3,
        }
        create_test_config(config_dict=config_dict,
                           config_filename=config_filename)
        config = load_config_file(config_filename)
                           config_filename=_config_filename)
        config = load_config_file(_config_filename)
        for key in keywords:
            grid_properties = get_grid_properties(grid_config=config['grid'],
                                                  keyword=key)
            self.assertListEqual(grid_properties, expected_result[key])

    def test_get_grid_properties(self):
        config_filename = os.path.join('test_files', 'config.yml')
        _config_filename = os.path.join('test_files', self.config_filename)
        config_dict = {
            'dim1': {"size": 0.1, "unit": "pc", "steps": 5, "refpix": 2},
            'dim2': {"size": 0.2, "unit": "pc", "steps": 10, "refpix": 4.5},
            'dim3': {"size": 0.3, "unit": "pc", "steps": 15, "refpix": 7},
            'dim1': {"size": 0.1, "unit": "pc", "shape": 5, "refpix": 2},
            'dim2': {"size": 0.2, "unit": "pc", "shape": 10, "refpix": 4.5},
            'dim3': {"size": 0.3, "unit": "pc", "shape": 15, "refpix": 7},
        }
        keywords = ['size', 'unit', 'steps', 'refpix']
        keywords = ['size', 'unit', 'shape', 'refpix']
        expected_result = {
            'size': [0.1, 0.2, 0.3],
            'unit': ['pc'] * 3,
            'steps': [5, 10, 15],
            'shape': [5, 10, 15],
            'refpix': [2, 4.5, 7],
        }
        create_test_config(config_dict=config_dict,
                           config_filename=config_filename)
                           config_filename=_config_filename)

        config = load_config_file(config_filename)
        config = load_config_file(_config_filename)
        for key in keywords:
            grid_properties = get_grid_properties(grid_config=config['grid'],
                                                  keyword=key)
@@ -88,7 +91,7 @@ class Test(TestCase):
        grid_metadata = {
            'grid_shape': [2, 2],
            'physical_px_size': [1 * u.cm, 1 * u.cm],
            'refpix': [0.5, 0.5]
            'grid_refpix': [0.5, 0.5]
        }
        expected_results_per_axis = np.array([-1, 0, 1])
        grid_edges = get_grid_edges(grid_metadata=grid_metadata)
@@ -96,7 +99,7 @@ class Test(TestCase):
            self.assertTrue(np.array_equal(expected_results_per_axis, grid_edges[:, axis_idx]))

    def test_get_physical_px_size(self):
        config_filename = os.path.join('test_files', 'config.yml')
        _config_filename = os.path.join('test_files', self.config_filename)
        grid_size = 0.2
        npix = 5
        config_dict = {
@@ -104,14 +107,14 @@ class Test(TestCase):
            'coordinate_system': 'cartesian',
            'central_density': '1e6',
            'density_unit': 'cm^-3',
            'dim1': {"size": grid_size, "unit": "pc", "steps": npix, "refpix": 2},
            'dim1': {"size": grid_size, "size_units": "pc", "shape": npix, "refpix": 2},
            'source_distance': '1000',
            'source_distance_unit': 'pc'
        }
        expected_result = [(grid_size * u.pc).to(u.cm) / npix] * 3
        create_test_config(config_dict=config_dict,
                           config_filename=config_filename)
        config = load_config_file(config_filename)
                           config_filename=_config_filename)
        config = load_config_file(_config_filename)
        grid_metadata = extract_grid_metadata(config=config)
        physical_px_size = get_physical_px_size(grid_metadata)
        self.assertListEqual(expected_result, physical_px_size)
@@ -124,7 +127,15 @@ class Test(TestCase):
                                                          physical_px_size=physical_px_size)
        self.assertTrue(np.array_equal(expected_results, computed_grid))

    def test_get_distance_matrix(self):
    def test_compute_cartesian_coordinate_grid(self):
        indices = np.indices([3,4])
        physical_px_size = [1 * u.cm, 1 * u.cm]
        expected_results = indices
        computed_grid = compute_cartesian_coordinate_grid(indices=indices,
                                                          physical_px_size=physical_px_size)
        self.assertTrue(np.array_equal(expected_results, computed_grid))

    def test_get_distance_matrix_2d(self):
        grid_metadata = {
            'grid_shape': [2, 2],
            'physical_px_size': [1 * u.cm, 1 * u.cm]
@@ -134,3 +145,33 @@ class Test(TestCase):
        distance_matrix = get_distance_matrix(grid_metadata=grid_metadata,
                                              indices=indices)
        self.assertTrue(np.array_equal(expected_results, distance_matrix))

    def test_get_distance_matrix(self):
        grid_metadata = {
            'grid_shape': [2, 2, 2],
            'physical_px_size': [1 * u.cm, 1 * u.cm, 1 * u.cm]
        }
        indices = np.indices(grid_metadata['grid_shape'])
        expected_results = np.array([[[0, 1], [1, np.sqrt(2)]],
                                     [[1, np.sqrt(2)], [np.sqrt(2), np.sqrt(3)]]])
        distance_matrix = get_distance_matrix(grid_metadata=grid_metadata,
                                              indices=indices)
        self.assertTrue(np.array_equal(expected_results, distance_matrix))

    def test_get_centered_indices(self):
        grid_metadata = {
            'grid_shape': [5, 3],
            'grid_refpix': [2, 1]
        }
        expected_indices = [[[-2, -2, -2],
                             [-1, -1, -1],
                             [0, 0, 0],
                             [1, 1, 1],
                             [2, 2, 2]],
                            [[-1, 0, 1],
                             [-1, 0, 1],
                             [-1, 0, 1],
                             [-1, 0, 1],
                             [-1, 0, 1]]]
        indices = get_centered_indices(grid_metadata=grid_metadata)
        self.assertTrue(np.array_equal(np.array(expected_indices), indices))