From 85df712d53ad424ec78e5c2c7e834ae81848ba8d Mon Sep 17 00:00:00 2001 From: Arjun Savel Date: Thu, 2 Jul 2020 08:29:55 -0700 Subject: [PATCH 01/21] Clean up error reporting formal integral --- tardis/montecarlo/formal_integral.py | 133 ++++++++++++++++----------- 1 file changed, 80 insertions(+), 53 deletions(-) diff --git a/tardis/montecarlo/formal_integral.py b/tardis/montecarlo/formal_integral.py index 872e0410e2c..2788f2309c0 100644 --- a/tardis/montecarlo/formal_integral.py +++ b/tardis/montecarlo/formal_integral.py @@ -313,7 +313,8 @@ def _formal_integral(self, iT, inu, inu_size, att_S_ul, Jred_lu, Jblue_lu, N): # find first contributing lines nu_start = nu * z[0] nu_end = nu * z[1] - self.line_search(nu_start, size_line, idx_nu_start) + idx_nu_start = line_search(self.runner.line_list_nu, + nu_start, size_line, idx_nu_start) offset = shell_id[0] * size_line # start tracking accumulated e-scattering optical depth @@ -456,63 +457,89 @@ def calculate_z(self, r, p, inv_t): else: return 0 - @njit(**njit_dict) - def line_search(self, nu, nu_insert, number_of_lines, result): - """ - Insert a value in to an array of line frequencies +class BoundsError(ValueError): + pass - Inputs: - :nu: (array) line frequencies - :nu_insert: (int) value of nu key - :number_of_lines: (int) number of lines in the line list - - Outputs: - index of the next line ot the red. - If the key value is redder - than the reddest line returns number_of_lines. - """ - # TODO: fix the TARDIS_ERROR_OK - # tardis_error_t ret_val = TARDIS_ERROR_OK # check - imin = 0 - imax = number_of_lines - 1 - if nu_insert > nu[imin]: - result = imin - elif nu_insert < nu[imax]: - result = imax + 1 +@njit(**njit_dict) +def line_search(nu, nu_insert, number_of_lines, result): + """ + Insert a value in to an array of line frequencies + + Inputs: + :nu: (array) line frequencies + :nu_insert: (int) value of nu key + :number_of_lines: (int) number of lines in the line list + + Outputs: + index of the next line ot the red. + If the key value is redder + than the reddest line returns number_of_lines. + """ + # TODO: fix the TARDIS_ERROR_OK + # tardis_error_t ret_val = TARDIS_ERROR_OK # check + imin = 0 + imax = number_of_lines - 1 + if nu_insert > nu[imin]: + result = imin + elif nu_insert < nu[imax]: + result = imax + 1 + else: + result = reverse_binary_search(nu, nu_insert, imin, imax, result) + result = result + 1 + return result + +@njit(**njit_dict) +def reverse_binary_search(x, x_insert, imin, imax, result): + """Look for a place to insert a value in an inversely sorted float array. + + Inputs: + :x: (array) an inversely (largest to lowest) sorted float array + :x_insert: (value) a value to insert + :imin: (int) lower bound + :imax: (int) upper bound + + Outputs: + index of the next boundary to the left + """ + # ret_val = TARDIS_ERROR_OK # check + if x_insert > x[imin] or x_insert < x[imax]: + raise BoundsError # check + else: + imid = (imin + imax) >> 1 + while imax - imin > 2: + if (x[imid] < x_insert): + imax = imid + 1 + else: + imin = imid + imid = (imin + imax) >> 1 + if (imax - imin == 2 and x_insert < x[imin + 1]): + result = imin + 1 else: - ret_val = self.reverse_binary_search(nu, nu_insert, imin, imax, result) - result = result + 1 - return ret_val + result = imin + return result - @njit(**njit_dict) - def reverse_binary_search(self, x, x_insert, imin, imax, result): - """Look for a place to insert a value in an inversely sorted float array. +@njit(**njit_dict) +def binary_search(x, x_insert, imin, imax, result): + # TODO: actually return result + if x_insert < x[imin] or x_insert > x[imax]: + raise BoundsError + else: + while imax >= imin: + imid = (imin + imax) / 2 + if x[imid] == x_insert: + result = imid + break + elif x[imid] < x_insert: + imin = imid + 1 + else: + imax = imid - 1 + if imax - imid == 2 and x_insert < x[imin + 1]: + result = imin + else: + result = imin # check + return result - Inputs: - :x: (array) an inversely (largest to lowest) sorted float array - :x_insert: (value) a value to insert - :imin: (int) lower bound - :imax: (int) upper bound - Outputs: - index of the next boundary to the left - """ - # ret_val = TARDIS_ERROR_OK # check - if x_insert > x[imin] or x_insert < x[imax]: - ret_val = TARDIS_ERROR_BOUNDS_ERROR # check - else: - imid = (imin + imax) >> 1 - while imax - imin > 2: - if (x[imid] < x_insert): - imax = imid + 1 - else: - imin = imid - imid = (imin + imax) >> 1 - if (imax - imin == 2 and x_insert < x[imin + 1]): - result = imin + 1 - else: - result = imin - return ret_val @njit(**njit_dict) def trapezoid_integration(array, h, N): From 6f43c4feba1d2f1fa3dfa6ac64067e6dd3acb334 Mon Sep 17 00:00:00 2001 From: Arjun Savel Date: Thu, 2 Jul 2020 08:33:17 -0700 Subject: [PATCH 02/21] Make first move to numba montecarlo tests --- tardis/montecarlo/tests/test_cmontecarlo.py | 1058 ------------------ tardis/montecarlo/tests/test_montecarlo.py | 1064 +++++++++++++++++++ 2 files changed, 1064 insertions(+), 1058 deletions(-) delete mode 100644 tardis/montecarlo/tests/test_cmontecarlo.py diff --git a/tardis/montecarlo/tests/test_cmontecarlo.py b/tardis/montecarlo/tests/test_cmontecarlo.py deleted file mode 100644 index ff7519e251a..00000000000 --- a/tardis/montecarlo/tests/test_cmontecarlo.py +++ /dev/null @@ -1,1058 +0,0 @@ -""" -Unit tests for methods in `tardis/montecarlo/src/cmontecarlo.c`. -* `ctypes` library is used to wrap C methods and expose them to python. - - -Probable Reasons for Failing Tests: ------------------------------------ - -1. Change made in C struct declarations: - - Reflect the changes done in C structs, into Python counterparts. - - Check **tardis/montecarlo/struct.py**. - -2. Return type of any method changed: - - Modify the `restype` parameter in the test method here. - - For example: - ``` - clib.rpacket_doppler_factor.restype = c_double - ``` - -3. Underlying logic modified: - - Check whether the changes made in C method are logically correct. - - If the changes made were correct and necessary, update the corresponding - test case. - - -General Test Design Procedure: ------------------------------- - -Please follow this design procedure while adding a new test: - -1. Parametrization as per desire of code coverage. - - C tests have different flows controlled by conditional statements. - Parameters checked in conditions can be provided in different testcases. - - Keep consistency with variable names as (in order): - - `packet_params` - - `model_params` - - `expected_params` (`expected` if only one value to be asserted.) - - Suggested variable names can be compromised if readability of the test - increases. - -2. Test Method body: - - Keep name as `test_` + `(name of C method)`. - - Refer to method `test_rpacket_doppler_factor` below for description. -""" - - -import os -import pytest -import numpy as np -import pandas as pd - - -from ctypes import ( - CDLL, - byref, - c_uint, - c_int64, - c_double, - c_ulong, - c_void_p, - cast, - POINTER, - pointer, - CFUNCTYPE - ) - -from numpy.testing import ( - assert_equal, - assert_almost_equal, - assert_array_equal, - assert_allclose - ) - -from tardis import __path__ as path -from tardis.montecarlo.struct import ( - RPacket, StorageModel, RKState, - PhotoXsect1level, - TARDIS_ERROR_OK, - TARDIS_ERROR_BOUNDS_ERROR, - TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE, - TARDIS_PACKET_STATUS_IN_PROCESS, - TARDIS_PACKET_STATUS_EMITTED, - TARDIS_PACKET_STATUS_REABSORBED, - CONTINUUM_OFF, - CONTINUUM_ON, - INVERSE_C, - BoundFreeTreatment -) - - -@pytest.fixture(scope='module') -def continuum_compare_data_fname(): - fname = 'continuum_compare_data.hdf' - return os.path.join(path[0], 'montecarlo', 'tests', 'data', fname) - - -@pytest.fixture(scope='module') -def continuum_compare_data(continuum_compare_data_fname, request): - compare_data = pd.HDFStore(continuum_compare_data_fname, mode='r') - - def fin(): - compare_data.close() - request.addfinalizer(fin) - - return compare_data - - -@pytest.fixture(scope="function") -def expected_ff_emissivity(continuum_compare_data): - emissivities = continuum_compare_data['ff_emissivity'] - - def ff_emissivity(t_electron): - emissivity = emissivities[t_electron] - nu_bins = emissivity['nu_bins'].values - emissivity_value = emissivity['emissivity'].dropna().values - - return nu_bins, emissivity_value - - return ff_emissivity - - -@pytest.fixture(scope='module') -def get_rkstate(continuum_compare_data): - data = continuum_compare_data['z2rkstate_key'] - pos_data = continuum_compare_data['z2rkstate_pos'] - - def z2rkstate(z_random): - key = (c_ulong * 624)(*data.loc[z_random].values) - pos = pos_data.loc[z_random] - return RKState( - key=key, - pos=pos, - has_gauss=0, - gauss=0.0 - ) - - return z2rkstate - - -@pytest.fixture(scope='function') -def model_w_edges(ion_edges, model): - photo_xsect = (POINTER(PhotoXsect1level) * len(ion_edges))() - - for i, edge in enumerate(ion_edges): - x_sect_1level = PhotoXsect1level() - for key, value in edge.items(): - if key in ['nu', 'x_sect']: - value = (c_double * len(value))(*value) - setattr(x_sect_1level, key, value) - photo_xsect[i] = pointer(x_sect_1level) - - no_of_edges = len(ion_edges) - continuum_list_nu = (c_double * no_of_edges)(*[edge['nu'][0] for edge in ion_edges]) - - model.photo_xsect = photo_xsect - model.continuum_list_nu = continuum_list_nu - model.no_of_edges = no_of_edges - - estimator_size = model.no_of_shells * no_of_edges - estims = ['photo_ion_estimator', 'stim_recomb_estimator', - 'bf_heating_estimator', 'stim_recomb_cooling_estimator' - ] - for estimator in estims: - setattr(model, estimator, (c_double * estimator_size)(*[0] * estimator_size)) - - model.photo_ion_estimator_statistics = (c_int64 * estimator_size)(*[0] * estimator_size) - return model - - -@pytest.fixture(scope='module') -def ion_edges(): - return [ - {'nu': [4.0e14, 4.1e14, 4.2e14, 4.3e14], - 'x_sect': [1.0, 0.9, 0.8, 0.7], 'no_of_points': 4}, - {'nu': [3.0e14, 3.1e14, 3.2e14, 3.3e14, 3.4e14], - 'x_sect': [1.0, 0.9, 0.8, 0.7, 0.6], 'no_of_points': 5}, - {'nu': [2.8e14, 3.0e14, 3.2e14, 3.4e14], - 'x_sect': [2.0, 1.8, 1.6, 1.4], 'no_of_points': 4} - ] - - -@pytest.fixture(scope='module') -def mock_sample_nu(): - SAMPLE_NUFUNC = CFUNCTYPE(c_double, POINTER(RPacket), - POINTER(StorageModel), POINTER(RKState)) - - def sample_nu_simple(packet, model, mt_state): - return packet.contents.nu - - return SAMPLE_NUFUNC(sample_nu_simple) - - -@pytest.fixture(scope='function') -def model_3lvlatom(model): - model.line2macro_level_upper = (c_int64 * 3)(*[2, 1, 2]) - model.macro_block_references = (c_int64 * 3)(*[0, 2, 5]) - - transition_probabilities = [ - 0.0, 0.0, 0.75, 0.25, 0.0, 0.25, 0.5, 0.25, 0.0, # shell_id = 0 - 0.0, 0.0, 1.00, 0.00, 0.0, 0.00, 0.0, 1.00, 0.0 # shell_id = 1 - ] - - nd = len(transition_probabilities)//2 - model.transition_type = (c_int64 * nd)(*[1, 1, -1, 1, 0, 0, -1, -1, 0]) - model.destination_level_id = (c_int64 * nd)(*[1, 2, 0, 2, 0, 1, 1, 0, 0]) - model.transition_line_id = (c_int64 * nd)(*[0, 1, 1, 2, 1, 2, 2, 0, 0]) - - model.transition_probabilities_nd = c_int64(nd) - model.transition_probabilities = (c_double * (nd * 2))(*transition_probabilities) - - model.last_line_interaction_out_id = (c_int64 * 1)(*[-5]) - - return model - - -def d_cont_setter(d_cont, model, packet): - model.inverse_electron_densities[packet.current_shell_id] = c_double(1.0) - model.inverse_sigma_thomson = c_double(1.0) - packet.tau_event = c_double(d_cont) - - -def d_line_setter(d_line, model, packet): - packet.mu = c_double(0.0) - scale = d_line * 1e1 - model.time_explosion = c_double(INVERSE_C * scale) - packet.nu = c_double(1.0) - nu_line = (1. - d_line/scale) - packet.nu_line = c_double(nu_line) - - -def d_boundary_setter(d_boundary, model, packet): - packet.mu = c_double(1e-16) - r_outer = 2. * d_boundary - model.r_outer[packet.current_shell_id] = r_outer - - r = np.sqrt(r_outer**2 - d_boundary**2) - packet.r = r - - - - -""" -Important Tests: ----------------- -The tests written further (till next block comment is encountered) have been -categorized as important tests, these tests correspond to methods which are -relatively old and stable code. -""" - - -@pytest.mark.parametrize( - ['x', 'x_insert', 'imin', 'imax', 'expected_params'], - [([5.0, 4.0, 3.0, 1.0], 2.0, 0, 3, - {'result': 2, 'ret_val': TARDIS_ERROR_OK}), - - ([5.0, 4.0, 3.0, 2.0], 0.0, 0, 3, - {'result': 0, 'ret_val': TARDIS_ERROR_BOUNDS_ERROR})] -) -def test_reverse_binary_search(clib, x, x_insert, imin, imax, expected_params): - x = (c_double * (imax - imin + 1))(*x) - x_insert = c_double(x_insert) - imin = c_int64(imin) - imax = c_int64(imax) - obtained_result = c_int64(0) - - clib.reverse_binary_search.restype = c_uint - obtained_tardis_error = clib.reverse_binary_search( - byref(x), x_insert, imin, imax, byref(obtained_result)) - - assert obtained_result.value == expected_params['result'] - assert obtained_tardis_error == expected_params['ret_val'] - - -@pytest.mark.parametrize( - ['nu', 'nu_insert', 'number_of_lines', 'expected_params'], - [([0.5, 0.4, 0.3, 0.1], 0.2, 4, - {'result': 3, 'ret_val': TARDIS_ERROR_OK}), - - ([0.5, 0.4, 0.3, 0.2], 0.1, 4, - {'result': 4, 'ret_val': TARDIS_ERROR_OK}), - - ([0.4, 0.3, 0.2, 0.1], 0.5, 4, - {'result': 0, 'ret_val': TARDIS_ERROR_OK})] -) -def test_line_search(clib, nu, nu_insert, number_of_lines, expected_params): - nu = (c_double * number_of_lines)(*nu) - nu_insert = c_double(nu_insert) - number_of_lines = c_int64(number_of_lines) - obtained_result = c_int64(0) - - clib.line_search.restype = c_uint - obtained_tardis_error = clib.line_search( - byref(nu), nu_insert, number_of_lines, byref(obtained_result)) - - assert obtained_result.value == expected_params['result'] - assert obtained_tardis_error == expected_params['ret_val'] - -@pytest.mark.parametrize( - ['x', 'x_insert', 'imin', 'imax', 'expected_params'], - [([2.0, 4.0, 6.0, 7.0], 5.0, 0, 3, - {'result': 2, 'ret_val': TARDIS_ERROR_OK}), - - ([2.0, 3.0, 5.0, 7.0], 8.0, 0, 3, - {'result': 0, 'ret_val': TARDIS_ERROR_BOUNDS_ERROR}), - - ([2.0, 4.0, 6.0, 7.0], 4.0, 0, 3, - {'result': 0, 'ret_val': TARDIS_ERROR_OK}) - ] -) -def test_binary_search(clib, x, x_insert, imin, imax, expected_params): - x = (c_double * (imax - imin + 1))(*x) - x_insert = c_double(x_insert) - imin = c_int64(imin) - imax = c_int64(imax) - obtained_result = c_int64(0) - - clib.binary_search.restype = c_uint - obtained_tardis_error = clib.binary_search( - byref(x), x_insert, imin, imax, byref(obtained_result)) - - assert obtained_result.value == expected_params['result'] - assert obtained_tardis_error == expected_params['ret_val'] - - -@pytest.mark.parametrize( - ['mu', 'r', 'inv_t_exp', 'expected'], - [(0.3, 7.5e14, 1 / 5.2e7, 0.9998556693818854), - (-.3, 8.1e14, 1 / 2.6e7, 1.0003117541351274)] -) -def test_rpacket_doppler_factor(clib, mu, r, inv_t_exp, expected, packet, model): - # Set the params from test cases here - packet.mu = mu - packet.r = r - model.inverse_time_explosion = inv_t_exp - - # Perform any other setups just before this, they can be additional calls - # to other methods or introduction of some temporary variables - - # Set `restype` attribute if returned quantity is used - clib.rpacket_doppler_factor.restype = c_double - # Call the C method (make sure to pass quantities as `ctypes` data types) - obtained = clib.rpacket_doppler_factor(byref(packet), byref(model)) - - # Perform required assertions - assert_almost_equal(obtained, expected) - - -@pytest.mark.parametrize( - ['packet_params', 'expected_params'], - [({'mu': 0.3, 'r': 7.5e14}, - {'d_boundary': 259376919351035.88}), - - ({'mu': -.3, 'r': 7.5e13}, - {'d_boundary': -664987228972291.5}), - - ({'mu': -.3, 'r': 7.5e14}, - {'d_boundary': 709376919351035.9})] -) -def test_compute_distance2boundary(clib, packet_params, expected_params, packet, model): - packet.mu = packet_params['mu'] - packet.r = packet_params['r'] - - clib.compute_distance2boundary(byref(packet), byref(model)) - - assert_almost_equal(packet.d_boundary, expected_params['d_boundary']) - - -# TODO: split this into two tests - one to assert errors and other for d_line -@pytest.mark.parametrize( - ['packet_params', 'expected_params'], - [({'nu_line': 0.1, 'next_line_id': 0, 'last_line': 1}, - {'tardis_error': TARDIS_ERROR_OK, 'd_line': 1e+99}), - - ({'nu_line': 0.2, 'next_line_id': 1, 'last_line': 0}, - {'tardis_error': TARDIS_ERROR_OK, 'd_line': 7.792353908000001e+17}), - - ({'nu_line': 0.5, 'next_line_id': 1, 'last_line': 0}, - {'tardis_error': TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE, 'd_line': 0.0}), - - ({'nu_line': 0.6, 'next_line_id': 0, 'last_line': 0}, - {'tardis_error': TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE, 'd_line': 0.0})] -) -def test_compute_distance2line(clib, packet_params, expected_params, packet, model): - packet.nu_line = packet_params['nu_line'] - packet.next_line_id = packet_params['next_line_id'] - packet.last_line = packet_params['last_line'] - - packet.d_line = 0.0 - clib.compute_distance2line.restype = c_uint - obtained_tardis_error = clib.compute_distance2line(byref(packet), byref(model)) - - assert_almost_equal(packet.d_line, expected_params['d_line']) - assert obtained_tardis_error == expected_params['tardis_error'] - - -@pytest.mark.parametrize( - ['packet_params', 'expected_params'], - [({'virtual_packet': 0}, - {'chi_cont': 6.652486e-16, 'd_cont': 4.359272608766106e+28}), - - ({'virtual_packet': 1}, - {'chi_cont': 6.652486e-16, 'd_cont': 1e+99})] -) -def test_compute_distance2continuum(clib, packet_params, expected_params, packet, model): - packet.virtual_packet = packet_params['virtual_packet'] - - clib.compute_distance2continuum(byref(packet), byref(model)) - - assert_almost_equal(packet.chi_cont, expected_params['chi_cont']) - assert_almost_equal(packet.d_cont, expected_params['d_cont']) - - -@pytest.mark.parametrize('full_relativity', [1, 0]) -@pytest.mark.parametrize( - ['packet_params', 'expected_params'], - [({'nu': 0.4, 'mu': 0.3, 'energy': 0.9, 'r': 7.5e14}, - {'mu': 0.3120599529139568, 'r': 753060422542573.9, - 'j': 8998701024436.969, 'nubar': 3598960894542.354}), - - ({'nu': 0.6, 'mu': -.5, 'energy': 0.5, 'r': 8.1e14}, - {'mu': -.4906548373534084, 'r': 805046582503149.2, - 'j': 5001298975563.031, 'nubar': 3001558973156.1387})] -) -def test_move_packet(clib, packet_params, expected_params, - packet, model, full_relativity): - packet.nu = packet_params['nu'] - packet.mu = packet_params['mu'] - packet.energy = packet_params['energy'] - packet.r = packet_params['r'] - model.full_relativity = full_relativity - - clib.rpacket_doppler_factor.restype = c_double - doppler_factor = clib.rpacket_doppler_factor(byref(packet), byref(model)) - clib.move_packet(byref(packet), byref(model), c_double(1.e13)) - - assert_almost_equal(packet.mu, expected_params['mu']) - assert_almost_equal(packet.r, expected_params['r']) - - expected_j = expected_params['j'] - expected_nubar = expected_params['nubar'] - if full_relativity: - expected_j *= doppler_factor - expected_nubar *= doppler_factor - - assert_allclose(model.js[packet.current_shell_id], - expected_j, rtol=5e-7) - assert_allclose(model.nubars[packet.current_shell_id], - expected_nubar, rtol=5e-7) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['packet_params', 'j_blue_idx', 'expected'], - [({'nu': 0.30, 'energy': 0.30}, 0, 1.0), - ({'nu': 0.20, 'energy': 1.e5}, 0, 5e5), - ({'nu': 2e15, 'energy': 0.50}, 1, 2.5e-16), - ({'nu': 0.40, 'energy': 1e-7}, 1, 2.5e-7)], -) -def test_increment_j_blue_estimator_full_relativity(clib, packet_params, - j_blue_idx, expected, - packet, model): - packet.nu = packet_params['nu'] - packet.energy = packet_params['energy'] - model.full_relativity = True - - clib.increment_j_blue_estimator(byref(packet), byref(model), - c_double(packet.d_line), - c_int64(j_blue_idx)) - - assert_almost_equal(model.line_lists_j_blues[j_blue_idx], expected) - - -@pytest.mark.parametrize( - ['packet_params', 'j_blue_idx', 'expected'], - [({'nu': 0.1, 'mu': 0.3, 'r': 7.5e14}, 0, 8.998643292289723), - ({'nu': 0.2, 'mu': -.3, 'r': 7.7e14}, 0, 4.499971133976377), - ({'nu': 0.5, 'mu': 0.5, 'r': 7.9e14}, 1, 0.719988453650551), - ({'nu': 0.6, 'mu': -.5, 'r': 8.1e14}, 1, 0.499990378058792)] -) -def test_increment_j_blue_estimator(clib, packet_params, j_blue_idx, expected, packet, model): - packet.nu = packet_params['nu'] - packet.mu = packet_params['mu'] - packet.r = packet_params['r'] - - clib.compute_distance2line(byref(packet), byref(model)) - clib.move_packet(byref(packet), byref(model), c_double(1.e13)) - clib.increment_j_blue_estimator(byref(packet), byref(model), - c_double(packet.d_line), c_int64(j_blue_idx)) - - assert_almost_equal(model.line_lists_j_blues[j_blue_idx], expected) - - -@pytest.mark.parametrize( - ['packet_params', 'expected_params'], - [({'virtual_packet': 0, 'current_shell_id': 0, 'next_shell_id': 1}, - {'status': TARDIS_PACKET_STATUS_IN_PROCESS, 'current_shell_id': 1}), - - ({'virtual_packet': 1, 'current_shell_id': 1, 'next_shell_id': 1}, - {'status': TARDIS_PACKET_STATUS_EMITTED, 'current_shell_id': 1, - 'tau_event': 29000000000000.008}), - - ({'virtual_packet': 1, 'current_shell_id': 0, 'next_shell_id': -1}, - {'status': TARDIS_PACKET_STATUS_REABSORBED, 'current_shell_id': 0, - 'tau_event': 29000000000000.008})] -) -def test_move_packet_across_shell_boundary(clib, packet_params, expected_params, - packet, model, mt_state): - packet.virtual_packet = packet_params['virtual_packet'] - packet.current_shell_id = packet_params['current_shell_id'] - packet.next_shell_id = packet_params['next_shell_id'] - - clib.move_packet_across_shell_boundary(byref(packet), byref(model), - c_double(1.e13), byref(mt_state)) - - if packet_params['virtual_packet'] == 1: - assert_almost_equal(packet.tau_event, expected_params['tau_event']) - assert packet.status == expected_params['status'] - assert packet.current_shell_id == expected_params['current_shell_id'] - - -@pytest.mark.parametrize( - ['packet_params', 'expected_params'], - [({'nu': 0.4, 'mu': 0.3, 'energy': 0.9, 'r': 7.5e14}, - {'nu': 0.39974659819356556, 'energy': 0.8994298459355226}), - - ({'nu': 0.6, 'mu': -.5, 'energy': 0.5, 'r': 8.1e14}, - {'nu': 0.5998422620533325, 'energy': 0.4998685517111104})] -) -def test_montecarlo_thomson_scatter(clib, packet_params, expected_params, packet, - model, mt_state): - packet.nu = packet_params['nu'] - packet.mu = packet_params['mu'] - packet.energy = packet_params['energy'] - packet.r = packet_params['r'] - - clib.montecarlo_thomson_scatter(byref(packet), byref(model), - c_double(1.e13), byref(mt_state)) - - assert_almost_equal(packet.nu, expected_params['nu']) - assert_almost_equal(packet.energy, expected_params['energy']) - - -@pytest.mark.parametrize( - ['packet_params', 'expected_params'], - # TODO: Add scientifically sound test cases. - [({'virtual_packet': 1, 'tau_event': 2.9e13, 'last_line': 0}, - {'tau_event': 2.9e13, 'next_line_id': 2}), - - ({'virtual_packet': 0, 'tau_event': 2.9e13, 'last_line': 0}, - {'tau_event': 2.9e13, 'next_line_id': 2}), - - ({'virtual_packet': 0, 'tau_event': 2.9e13, 'last_line': 0}, - {'tau_event': 2.9e13, 'next_line_id': 2}), - ] -) -def test_montecarlo_line_scatter(clib, packet_params, expected_params, packet, model, mt_state): - packet.virtual_packet = packet_params['virtual_packet'] - packet.tau_event = packet_params['tau_event'] - packet.last_line = packet_params['last_line'] - - clib.montecarlo_line_scatter(byref(packet), byref(model), - c_double(1.e13), byref(mt_state)) - - assert_almost_equal(packet.tau_event, expected_params['tau_event']) - assert_almost_equal(packet.next_line_id, expected_params['next_line_id']) - - -@pytest.mark.parametrize( - ['distances', 'expected'], - [({'boundary': 1.3e13, 'continuum': 1e14, 'line': 1e15}, - {'handler': 'move_packet_across_shell_boundary', 'distance': 1.3e13}), - - ({'boundary': 1.3e13, 'continuum': 1e14, 'line': 2.5e12}, - {'handler': 'montecarlo_line_scatter', 'distance': 2.5e12}), - - ({'boundary': 1.3e13, 'continuum': 1e11, 'line': 2.5e12}, - {'handler': 'montecarlo_thomson_scatter', 'distance': 1e11})] -) -def test_get_event_handler(clib, packet, model, mt_state, distances, expected): - d_cont_setter(distances['continuum'], model, packet) - d_line_setter(distances['line'], model, packet) - d_boundary_setter(distances['boundary'], model, packet) - obtained_distance = c_double() - - clib.get_event_handler.restype = c_void_p - obtained_handler = clib.get_event_handler(byref(packet), byref(model), - byref(obtained_distance), - byref(mt_state)) - - expected_handler = getattr(clib, expected['handler']) - expected_handler = cast(expected_handler, c_void_p).value - - assert_equal(obtained_handler, expected_handler) - assert_allclose(obtained_distance.value, expected['distance'], rtol=1e-10) - - -@pytest.mark.parametrize( - ['z_random', 'packet_params', 'expected'], - [(0.22443743797312765, - {'activation_level': 1, 'shell_id': 0}, 1), # Direct deactivation - - (0.78961460371187597, # next z_random = 0.818455414618 - {'activation_level': 1, 'shell_id': 0}, 0), # Upwards jump, then deactivation - - (0.22443743797312765, # next z_random = 0.545678896748 - {'activation_level': 2, 'shell_id': 0}, 1), # Downwards jump, then deactivation - - (0.765958602560605, # next z_random = 0.145914243888, 0.712382380384 - {'activation_level': 1, 'shell_id': 0}, 1), # Upwards jump, downwards jump, then deactivation - - (0.22443743797312765, - {'activation_level': 2, 'shell_id': 1}, 0)] # Direct deactivation -) -def test_macro_atom(clib, model_3lvlatom, packet, z_random, packet_params, get_rkstate, expected): - packet.macro_atom_activation_level = packet_params['activation_level'] - packet.current_shell_id = packet_params['shell_id'] - rkstate = get_rkstate(z_random) - - clib.macro_atom(byref(packet), byref(model_3lvlatom), byref(rkstate)) - obtained_line_id = model_3lvlatom.last_line_interaction_out_id[packet.id] - - assert_equal(obtained_line_id, expected) - - - -""" -Simple Tests: ----------------- -These test check very simple pices of code still work. -""" - -@pytest.mark.parametrize( - ['packet_params', 'line_idx', 'expected'], - [({'energy': 0.0}, 0, 0), - ({'energy': 1.0}, 1, 1), - ({'energy': 0.5}, 2, 1.5)] -) -@pytest.mark.skipif(True, reason="Needs rewrite to be relevant.") -def test_increment_Edotlu_estimator(clib, packet_params, line_idx, expected, packet, model): - packet.energy = packet_params['energy'] - - clib.increment_Edotlu_estimator(byref(packet), byref(model), c_int64(line_idx)) - - assert_almost_equal(model.line_lists_Edotlu[line_idx], expected) - - -""" -Difficult Tests: ----------------- -The tests written further are more complex than previous tests. They require -proper design procedure. They are not taken up yet and intended to be -completed together in future. -""" - - - -@pytest.mark.skipif(True, reason="Yet to be written.") -def test_montecarlo_one_packet(packet, model, mt_state): - pass - - -@pytest.mark.skipif(True, reason="Yet to be written.") -def test_montecarlo_one_packet_loop(packet, model, mt_state): - pass - - -@pytest.mark.skipif(True, reason="Yet to be written.") -def test_montecarlo_main_loop(packet, model, mt_state): - pass - - - -""" -Continuum Tests: ----------------- -The tests written further (till next block comment is encountered) are for the -methods related to continuum interactions. -""" - - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - 't_electron', [2500., 15000.] -) -def test_sample_nu_free_free(clib, t_electron, packet, model, mt_state_seeded, expected_ff_emissivity): - model.t_electrons[packet.current_shell_id] = t_electron - clib.sample_nu_free_free.restype = c_double - - nu_bins, expected_emissivity = expected_ff_emissivity(t_electron) - - nus = [] - for _ in range(int(1e5)): - nu = clib.sample_nu_free_free(byref(packet), byref(model), byref(mt_state_seeded)) - nus.append(nu) - - obtained_emissivity, _ = np.histogram(nus, density=True, bins=nu_bins) - - assert_allclose(obtained_emissivity, expected_emissivity, rtol=1e-10) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['packet_params', 't_electrons', 'chi_ff_factor', 'expected'], - [({'nu': 4.5e14, 'mu': 0.0, 'current_shell_id': 1}, 15000, 2.0, 1.6746639430359494e-44), - ({'nu': 3.0e15, 'mu': 0.0, 'current_shell_id': 0}, 5000, 3.0, 1.1111111111107644e-46), - ({'nu': 3.0e15, 'mu': 0.4, 'current_shell_id': 0}, 10000, 4.0, 1.5638286016098277e-46)] -) -def test_calculate_chi_ff(clib, packet, model, packet_params, t_electrons, chi_ff_factor, expected): - packet.mu = packet_params['mu'] - packet.nu = packet_params['nu'] - packet.current_shell_id = packet_params['current_shell_id'] - packet.r = 1.04e17 - - model.t_electrons[packet_params['current_shell_id']] = t_electrons - model.chi_ff_factor[packet_params['current_shell_id']] = chi_ff_factor - - clib.calculate_chi_ff(byref(packet), byref(model)) - obtained = packet.chi_ff - - assert_equal(obtained, expected) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['continuum_status', 'z_random', 'packet_params', 'expected'], - [(CONTINUUM_OFF, 0.94183547596539363, - {'chi_c': 1.0, 'chi_th': 0.4, 'chi_bf': 0.5}, - 'montecarlo_thomson_scatter'), - - (CONTINUUM_ON, 0.22443743797312765, - {'chi_c': 1.0, 'chi_th': 0.4, 'chi_bf': 0.5}, - 'montecarlo_thomson_scatter'), - - (CONTINUUM_ON, 0.54510721066252377, - {'chi_c': 1.0, 'chi_th': 0.4, 'chi_bf': 0.5}, - 'montecarlo_bound_free_scatter'), - - (CONTINUUM_ON, 0.94183547596539363, - {'chi_c': 1.0, 'chi_th': 0.4, 'chi_bf': 0.5}, - 'montecarlo_free_free_scatter'), - - (CONTINUUM_ON, 0.22443743797312765, - {'chi_c': 1e2, 'chi_th': 1e1, 'chi_bf': 2e1}, - 'montecarlo_bound_free_scatter')] -) -def test_montecarlo_continuum_event_handler(clib, continuum_status, expected, z_random, - packet_params, packet, model, get_rkstate): - packet.chi_cont = packet_params['chi_c'] - packet.chi_th = packet_params['chi_th'] - packet.chi_bf = packet_params['chi_bf'] - model.cont_status = continuum_status - - rkstate = get_rkstate(z_random) - - clib.montecarlo_continuum_event_handler.restype = c_void_p - obtained = clib.montecarlo_continuum_event_handler(byref(packet), - byref(model), byref(rkstate)) - expected = getattr(clib, expected) - expected = cast(expected, c_void_p).value - - assert_equal(obtained, expected) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['nu', 'continuum_id', 'expected', 'bf_treatment'], - [(4.40e14, 1, 0.00, BoundFreeTreatment.LIN_INTERPOLATION), - (3.25e14, 1, 0.75, BoundFreeTreatment.LIN_INTERPOLATION), - (4.03e14, 0, 0.97, BoundFreeTreatment.LIN_INTERPOLATION), - (4.10e14 + 1e-1, 0, 0.90, BoundFreeTreatment.LIN_INTERPOLATION), - pytest.param(4.1e14, 0, 0.90, BoundFreeTreatment.LIN_INTERPOLATION, - marks=pytest.mark.xfail), - (6.50e14, 0, 0.23304506144742834, BoundFreeTreatment.HYDROGENIC), - (3.40e14, 2, 1.1170364339507428, BoundFreeTreatment.HYDROGENIC)] -) -def test_bf_cross_section(clib, nu, continuum_id, model_w_edges, expected, bf_treatment): - model_w_edges.bf_treatment = bf_treatment.value - - clib.bf_cross_section.restype = c_double - obtained = clib.bf_cross_section(byref(model_w_edges), continuum_id, c_double(nu)) - - assert_almost_equal(obtained, expected) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['packet_params', 'expected'], - [({'nu': 4.13e14, 'mu': 0.0, 'current_shell_id': 1}, - [3.2882087455641473, 0.0, 0.0]), - - ({'nu': 3.27e14, 'mu': 0.0, 'current_shell_id': 0}, - [0.0, 1.3992114634681028, 5.702548202131454]), - - ({'nu': 3.27e14, 'mu': -0.4, 'current_shell_id': 0}, - [0.0, 1.2670858, 5.4446587])] -) -def test_calculate_chi_bf(clib, packet_params, expected, packet, model_w_edges): - model_w_edges.l_pop = (c_double * 6)(*range(1, 7)) - model_w_edges.l_pop_r = (c_double * 6)(*np.linspace(0.1, 0.6, 6)) - model_w_edges.t_electrons[packet_params['current_shell_id']] = 1e4 - - packet.mu = packet_params['mu'] - packet.nu = packet_params['nu'] - packet.r = 1.04e17 - packet.current_shell_id = packet_params['current_shell_id'] - packet.chi_bf_tmp_partial = (c_double * model_w_edges.no_of_edges)() - - clib.calculate_chi_bf(byref(packet), byref(model_w_edges)) - - obtained_chi_bf_tmp = np.ctypeslib.as_array(packet.chi_bf_tmp_partial, shape=(model_w_edges.no_of_edges,)) - expected_chi_bf_tmp = np.array(expected) - expected_chi_bf = expected_chi_bf_tmp[expected_chi_bf_tmp > 0][-1] - - assert_almost_equal(obtained_chi_bf_tmp, expected_chi_bf_tmp) - assert_almost_equal(packet.chi_bf, expected_chi_bf) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['comov_energy', 'distance', 'chi_ff', 'no_of_updates', 'expected'], - [(1.3, 1.3e14, 3e-12, 1, 507.), - (0.9, 0.7e15, 2e-12, 10, 1.260e4), - (0.8, 0.8e13, 1.5e-12, 35, 336.)] -) -def test_increment_continuum_estimators_ff_heating_estimator(clib, packet, model_w_edges, comov_energy, distance, - chi_ff, no_of_updates, expected): - packet.chi_ff = chi_ff - - for _ in range(no_of_updates): - clib.increment_continuum_estimators(byref(packet), byref(model_w_edges), c_double(distance), - c_double(0), c_double(comov_energy)) - obtained = model_w_edges.ff_heating_estimator[packet.current_shell_id] - - assert_almost_equal(obtained, expected) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['comov_nus', 'expected'], - [([4.05e14, 4.17e14, 3.3e14, 3.2e14, 2.9e14], [2, 2, 3]), - ([4.15e15, 3.25e14, 3.3e14, 2.85e14, 2.9e14], [0, 2, 4])] -) -def test_increment_continuum_estimators_photo_ion_estimator_statistics(clib, packet, model_w_edges, comov_nus, expected): - for comov_nu in comov_nus: - clib.increment_continuum_estimators(byref(packet), byref(model_w_edges), c_double(1e13), - c_double(comov_nu), c_double(1.0)) - - no_of_edges = model_w_edges.no_of_edges - no_of_shells = model_w_edges.no_of_shells - - obtained = np.ctypeslib.as_array(model_w_edges.photo_ion_estimator_statistics, - shape=(no_of_edges * no_of_shells,)) - obtained = np.reshape(obtained, newshape=(no_of_shells, no_of_edges), order='F') - obtained = obtained[packet.current_shell_id] - expected = np.array(expected) - - assert_array_equal(obtained, expected) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['comov_energy', 'distance', 'comov_nus', 'expected'], - [(1.3, 1.3e14, [4.05e14, 2.65e14], - {"photo_ion": [0.39641975308641975, 0., 0.], - "stim_recomb": [0.056757061269242064, 0., 0.], - "bf_heating": [1.9820987654321076e12, 0., 0.], - "stim_recomb_cooling": [283784812699.75476, 0., 0.]}), - - (0.9, 0.7e15, [3.25e14, 2.85e14], - {"photo_ion": [0., 1.4538461538461538, 7.315141700404858], - "stim_recomb": [0., 0.3055802, 1.7292954], - "bf_heating": [0., 36346153846153.82, 156760323886639.69], - "stim_recomb_cooling": [0., 7639505724285.9746, 33907776077426.875]})] -) -def test_increment_continuum_estimators_bf_estimators(clib, packet, model_w_edges, comov_energy, - distance, comov_nus, expected): - for comov_nu in comov_nus: - clib.increment_continuum_estimators(byref(packet), byref(model_w_edges), c_double(distance), - c_double(comov_nu), c_double(comov_energy)) - - no_of_edges = model_w_edges.no_of_edges - no_of_shells = model_w_edges.no_of_shells - - for estim_name, expected_value in expected.items(): - obtained = np.ctypeslib.as_array(getattr(model_w_edges, estim_name + "_estimator"), - shape=(no_of_edges * no_of_shells,)) - obtained = np.reshape(obtained, newshape=(no_of_shells, no_of_edges), order='F') - obtained = obtained[packet.current_shell_id] - - assert_almost_equal(obtained, np.array(expected_value)) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['packet_params', 'expected_params'], - [({'nu_comov': 1.1e16, 'mu': 0.0, 'r': 1.4e14}, - {'next_line_id': 5, 'last_line': True, 'nu': 1.1e16, 'type_id': 3}), - - ({'nu_comov': 1.3e16, 'mu': 0.3, 'r': 7.5e14}, - {'next_line_id': 0, 'last_line': False, 'nu': 1.30018766e+16, 'type_id': 3}), - - ({'nu_comov': 1.24e16, 'mu': -0.3, 'r': 7.5e14}, - {'next_line_id': 2, 'last_line': False, 'nu': 1.23982106e+16, 'type_id': 4})] -) -def test_continuum_emission(clib, packet, model, mock_sample_nu, packet_params, expected_params, mt_state): - packet.nu = packet_params['nu_comov'] # Is returned by mock function mock_sample_nu - packet.mu = packet_params['mu'] - packet.r = packet_params['r'] - expected_interaction_out_type = expected_params['type_id'] - - clib.continuum_emission(byref(packet), byref(model), byref(mt_state), - mock_sample_nu, expected_interaction_out_type) - - obtained_next_line_id = packet.next_line_id - obtained_last_interaction_out_type = model.last_interaction_out_type[0] - - assert_equal(obtained_next_line_id, expected_params['next_line_id']) - assert_equal(packet.last_line, expected_params['last_line']) - assert_equal(expected_interaction_out_type, obtained_last_interaction_out_type) - assert_allclose(packet.nu, expected_params['nu'], rtol=1e-7) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['packet_params', 'expected'], - [({'next_line_id': 3, 'last_line': 0}, 1), - ({'next_line_id': 5, 'last_line': 1}, 0), - ({'next_line_id': 2, 'last_line': 0}, 0), - ({'next_line_id': 1, 'last_line': 0}, 1)] -) -def test_test_for_close_line(clib, packet, model, packet_params, expected): - packet.nu_line = model.line_list_nu[packet_params['next_line_id'] - 1] - packet.next_line_id = packet_params['next_line_id'] - - clib.test_for_close_line(byref(packet), byref(model)) - - assert_equal(expected, packet.close_line) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['packet_params', 'z_random', 'expected'], - [({'current_continuum_id': 1, 'chi_bf_tmp_partial': [0.0, 0.23e13, 1.0e13]}, - 0.22443743797312765, 1), - - ({'current_continuum_id': 1, 'chi_bf_tmp_partial': [0.0, 0.23e10, 1.0e10]}, - 0.78961460371187597, 2), - - ({'current_continuum_id': 0, 'chi_bf_tmp_partial': [0.2e5, 0.5e5, 0.6e5, 1.0e5, 1.0e5]}, - 0.78961460371187597, 3)] -) -def test_montecarlo_bound_free_scatter_continuum_selection(clib, packet, model_3lvlatom, packet_params, - get_rkstate, z_random, expected): - rkstate = get_rkstate(z_random) - packet.current_continuum_id = packet_params['current_continuum_id'] - - chi_bf_tmp = packet_params['chi_bf_tmp_partial'] - packet.chi_bf_tmp_partial = (c_double * len(chi_bf_tmp))(*chi_bf_tmp) - packet.chi_bf = chi_bf_tmp[-1] - model_3lvlatom.no_of_edges = len(chi_bf_tmp) - - clib.montecarlo_bound_free_scatter(byref(packet), byref(model_3lvlatom), - c_double(1.e13), byref(rkstate)) - - assert_equal(packet.current_continuum_id, expected) - assert_equal(model_3lvlatom.last_line_interaction_in_id[packet.id], expected) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['mu', 'r', 'inv_t_exp', 'full_relativity'], - [(0.8, 7.5e14, 1 / 5.2e5, 1), - (-0.7, 7.5e14, 1 / 5.2e5, 1), - (0.3, 7.5e14, 1 / 2.2e5, 1), - (0.0, 7.5e14, 1 / 2.2e5, 1), - (-0.7, 7.5e14, 1 / 5.2e5, 0)] -) -def test_frame_transformations(clib, packet, model, mu, r, - inv_t_exp, full_relativity): - packet.r = r - packet.mu = mu - model.inverse_time_explosion = inv_t_exp - model.full_relativity = full_relativity - clib.rpacket_doppler_factor.restype = c_double - clib.rpacket_inverse_doppler_factor.restype = c_double - - inverse_doppler_factor = clib.rpacket_inverse_doppler_factor(byref(packet), byref(model)) - clib.angle_aberration_CMF_to_LF(byref(packet), byref(model)) - - doppler_factor = clib.rpacket_doppler_factor(byref(packet), byref(model)) - - assert_almost_equal(doppler_factor * inverse_doppler_factor, 1.0) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['mu', 'r', 'inv_t_exp'], - [(0.8, 7.5e14, 1 / 5.2e5), - (-0.7, 7.5e14, 1 / 5.2e5), - (0.3, 7.5e14, 1 / 2.2e5), - (0.0, 7.5e14, 1 / 2.2e5), - (-0.7, 7.5e14, 1 / 5.2e5)] -) -def test_angle_transformation_invariance(clib, packet, model, - mu, r, inv_t_exp): - packet.r = r - packet.mu = mu - model.inverse_time_explosion = inv_t_exp - model.full_relativity = 1 - clib.angle_aberration_LF_to_CMF.restype = c_double - - clib.angle_aberration_CMF_to_LF(byref(packet), byref(model)) - mu_obtained = clib.angle_aberration_LF_to_CMF( - byref(packet), byref(model), c_double(packet.mu)) - - assert_almost_equal(mu_obtained, mu) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - 'full_relativity', - [1, 0] -) -@pytest.mark.parametrize( - ['mu', 'r', 't_exp', 'nu', 'nu_line'], - [(0.8, 7.5e14, 5.2e5, 1.0e15, 9.4e14), - (0.0, 6.3e14, 2.2e5, 6.0e12, 5.8e12), - (1.0, 9.0e14, 2.2e5, 4.0e8, 3.4e8), - (0.9, 9.0e14, 0.5e5, 1.0e15, 4.5e14), - (-0.7, 7.5e14, 5.2e5, 1.0e15, 9.8e14), - (-1.0, 6.3e14, 2.2e5, 6.0e12, 6.55e12)] -) -def test_compute_distance2line_relativistic(clib, mu, r, t_exp, nu, nu_line, - full_relativity, packet, model): - packet.r = r - packet.mu = mu - packet.nu = nu - packet.nu_line = nu_line - model.inverse_time_explosion = 1 / t_exp - model.time_explosion = t_exp - model.full_relativity = full_relativity - - clib.rpacket_doppler_factor.restype = c_double - - clib.compute_distance2line(byref(packet), byref(model)) - clib.move_packet(byref(packet), byref(model), c_double(packet.d_line)) - - doppler_factor = clib.rpacket_doppler_factor(byref(packet), byref(model)) - comov_nu = packet.nu * doppler_factor - - assert_allclose(comov_nu, nu_line, rtol=1e-14) - - -@pytest.mark.continuumtest -@pytest.mark.skipif(True, reason="Yet to be written.") -def test_montecarlo_free_free_scatter(packet, model, mt_state): - pass diff --git a/tardis/montecarlo/tests/test_montecarlo.py b/tardis/montecarlo/tests/test_montecarlo.py index e69de29bb2d..86bc15d72d7 100644 --- a/tardis/montecarlo/tests/test_montecarlo.py +++ b/tardis/montecarlo/tests/test_montecarlo.py @@ -0,0 +1,1064 @@ +""" +Unit tests for methods in `tardis/montecarlo/src/cmontecarlo.c`. +* `ctypes` library is used to wrap C methods and expose them to python. + + +Probable Reasons for Failing Tests: +----------------------------------- + +1. Change made in C struct declarations: + - Reflect the changes done in C structs, into Python counterparts. + - Check **tardis/montecarlo/struct.py**. + +2. Return type of any method changed: + - Modify the `restype` parameter in the test method here. + - For example: + ``` + clib.rpacket_doppler_factor.restype = c_double + ``` + +3. Underlying logic modified: + - Check whether the changes made in C method are logically correct. + - If the changes made were correct and necessary, update the corresponding + test case. + + +General Test Design Procedure: +------------------------------ + +Please follow this design procedure while adding a new test: + +1. Parametrization as per desire of code coverage. + - C tests have different flows controlled by conditional statements. + Parameters checked in conditions can be provided in different testcases. + - Keep consistency with variable names as (in order): + - `packet_params` + - `model_params` + - `expected_params` (`expected` if only one value to be asserted.) + - Suggested variable names can be compromised if readability of the test + increases. + +2. Test Method body: + - Keep name as `test_` + `(name of C method)`. + - Refer to method `test_rpacket_doppler_factor` below for description. +""" + + +import os +import pytest +import numpy as np +import pandas as pd +import tardis.montecarlo.formal_integral as formal_integral +import tardis.montecarlo.montecarlo_numba.r_packet as r_packet + + + +from ctypes import ( + CDLL, + byref, + c_uint, + c_int64, + c_double, + c_ulong, + c_void_p, + cast, + POINTER, + pointer, + CFUNCTYPE + ) + +from numpy.testing import ( + assert_equal, + assert_almost_equal, + assert_array_equal, + assert_allclose + ) + +from tardis import __path__ as path +from tardis.montecarlo.struct import ( + RPacket, StorageModel, RKState, + PhotoXsect1level, + TARDIS_ERROR_OK, + TARDIS_ERROR_BOUNDS_ERROR, + TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE, + TARDIS_PACKET_STATUS_IN_PROCESS, + TARDIS_PACKET_STATUS_EMITTED, + TARDIS_PACKET_STATUS_REABSORBED, + CONTINUUM_OFF, + CONTINUUM_ON, + INVERSE_C, + BoundFreeTreatment +) + + +@pytest.fixture(scope='module') +def continuum_compare_data_fname(): + fname = 'continuum_compare_data.hdf' + return os.path.join(path[0], 'montecarlo', 'tests', 'data', fname) + + +@pytest.fixture(scope='module') +def continuum_compare_data(continuum_compare_data_fname, request): + compare_data = pd.HDFStore(continuum_compare_data_fname, mode='r') + + def fin(): + compare_data.close() + request.addfinalizer(fin) + + return compare_data + + +@pytest.fixture(scope="function") +def expected_ff_emissivity(continuum_compare_data): + emissivities = continuum_compare_data['ff_emissivity'] + + def ff_emissivity(t_electron): + emissivity = emissivities[t_electron] + nu_bins = emissivity['nu_bins'].values + emissivity_value = emissivity['emissivity'].dropna().values + + return nu_bins, emissivity_value + + return ff_emissivity + + +@pytest.fixture(scope='module') +def get_rkstate(continuum_compare_data): + data = continuum_compare_data['z2rkstate_key'] + pos_data = continuum_compare_data['z2rkstate_pos'] + + def z2rkstate(z_random): + key = (c_ulong * 624)(*data.loc[z_random].values) + pos = pos_data.loc[z_random] + return RKState( + key=key, + pos=pos, + has_gauss=0, + gauss=0.0 + ) + + return z2rkstate + + +@pytest.fixture(scope='function') +def model_w_edges(ion_edges, model): + photo_xsect = (POINTER(PhotoXsect1level) * len(ion_edges))() + + for i, edge in enumerate(ion_edges): + x_sect_1level = PhotoXsect1level() + for key, value in edge.items(): + if key in ['nu', 'x_sect']: + value = (c_double * len(value))(*value) + setattr(x_sect_1level, key, value) + photo_xsect[i] = pointer(x_sect_1level) + + no_of_edges = len(ion_edges) + continuum_list_nu = (c_double * no_of_edges)(*[edge['nu'][0] for edge in ion_edges]) + + model.photo_xsect = photo_xsect + model.continuum_list_nu = continuum_list_nu + model.no_of_edges = no_of_edges + + estimator_size = model.no_of_shells * no_of_edges + estims = ['photo_ion_estimator', 'stim_recomb_estimator', + 'bf_heating_estimator', 'stim_recomb_cooling_estimator' + ] + for estimator in estims: + setattr(model, estimator, (c_double * estimator_size)(*[0] * estimator_size)) + + model.photo_ion_estimator_statistics = (c_int64 * estimator_size)(*[0] * estimator_size) + return model + + +@pytest.fixture(scope='module') +def ion_edges(): + return [ + {'nu': [4.0e14, 4.1e14, 4.2e14, 4.3e14], + 'x_sect': [1.0, 0.9, 0.8, 0.7], 'no_of_points': 4}, + {'nu': [3.0e14, 3.1e14, 3.2e14, 3.3e14, 3.4e14], + 'x_sect': [1.0, 0.9, 0.8, 0.7, 0.6], 'no_of_points': 5}, + {'nu': [2.8e14, 3.0e14, 3.2e14, 3.4e14], + 'x_sect': [2.0, 1.8, 1.6, 1.4], 'no_of_points': 4} + ] + + +@pytest.fixture(scope='module') +def mock_sample_nu(): + SAMPLE_NUFUNC = CFUNCTYPE(c_double, POINTER(RPacket), + POINTER(StorageModel), POINTER(RKState)) + + def sample_nu_simple(packet, model, mt_state): + return packet.contents.nu + + return SAMPLE_NUFUNC(sample_nu_simple) + + +@pytest.fixture(scope='function') +def model_3lvlatom(model): + model.line2macro_level_upper = (c_int64 * 3)(*[2, 1, 2]) + model.macro_block_references = (c_int64 * 3)(*[0, 2, 5]) + + transition_probabilities = [ + 0.0, 0.0, 0.75, 0.25, 0.0, 0.25, 0.5, 0.25, 0.0, # shell_id = 0 + 0.0, 0.0, 1.00, 0.00, 0.0, 0.00, 0.0, 1.00, 0.0 # shell_id = 1 + ] + + nd = len(transition_probabilities)//2 + model.transition_type = (c_int64 * nd)(*[1, 1, -1, 1, 0, 0, -1, -1, 0]) + model.destination_level_id = (c_int64 * nd)(*[1, 2, 0, 2, 0, 1, 1, 0, 0]) + model.transition_line_id = (c_int64 * nd)(*[0, 1, 1, 2, 1, 2, 2, 0, 0]) + + model.transition_probabilities_nd = c_int64(nd) + model.transition_probabilities = (c_double * (nd * 2))(*transition_probabilities) + + model.last_line_interaction_out_id = (c_int64 * 1)(*[-5]) + + return model + + +def d_cont_setter(d_cont, model, packet): + model.inverse_electron_densities[packet.current_shell_id] = c_double(1.0) + model.inverse_sigma_thomson = c_double(1.0) + packet.tau_event = c_double(d_cont) + + +def d_line_setter(d_line, model, packet): + packet.mu = c_double(0.0) + scale = d_line * 1e1 + model.time_explosion = c_double(INVERSE_C * scale) + packet.nu = c_double(1.0) + nu_line = (1. - d_line/scale) + packet.nu_line = c_double(nu_line) + + +def d_boundary_setter(d_boundary, model, packet): + packet.mu = c_double(1e-16) + r_outer = 2. * d_boundary + model.r_outer[packet.current_shell_id] = r_outer + + r = np.sqrt(r_outer**2 - d_boundary**2) + packet.r = r + + + + +""" +Important Tests: +---------------- +The tests written further (till next block comment is encountered) have been +categorized as important tests, these tests correspond to methods which are +relatively old and stable code. +""" + + +@pytest.mark.parametrize( + ['x', 'x_insert', 'imin', 'imax', 'expected_params'], + [([5.0, 4.0, 3.0, 1.0], 2.0, 0, 3, + {'result': 2, 'ret_val': TARDIS_ERROR_OK}), + + ([5.0, 4.0, 3.0, 2.0], 0.0, 0, 3, + {'result': 0, 'ret_val': TARDIS_ERROR_BOUNDS_ERROR})] +) +def test_reverse_binary_search(x, x_insert, imin, imax, expected_params): + x = (c_double * (imax - imin + 1))(*x) + obtained_result = 0 + + obtained_tardis_error = formal_integral.reverse_binary_search( + x, x_insert, imin, imax, obtained_result) + + assert obtained_result.value == expected_params['result'] + assert obtained_tardis_error == expected_params['ret_val'] + + +@pytest.mark.parametrize( + ['nu', 'nu_insert', 'number_of_lines', 'expected_params'], + [([0.5, 0.4, 0.3, 0.1], 0.2, 4, + {'result': 3, 'ret_val': TARDIS_ERROR_OK}), + + ([0.5, 0.4, 0.3, 0.2], 0.1, 4, + {'result': 4, 'ret_val': TARDIS_ERROR_OK}), + + ([0.4, 0.3, 0.2, 0.1], 0.5, 4, + {'result': 0, 'ret_val': TARDIS_ERROR_OK})] +) +def test_line_search(nu, nu_insert, number_of_lines, expected_params): + nu = (c_double * number_of_lines)(*nu) + nu_insert = c_double(nu_insert) + obtained_result = 0 + + obtained_tardis_error = formal_integral.line_search( + nu, nu_insert, number_of_lines, obtained_result) + + assert obtained_result.value == expected_params['result'] + assert obtained_tardis_error == expected_params['ret_val'] + +@pytest.mark.parametrize( + ['x', 'x_insert', 'imin', 'imax', 'expected_params'], + [([2.0, 4.0, 6.0, 7.0], 5.0, 0, 3, + {'result': 2, 'ret_val': TARDIS_ERROR_OK}), + + ([2.0, 3.0, 5.0, 7.0], 8.0, 0, 3, + {'result': 0, 'ret_val': TARDIS_ERROR_BOUNDS_ERROR}), + + ([2.0, 4.0, 6.0, 7.0], 4.0, 0, 3, + {'result': 0, 'ret_val': TARDIS_ERROR_OK}) + ] +) +def test_binary_search(x, x_insert, imin, imax, expected_params): + obtained_result = 0 + + obtained_tardis_error = formal_integral.binary_search( + x, x_insert, imin, imax, obtained_result) + + assert obtained_result.value == expected_params['result'] + assert obtained_tardis_error == expected_params['ret_val'] + + +@pytest.mark.parametrize( + ['mu', 'r', 'inv_t_exp', 'expected'], + [(0.3, 7.5e14, 1 / 5.2e7, 0.9998556693818854), + (-.3, 8.1e14, 1 / 2.6e7, 1.0003117541351274)] +) +def test_rpacket_doppler_factor(mu, r, inv_t_exp, expected, packet, model): + # Set the params from test cases here + # TODO: add relativity tests + time_explosion = 1/inv_t_exp + + # Perform any other setups just before this, they can be additional calls + # to other methods or introduction of some temporary variables + + obtained = r_packet.get_doppler_factor(r, mu, time_explosion) + + # Perform required assertions + assert_almost_equal(obtained, expected) + + +@pytest.mark.parametrize( + ['packet_params', 'expected_params'], + [({'mu': 0.3, 'r': 7.5e14}, + {'d_boundary': 259376919351035.88}), + + ({'mu': -.3, 'r': 7.5e13}, + {'d_boundary': -664987228972291.5}), + + ({'mu': -.3, 'r': 7.5e14}, + {'d_boundary': 709376919351035.9})] +) +def test_compute_distance2boundary(packet_params, expected_params, packet, model): + mu = packet_params['mu'] + r = packet_params['r'] + + d_boundary = r_packet.calculate_distance_boundary( + r, mu, model.r_inner[0], model.r_outer[0]) + + assert_almost_equal(d_boundary[0], expected_params['d_boundary']) +# +# +# TODO: split this into two tests - one to assert errors and other for d_line +@pytest.mark.parametrize( + ['packet_params', 'expected_params'], + [({'nu_line': 0.1, 'next_line_id': 0, 'last_line': 1}, + {'tardis_error': None, 'd_line': 1e+99}), + + ({'nu_line': 0.2, 'next_line_id': 1, 'last_line': 0}, + {'tardis_error': None, 'd_line': 7.792353908000001e+17}), + + ({'nu_line': 0.5, 'next_line_id': 1, 'last_line': 0}, + {'tardis_error': r_packet.MonteCarloException, 'd_line': 0.0}), + + ({'nu_line': 0.6, 'next_line_id': 0, 'last_line': 0}, + {'tardis_error': r_packet.MonteCarloException, 'd_line': 0.0})] +) +def test_compute_distance2line(packet_params, expected_params, packet, model): + packet = r_packet.RPacket(packet.r, packet.mu, packet.nu, packet.energy) + nu_line = packet_params['nu_line'] + # packet.next_line_id = packet_params['next_line_id'] + # packet.last_line = packet_params['last_line'] + + time_explosion = model.time_explosion + + doppler_factor = r_packet.get_doppler_factor(packet.r, + packet.mu, + time_explosion) + comov_nu = packet.nu * doppler_factor + + d_line = 0 + obtained_tardis_error = None + try: + d_line = r_packet.calculate_distance_line(packet, + comov_nu, + nu_line, + time_explosion) + except r_packet.MonteCarloException: + obtained_tardis_error = r_packet.MonteCarloException + + assert_almost_equal(d_line, expected_params['d_line']) + assert obtained_tardis_error == expected_params['tardis_error'] + + +# @pytest.mark.parametrize( +# ['packet_params', 'expected_params'], +# [({'virtual_packet': 0}, +# {'chi_cont': 6.652486e-16, 'd_cont': 4.359272608766106e+28}), +# +# ({'virtual_packet': 1}, +# {'chi_cont': 6.652486e-16, 'd_cont': 1e+99})] +# ) +# def test_compute_distance2continuum(clib, packet_params, expected_params, packet, model): +# packet.virtual_packet = packet_params['virtual_packet'] +# +# clib.compute_distance2continuum(byref(packet), byref(model)) +# +# assert_almost_equal(packet.chi_cont, expected_params['chi_cont']) +# assert_almost_equal(packet.d_cont, expected_params['d_cont']) + + +# @pytest.mark.parametrize('full_relativity', [1, 0]) +@pytest.mark.parametrize( + ['packet_params', 'expected_params'], + [({'nu': 0.4, 'mu': 0.3, 'energy': 0.9, 'r': 7.5e14}, + {'mu': 0.3120599529139568, 'r': 753060422542573.9, + 'j': 8998701024436.969, 'nubar': 3598960894542.354}), + + ({'nu': 0.6, 'mu': -.5, 'energy': 0.5, 'r': 8.1e14}, + {'mu': -.4906548373534084, 'r': 805046582503149.2, + 'j': 5001298975563.031, 'nubar': 3001558973156.1387})] +) +def test_move_packet(packet_params, expected_params, + packet, model): + # TODO: test relativity + packet.nu = packet_params['nu'] + packet.mu = packet_params['mu'] + packet.energy = packet_params['energy'] + packet.r = packet_params['r'] + # model.full_relativity = full_relativity + + doppler_factor = r_packet.get_doppler_factor(packet.r, + packet.mu, + model.time_explosion) + r_packet.move_r_packet( + packet, distance, time_explosion, numba_estimator) + + assert_almost_equal(packet.mu, expected_params['mu']) + assert_almost_equal(packet.r, expected_params['r']) + + expected_j = expected_params['j'] + expected_nubar = expected_params['nubar'] + # if full_relativity: + # expected_j *= doppler_factor + # expected_nubar *= doppler_factor + + assert_allclose(model.js[packet.current_shell_id], + expected_j, rtol=5e-7) + assert_allclose(model.nubars[packet.current_shell_id], + expected_nubar, rtol=5e-7) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['packet_params', 'j_blue_idx', 'expected'], + [({'nu': 0.30, 'energy': 0.30}, 0, 1.0), + ({'nu': 0.20, 'energy': 1.e5}, 0, 5e5), + ({'nu': 2e15, 'energy': 0.50}, 1, 2.5e-16), + ({'nu': 0.40, 'energy': 1e-7}, 1, 2.5e-7)], +) +def test_increment_j_blue_estimator_full_relativity(clib, packet_params, + j_blue_idx, expected, + packet, model): + packet.nu = packet_params['nu'] + packet.energy = packet_params['energy'] + model.full_relativity = True + + clib.increment_j_blue_estimator(byref(packet), byref(model), + c_double(packet.d_line), + c_int64(j_blue_idx)) + + assert_almost_equal(model.line_lists_j_blues[j_blue_idx], expected) + + +@pytest.mark.parametrize( + ['packet_params', 'j_blue_idx', 'expected'], + [({'nu': 0.1, 'mu': 0.3, 'r': 7.5e14}, 0, 8.998643292289723), + ({'nu': 0.2, 'mu': -.3, 'r': 7.7e14}, 0, 4.499971133976377), + ({'nu': 0.5, 'mu': 0.5, 'r': 7.9e14}, 1, 0.719988453650551), + ({'nu': 0.6, 'mu': -.5, 'r': 8.1e14}, 1, 0.499990378058792)] +) +def test_increment_j_blue_estimator(clib, packet_params, j_blue_idx, expected, packet, model): + packet.nu = packet_params['nu'] + packet.mu = packet_params['mu'] + packet.r = packet_params['r'] + + clib.compute_distance2line(byref(packet), byref(model)) + clib.move_packet(byref(packet), byref(model), c_double(1.e13)) + clib.increment_j_blue_estimator(byref(packet), byref(model), + c_double(packet.d_line), c_int64(j_blue_idx)) + + assert_almost_equal(model.line_lists_j_blues[j_blue_idx], expected) + + +@pytest.mark.parametrize( + ['packet_params', 'expected_params'], + [({'virtual_packet': 0, 'current_shell_id': 0, 'next_shell_id': 1}, + {'status': TARDIS_PACKET_STATUS_IN_PROCESS, 'current_shell_id': 1}), + + ({'virtual_packet': 1, 'current_shell_id': 1, 'next_shell_id': 1}, + {'status': TARDIS_PACKET_STATUS_EMITTED, 'current_shell_id': 1, + 'tau_event': 29000000000000.008}), + + ({'virtual_packet': 1, 'current_shell_id': 0, 'next_shell_id': -1}, + {'status': TARDIS_PACKET_STATUS_REABSORBED, 'current_shell_id': 0, + 'tau_event': 29000000000000.008})] +) +def test_move_packet_across_shell_boundary(clib, packet_params, expected_params, + packet, model, mt_state): + packet.virtual_packet = packet_params['virtual_packet'] + packet.current_shell_id = packet_params['current_shell_id'] + packet.next_shell_id = packet_params['next_shell_id'] + + clib.move_packet_across_shell_boundary(byref(packet), byref(model), + c_double(1.e13), byref(mt_state)) + + if packet_params['virtual_packet'] == 1: + assert_almost_equal(packet.tau_event, expected_params['tau_event']) + assert packet.status == expected_params['status'] + assert packet.current_shell_id == expected_params['current_shell_id'] + + +@pytest.mark.parametrize( + ['packet_params', 'expected_params'], + [({'nu': 0.4, 'mu': 0.3, 'energy': 0.9, 'r': 7.5e14}, + {'nu': 0.39974659819356556, 'energy': 0.8994298459355226}), + + ({'nu': 0.6, 'mu': -.5, 'energy': 0.5, 'r': 8.1e14}, + {'nu': 0.5998422620533325, 'energy': 0.4998685517111104})] +) +def test_montecarlo_thomson_scatter(clib, packet_params, expected_params, packet, + model, mt_state): + packet.nu = packet_params['nu'] + packet.mu = packet_params['mu'] + packet.energy = packet_params['energy'] + packet.r = packet_params['r'] + + clib.montecarlo_thomson_scatter(byref(packet), byref(model), + c_double(1.e13), byref(mt_state)) + + assert_almost_equal(packet.nu, expected_params['nu']) + assert_almost_equal(packet.energy, expected_params['energy']) + + +@pytest.mark.parametrize( + ['packet_params', 'expected_params'], + # TODO: Add scientifically sound test cases. + [({'virtual_packet': 1, 'tau_event': 2.9e13, 'last_line': 0}, + {'tau_event': 2.9e13, 'next_line_id': 2}), + + ({'virtual_packet': 0, 'tau_event': 2.9e13, 'last_line': 0}, + {'tau_event': 2.9e13, 'next_line_id': 2}), + + ({'virtual_packet': 0, 'tau_event': 2.9e13, 'last_line': 0}, + {'tau_event': 2.9e13, 'next_line_id': 2}), + ] +) +def test_montecarlo_line_scatter(clib, packet_params, expected_params, packet, model, mt_state): + packet.virtual_packet = packet_params['virtual_packet'] + packet.tau_event = packet_params['tau_event'] + packet.last_line = packet_params['last_line'] + + clib.montecarlo_line_scatter(byref(packet), byref(model), + c_double(1.e13), byref(mt_state)) + + assert_almost_equal(packet.tau_event, expected_params['tau_event']) + assert_almost_equal(packet.next_line_id, expected_params['next_line_id']) + + +@pytest.mark.parametrize( + ['distances', 'expected'], + [({'boundary': 1.3e13, 'continuum': 1e14, 'line': 1e15}, + {'handler': 'move_packet_across_shell_boundary', 'distance': 1.3e13}), + + ({'boundary': 1.3e13, 'continuum': 1e14, 'line': 2.5e12}, + {'handler': 'montecarlo_line_scatter', 'distance': 2.5e12}), + + ({'boundary': 1.3e13, 'continuum': 1e11, 'line': 2.5e12}, + {'handler': 'montecarlo_thomson_scatter', 'distance': 1e11})] +) +def test_get_event_handler(clib, packet, model, mt_state, distances, expected): + d_cont_setter(distances['continuum'], model, packet) + d_line_setter(distances['line'], model, packet) + d_boundary_setter(distances['boundary'], model, packet) + obtained_distance = c_double() + + clib.get_event_handler.restype = c_void_p + obtained_handler = clib.get_event_handler(byref(packet), byref(model), + byref(obtained_distance), + byref(mt_state)) + + expected_handler = getattr(clib, expected['handler']) + expected_handler = cast(expected_handler, c_void_p).value + + assert_equal(obtained_handler, expected_handler) + assert_allclose(obtained_distance.value, expected['distance'], rtol=1e-10) + + +@pytest.mark.parametrize( + ['z_random', 'packet_params', 'expected'], + [(0.22443743797312765, + {'activation_level': 1, 'shell_id': 0}, 1), # Direct deactivation + + (0.78961460371187597, # next z_random = 0.818455414618 + {'activation_level': 1, 'shell_id': 0}, 0), # Upwards jump, then deactivation + + (0.22443743797312765, # next z_random = 0.545678896748 + {'activation_level': 2, 'shell_id': 0}, 1), # Downwards jump, then deactivation + + (0.765958602560605, # next z_random = 0.145914243888, 0.712382380384 + {'activation_level': 1, 'shell_id': 0}, 1), # Upwards jump, downwards jump, then deactivation + + (0.22443743797312765, + {'activation_level': 2, 'shell_id': 1}, 0)] # Direct deactivation +) +def test_macro_atom(clib, model_3lvlatom, packet, z_random, packet_params, get_rkstate, expected): + packet.macro_atom_activation_level = packet_params['activation_level'] + packet.current_shell_id = packet_params['shell_id'] + rkstate = get_rkstate(z_random) + + clib.macro_atom(byref(packet), byref(model_3lvlatom), byref(rkstate)) + obtained_line_id = model_3lvlatom.last_line_interaction_out_id[packet.id] + + assert_equal(obtained_line_id, expected) + + + +""" +Simple Tests: +---------------- +These test check very simple pices of code still work. +""" + +@pytest.mark.parametrize( + ['packet_params', 'line_idx', 'expected'], + [({'energy': 0.0}, 0, 0), + ({'energy': 1.0}, 1, 1), + ({'energy': 0.5}, 2, 1.5)] +) +@pytest.mark.skipif(True, reason="Needs rewrite to be relevant.") +def test_increment_Edotlu_estimator(clib, packet_params, line_idx, expected, packet, model): + packet.energy = packet_params['energy'] + + clib.increment_Edotlu_estimator(byref(packet), byref(model), c_int64(line_idx)) + + assert_almost_equal(model.line_lists_Edotlu[line_idx], expected) + + +""" +Difficult Tests: +---------------- +The tests written further are more complex than previous tests. They require +proper design procedure. They are not taken up yet and intended to be +completed together in future. +""" + + + +@pytest.mark.skipif(True, reason="Yet to be written.") +def test_montecarlo_one_packet(packet, model, mt_state): + pass + + +@pytest.mark.skipif(True, reason="Yet to be written.") +def test_montecarlo_one_packet_loop(packet, model, mt_state): + pass + + +@pytest.mark.skipif(True, reason="Yet to be written.") +def test_montecarlo_main_loop(packet, model, mt_state): + pass + + + +""" +Continuum Tests: +---------------- +The tests written further (till next block comment is encountered) are for the +methods related to continuum interactions. +""" + + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + 't_electron', [2500., 15000.] +) +def test_sample_nu_free_free(clib, t_electron, packet, model, mt_state_seeded, expected_ff_emissivity): + model.t_electrons[packet.current_shell_id] = t_electron + clib.sample_nu_free_free.restype = c_double + + nu_bins, expected_emissivity = expected_ff_emissivity(t_electron) + + nus = [] + for _ in range(int(1e5)): + nu = clib.sample_nu_free_free(byref(packet), byref(model), byref(mt_state_seeded)) + nus.append(nu) + + obtained_emissivity, _ = np.histogram(nus, density=True, bins=nu_bins) + + assert_allclose(obtained_emissivity, expected_emissivity, rtol=1e-10) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['packet_params', 't_electrons', 'chi_ff_factor', 'expected'], + [({'nu': 4.5e14, 'mu': 0.0, 'current_shell_id': 1}, 15000, 2.0, 1.6746639430359494e-44), + ({'nu': 3.0e15, 'mu': 0.0, 'current_shell_id': 0}, 5000, 3.0, 1.1111111111107644e-46), + ({'nu': 3.0e15, 'mu': 0.4, 'current_shell_id': 0}, 10000, 4.0, 1.5638286016098277e-46)] +) +def test_calculate_chi_ff(clib, packet, model, packet_params, t_electrons, chi_ff_factor, expected): + packet.mu = packet_params['mu'] + packet.nu = packet_params['nu'] + packet.current_shell_id = packet_params['current_shell_id'] + packet.r = 1.04e17 + + model.t_electrons[packet_params['current_shell_id']] = t_electrons + model.chi_ff_factor[packet_params['current_shell_id']] = chi_ff_factor + + clib.calculate_chi_ff(byref(packet), byref(model)) + obtained = packet.chi_ff + + assert_equal(obtained, expected) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['continuum_status', 'z_random', 'packet_params', 'expected'], + [(CONTINUUM_OFF, 0.94183547596539363, + {'chi_c': 1.0, 'chi_th': 0.4, 'chi_bf': 0.5}, + 'montecarlo_thomson_scatter'), + + (CONTINUUM_ON, 0.22443743797312765, + {'chi_c': 1.0, 'chi_th': 0.4, 'chi_bf': 0.5}, + 'montecarlo_thomson_scatter'), + + (CONTINUUM_ON, 0.54510721066252377, + {'chi_c': 1.0, 'chi_th': 0.4, 'chi_bf': 0.5}, + 'montecarlo_bound_free_scatter'), + + (CONTINUUM_ON, 0.94183547596539363, + {'chi_c': 1.0, 'chi_th': 0.4, 'chi_bf': 0.5}, + 'montecarlo_free_free_scatter'), + + (CONTINUUM_ON, 0.22443743797312765, + {'chi_c': 1e2, 'chi_th': 1e1, 'chi_bf': 2e1}, + 'montecarlo_bound_free_scatter')] +) +def test_montecarlo_continuum_event_handler(clib, continuum_status, expected, z_random, + packet_params, packet, model, get_rkstate): + packet.chi_cont = packet_params['chi_c'] + packet.chi_th = packet_params['chi_th'] + packet.chi_bf = packet_params['chi_bf'] + model.cont_status = continuum_status + + rkstate = get_rkstate(z_random) + + clib.montecarlo_continuum_event_handler.restype = c_void_p + obtained = clib.montecarlo_continuum_event_handler(byref(packet), + byref(model), byref(rkstate)) + expected = getattr(clib, expected) + expected = cast(expected, c_void_p).value + + assert_equal(obtained, expected) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['nu', 'continuum_id', 'expected', 'bf_treatment'], + [(4.40e14, 1, 0.00, BoundFreeTreatment.LIN_INTERPOLATION), + (3.25e14, 1, 0.75, BoundFreeTreatment.LIN_INTERPOLATION), + (4.03e14, 0, 0.97, BoundFreeTreatment.LIN_INTERPOLATION), + (4.10e14 + 1e-1, 0, 0.90, BoundFreeTreatment.LIN_INTERPOLATION), + pytest.param(4.1e14, 0, 0.90, BoundFreeTreatment.LIN_INTERPOLATION, + marks=pytest.mark.xfail), + (6.50e14, 0, 0.23304506144742834, BoundFreeTreatment.HYDROGENIC), + (3.40e14, 2, 1.1170364339507428, BoundFreeTreatment.HYDROGENIC)] +) +def test_bf_cross_section(clib, nu, continuum_id, model_w_edges, expected, bf_treatment): + model_w_edges.bf_treatment = bf_treatment.value + + clib.bf_cross_section.restype = c_double + obtained = clib.bf_cross_section(byref(model_w_edges), continuum_id, c_double(nu)) + + assert_almost_equal(obtained, expected) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['packet_params', 'expected'], + [({'nu': 4.13e14, 'mu': 0.0, 'current_shell_id': 1}, + [3.2882087455641473, 0.0, 0.0]), + + ({'nu': 3.27e14, 'mu': 0.0, 'current_shell_id': 0}, + [0.0, 1.3992114634681028, 5.702548202131454]), + + ({'nu': 3.27e14, 'mu': -0.4, 'current_shell_id': 0}, + [0.0, 1.2670858, 5.4446587])] +) +def test_calculate_chi_bf(clib, packet_params, expected, packet, model_w_edges): + model_w_edges.l_pop = (c_double * 6)(*range(1, 7)) + model_w_edges.l_pop_r = (c_double * 6)(*np.linspace(0.1, 0.6, 6)) + model_w_edges.t_electrons[packet_params['current_shell_id']] = 1e4 + + packet.mu = packet_params['mu'] + packet.nu = packet_params['nu'] + packet.r = 1.04e17 + packet.current_shell_id = packet_params['current_shell_id'] + packet.chi_bf_tmp_partial = (c_double * model_w_edges.no_of_edges)() + + clib.calculate_chi_bf(byref(packet), byref(model_w_edges)) + + obtained_chi_bf_tmp = np.ctypeslib.as_array(packet.chi_bf_tmp_partial, shape=(model_w_edges.no_of_edges,)) + expected_chi_bf_tmp = np.array(expected) + expected_chi_bf = expected_chi_bf_tmp[expected_chi_bf_tmp > 0][-1] + + assert_almost_equal(obtained_chi_bf_tmp, expected_chi_bf_tmp) + assert_almost_equal(packet.chi_bf, expected_chi_bf) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['comov_energy', 'distance', 'chi_ff', 'no_of_updates', 'expected'], + [(1.3, 1.3e14, 3e-12, 1, 507.), + (0.9, 0.7e15, 2e-12, 10, 1.260e4), + (0.8, 0.8e13, 1.5e-12, 35, 336.)] +) +def test_increment_continuum_estimators_ff_heating_estimator(clib, packet, model_w_edges, comov_energy, distance, + chi_ff, no_of_updates, expected): + packet.chi_ff = chi_ff + + for _ in range(no_of_updates): + clib.increment_continuum_estimators(byref(packet), byref(model_w_edges), c_double(distance), + c_double(0), c_double(comov_energy)) + obtained = model_w_edges.ff_heating_estimator[packet.current_shell_id] + + assert_almost_equal(obtained, expected) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['comov_nus', 'expected'], + [([4.05e14, 4.17e14, 3.3e14, 3.2e14, 2.9e14], [2, 2, 3]), + ([4.15e15, 3.25e14, 3.3e14, 2.85e14, 2.9e14], [0, 2, 4])] +) +def test_increment_continuum_estimators_photo_ion_estimator_statistics(clib, packet, model_w_edges, comov_nus, expected): + for comov_nu in comov_nus: + clib.increment_continuum_estimators(byref(packet), byref(model_w_edges), c_double(1e13), + c_double(comov_nu), c_double(1.0)) + + no_of_edges = model_w_edges.no_of_edges + no_of_shells = model_w_edges.no_of_shells + + obtained = np.ctypeslib.as_array(model_w_edges.photo_ion_estimator_statistics, + shape=(no_of_edges * no_of_shells,)) + obtained = np.reshape(obtained, newshape=(no_of_shells, no_of_edges), order='F') + obtained = obtained[packet.current_shell_id] + expected = np.array(expected) + + assert_array_equal(obtained, expected) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['comov_energy', 'distance', 'comov_nus', 'expected'], + [(1.3, 1.3e14, [4.05e14, 2.65e14], + {"photo_ion": [0.39641975308641975, 0., 0.], + "stim_recomb": [0.056757061269242064, 0., 0.], + "bf_heating": [1.9820987654321076e12, 0., 0.], + "stim_recomb_cooling": [283784812699.75476, 0., 0.]}), + + (0.9, 0.7e15, [3.25e14, 2.85e14], + {"photo_ion": [0., 1.4538461538461538, 7.315141700404858], + "stim_recomb": [0., 0.3055802, 1.7292954], + "bf_heating": [0., 36346153846153.82, 156760323886639.69], + "stim_recomb_cooling": [0., 7639505724285.9746, 33907776077426.875]})] +) +def test_increment_continuum_estimators_bf_estimators(clib, packet, model_w_edges, comov_energy, + distance, comov_nus, expected): + for comov_nu in comov_nus: + clib.increment_continuum_estimators(byref(packet), byref(model_w_edges), c_double(distance), + c_double(comov_nu), c_double(comov_energy)) + + no_of_edges = model_w_edges.no_of_edges + no_of_shells = model_w_edges.no_of_shells + + for estim_name, expected_value in expected.items(): + obtained = np.ctypeslib.as_array(getattr(model_w_edges, estim_name + "_estimator"), + shape=(no_of_edges * no_of_shells,)) + obtained = np.reshape(obtained, newshape=(no_of_shells, no_of_edges), order='F') + obtained = obtained[packet.current_shell_id] + + assert_almost_equal(obtained, np.array(expected_value)) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['packet_params', 'expected_params'], + [({'nu_comov': 1.1e16, 'mu': 0.0, 'r': 1.4e14}, + {'next_line_id': 5, 'last_line': True, 'nu': 1.1e16, 'type_id': 3}), + + ({'nu_comov': 1.3e16, 'mu': 0.3, 'r': 7.5e14}, + {'next_line_id': 0, 'last_line': False, 'nu': 1.30018766e+16, 'type_id': 3}), + + ({'nu_comov': 1.24e16, 'mu': -0.3, 'r': 7.5e14}, + {'next_line_id': 2, 'last_line': False, 'nu': 1.23982106e+16, 'type_id': 4})] +) +def test_continuum_emission(clib, packet, model, mock_sample_nu, packet_params, expected_params, mt_state): + packet.nu = packet_params['nu_comov'] # Is returned by mock function mock_sample_nu + packet.mu = packet_params['mu'] + packet.r = packet_params['r'] + expected_interaction_out_type = expected_params['type_id'] + + clib.continuum_emission(byref(packet), byref(model), byref(mt_state), + mock_sample_nu, expected_interaction_out_type) + + obtained_next_line_id = packet.next_line_id + obtained_last_interaction_out_type = model.last_interaction_out_type[0] + + assert_equal(obtained_next_line_id, expected_params['next_line_id']) + assert_equal(packet.last_line, expected_params['last_line']) + assert_equal(expected_interaction_out_type, obtained_last_interaction_out_type) + assert_allclose(packet.nu, expected_params['nu'], rtol=1e-7) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['packet_params', 'expected'], + [({'next_line_id': 3, 'last_line': 0}, 1), + ({'next_line_id': 5, 'last_line': 1}, 0), + ({'next_line_id': 2, 'last_line': 0}, 0), + ({'next_line_id': 1, 'last_line': 0}, 1)] +) +def test_test_for_close_line(clib, packet, model, packet_params, expected): + packet.nu_line = model.line_list_nu[packet_params['next_line_id'] - 1] + packet.next_line_id = packet_params['next_line_id'] + + clib.test_for_close_line(byref(packet), byref(model)) + + assert_equal(expected, packet.close_line) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['packet_params', 'z_random', 'expected'], + [({'current_continuum_id': 1, 'chi_bf_tmp_partial': [0.0, 0.23e13, 1.0e13]}, + 0.22443743797312765, 1), + + ({'current_continuum_id': 1, 'chi_bf_tmp_partial': [0.0, 0.23e10, 1.0e10]}, + 0.78961460371187597, 2), + + ({'current_continuum_id': 0, 'chi_bf_tmp_partial': [0.2e5, 0.5e5, 0.6e5, 1.0e5, 1.0e5]}, + 0.78961460371187597, 3)] +) +def test_montecarlo_bound_free_scatter_continuum_selection(clib, packet, model_3lvlatom, packet_params, + get_rkstate, z_random, expected): + rkstate = get_rkstate(z_random) + packet.current_continuum_id = packet_params['current_continuum_id'] + + chi_bf_tmp = packet_params['chi_bf_tmp_partial'] + packet.chi_bf_tmp_partial = (c_double * len(chi_bf_tmp))(*chi_bf_tmp) + packet.chi_bf = chi_bf_tmp[-1] + model_3lvlatom.no_of_edges = len(chi_bf_tmp) + + clib.montecarlo_bound_free_scatter(byref(packet), byref(model_3lvlatom), + c_double(1.e13), byref(rkstate)) + + assert_equal(packet.current_continuum_id, expected) + assert_equal(model_3lvlatom.last_line_interaction_in_id[packet.id], expected) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['mu', 'r', 'inv_t_exp', 'full_relativity'], + [(0.8, 7.5e14, 1 / 5.2e5, 1), + (-0.7, 7.5e14, 1 / 5.2e5, 1), + (0.3, 7.5e14, 1 / 2.2e5, 1), + (0.0, 7.5e14, 1 / 2.2e5, 1), + (-0.7, 7.5e14, 1 / 5.2e5, 0)] +) +def test_frame_transformations(clib, packet, model, mu, r, + inv_t_exp, full_relativity): + packet.r = r + packet.mu = mu + model.inverse_time_explosion = inv_t_exp + model.full_relativity = full_relativity + clib.rpacket_doppler_factor.restype = c_double + clib.rpacket_inverse_doppler_factor.restype = c_double + + inverse_doppler_factor = clib.rpacket_inverse_doppler_factor(byref(packet), byref(model)) + clib.angle_aberration_CMF_to_LF(byref(packet), byref(model)) + + doppler_factor = clib.rpacket_doppler_factor(byref(packet), byref(model)) + + assert_almost_equal(doppler_factor * inverse_doppler_factor, 1.0) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + ['mu', 'r', 'inv_t_exp'], + [(0.8, 7.5e14, 1 / 5.2e5), + (-0.7, 7.5e14, 1 / 5.2e5), + (0.3, 7.5e14, 1 / 2.2e5), + (0.0, 7.5e14, 1 / 2.2e5), + (-0.7, 7.5e14, 1 / 5.2e5)] +) +def test_angle_transformation_invariance(clib, packet, model, + mu, r, inv_t_exp): + packet.r = r + packet.mu = mu + model.inverse_time_explosion = inv_t_exp + model.full_relativity = 1 + clib.angle_aberration_LF_to_CMF.restype = c_double + + clib.angle_aberration_CMF_to_LF(byref(packet), byref(model)) + mu_obtained = clib.angle_aberration_LF_to_CMF( + byref(packet), byref(model), c_double(packet.mu)) + + assert_almost_equal(mu_obtained, mu) + + +@pytest.mark.continuumtest +@pytest.mark.parametrize( + 'full_relativity', + [1, 0] +) +@pytest.mark.parametrize( + ['mu', 'r', 't_exp', 'nu', 'nu_line'], + [(0.8, 7.5e14, 5.2e5, 1.0e15, 9.4e14), + (0.0, 6.3e14, 2.2e5, 6.0e12, 5.8e12), + (1.0, 9.0e14, 2.2e5, 4.0e8, 3.4e8), + (0.9, 9.0e14, 0.5e5, 1.0e15, 4.5e14), + (-0.7, 7.5e14, 5.2e5, 1.0e15, 9.8e14), + (-1.0, 6.3e14, 2.2e5, 6.0e12, 6.55e12)] +) +def test_compute_distance2line_relativistic(clib, mu, r, t_exp, nu, nu_line, + full_relativity, packet, model): + packet.r = r + packet.mu = mu + packet.nu = nu + packet.nu_line = nu_line + model.inverse_time_explosion = 1 / t_exp + model.time_explosion = t_exp + model.full_relativity = full_relativity + + clib.rpacket_doppler_factor.restype = c_double + + clib.compute_distance2line(byref(packet), byref(model)) + clib.move_packet(byref(packet), byref(model), c_double(packet.d_line)) + + doppler_factor = clib.rpacket_doppler_factor(byref(packet), byref(model)) + comov_nu = packet.nu * doppler_factor + + assert_allclose(comov_nu, nu_line, rtol=1e-14) + + +@pytest.mark.continuumtest +@pytest.mark.skipif(True, reason="Yet to be written.") +def test_montecarlo_free_free_scatter(packet, model, mt_state): + pass From 5300b7aea22d12414f40353c0a19148bda1070e2 Mon Sep 17 00:00:00 2001 From: Arjun Savel Date: Thu, 2 Jul 2020 13:43:25 -0700 Subject: [PATCH 03/21] Add doppler factor tests --- .../montecarlo/montecarlo_numba/r_packet.py | 9 ++ tardis/montecarlo/tests/test_montecarlo.py | 109 +++++++++++++++++- 2 files changed, 115 insertions(+), 3 deletions(-) diff --git a/tardis/montecarlo/montecarlo_numba/r_packet.py b/tardis/montecarlo/montecarlo_numba/r_packet.py index 70fcb0bed50..9ebd4716267 100644 --- a/tardis/montecarlo/montecarlo_numba/r_packet.py +++ b/tardis/montecarlo/montecarlo_numba/r_packet.py @@ -12,6 +12,9 @@ class MonteCarloException(ValueError): pass +class SuperluminalError(ValueError): + pass + CLOSE_LINE_THRESHOLD = 1e-7 C_SPEED_OF_LIGHT = const.c.to('cm/s').value MISS_DISTANCE = 1e99 @@ -129,6 +132,8 @@ def calculate_tau_electron(electron_density, distance): @njit(**njit_dict) def get_doppler_factor(r, mu, time_explosion): beta = r / (time_explosion * C_SPEED_OF_LIGHT) + if beta > 1.: # TODO: check speed penalty + raise SuperluminalError if not montecarlo_configuration.full_relativity: return get_doppler_factor_partial_relativity(mu, beta) else: @@ -140,12 +145,15 @@ def get_doppler_factor_partial_relativity(mu, beta): @njit(**njit_dict) def get_doppler_factor_full_relativity(mu, beta): + # todo: check numpy vs other sqrt return (1.0 - mu * beta) / np.sqrt(1 - beta * beta) @njit(**njit_dict) def get_inverse_doppler_factor(r, mu, time_explosion): beta = (r / time_explosion) / C_SPEED_OF_LIGHT + if beta > 1.: + raise SuperluminalError if not montecarlo_configuration.full_relativity: return get_inverse_doppler_factor_partial_relativity(mu, beta) else: @@ -157,6 +165,7 @@ def get_inverse_doppler_factor_partial_relativity(mu, beta): @njit(**njit_dict) def get_inverse_doppler_factor_full_relativity(mu, beta): + # TODO: check numpy vs. other sqrt return (1.0 + mu * beta) / np.sqrt(1 - beta * beta) @njit(**njit_dict) diff --git a/tardis/montecarlo/tests/test_montecarlo.py b/tardis/montecarlo/tests/test_montecarlo.py index 86bc15d72d7..e0db0f6a752 100644 --- a/tardis/montecarlo/tests/test_montecarlo.py +++ b/tardis/montecarlo/tests/test_montecarlo.py @@ -50,7 +50,7 @@ import pandas as pd import tardis.montecarlo.formal_integral as formal_integral import tardis.montecarlo.montecarlo_numba.r_packet as r_packet - +import tardis.montecarlo.montecarlo_configuration as mc from ctypes import ( @@ -317,18 +317,121 @@ def test_binary_search(x, x_insert, imin, imax, expected_params): @pytest.mark.parametrize( ['mu', 'r', 'inv_t_exp', 'expected'], [(0.3, 7.5e14, 1 / 5.2e7, 0.9998556693818854), - (-.3, 8.1e14, 1 / 2.6e7, 1.0003117541351274)] + (-.3, 0, 1 / 2.6e7, 1.0), + (0, 1, 1 / 2.6e7, 1.0)] +) +def test_get_doppler_factor(mu, r, inv_t_exp, expected, packet, model): + # Set the params from test cases here + # TODO: add relativity tests + time_explosion = 1/inv_t_exp + + # Perform any other setups just before this, they can be additional calls + # to other methods or introduction of some temporary variables + + obtained = r_packet.get_doppler_factor(r, mu, time_explosion) + + + # Perform required assertions + assert_almost_equal(obtained, expected) + + +@pytest.mark.parametrize( + ['mu', 'r', 'inv_t_exp'], + [(-0.3, 5, 1e10)] +) +def test_unphysical_doppler_factor(mu, r, inv_t_exp, expected, packet, model): + # Set the params from test cases here + # TODO: add relativity tests + time_explosion = 1 / inv_t_exp + + # Perform any other setups just before this, they can be additional calls + # to other methods or introduction of some temporary variables + with pytest.raises(r_packet.SuperluminalError): + obtained = r_packet.get_doppler_factor(r, mu, time_explosion) + + +@pytest.mark.parametrize( + ['mu', 'r', 'inv_t_exp', 'expected'], + [(0.3, 7.5e14, 1 / 5.2e7, 1/0.9998556693818854), + (-.3, 0, 1 / 2.6e7, 1.0), + (0, 1, 1 / 2.6e7, 1.0)] +) +def test_get_inverse_doppler_factor(mu, r, inv_t_exp, expected, packet, model): + # Set the params from test cases here + # TODO: add relativity tests + time_explosion = 1/inv_t_exp + + # Perform any other setups just before this, they can be additional calls + # to other methods or introduction of some temporary variables + + obtained = r_packet.get_inverse_doppler_factor(r, mu, time_explosion) + + # Perform required assertions + assert_almost_equal(obtained, expected) + +@pytest.mark.parametrize( + ['mu', 'r', 'inv_t_exp'], + [(-0.3, 5, 1e10)] ) -def test_rpacket_doppler_factor(mu, r, inv_t_exp, expected, packet, model): +def test_unphysical_inverse_doppler_factor(mu, r, inv_t_exp, packet, model): # Set the params from test cases here # TODO: add relativity tests time_explosion = 1/inv_t_exp + # Perform any other setups just before this, they can be additional calls + # to other methods or introduction of some temporary variables + with pytest.raises(r_packet.SuperluminalError): + obtained = r_packet.get_inverse_doppler_factor(r, mu, time_explosion) + + +@pytest.mark.parametrize( + ['mu', 'r', 'inv_t_exp', 'expected'], + [(-0.3, 10000000, 0.001, 1.0000001000692842), + (-.3, 0, 1 / 2.6e7, 1.0), + (0, 1, 1 / 2.6e7, 1.0)] +) +def test_get_doppler_factor_full_relativity(mu, + r, + inv_t_exp, + expected, + packet, + model): + # Set the params from test cases here + # TODO: add relativity tests + mc.full_relativity = True + time_explosion = 1/inv_t_exp + # Perform any other setups just before this, they can be additional calls # to other methods or introduction of some temporary variables obtained = r_packet.get_doppler_factor(r, mu, time_explosion) + mc.full_relativity = False + # Perform required assertions + assert_almost_equal(obtained, expected) + + +@pytest.mark.parametrize( + ['mu', 'r', 'inv_t_exp', 'expected'], + [(-0.3, 10000000, 0.001, 0.999999899930827), + (-.3, 0, 1 / 2.6e7, 1.0), + (0, 1, 1 / 2.6e7, 1.0)] +) +def test_get_inverse_doppler_factor_full_relativity(mu, + r, + inv_t_exp, + expected, + packet, + model): + # Set the params from test cases here + # TODO: add relativity tests + mc.full_relativity = True + time_explosion = 1/inv_t_exp + + # Perform any other setups just before this, they can be additional calls + # to other methods or introduction of some temporary variables + obtained = r_packet.get_inverse_doppler_factor(r, mu, time_explosion) + mc.full_relativity = False # Perform required assertions assert_almost_equal(obtained, expected) From 6155e15c82aeeec19b159fef2512e38a8623a396 Mon Sep 17 00:00:00 2001 From: Arjun Savel Date: Sat, 4 Jul 2020 20:28:22 -0700 Subject: [PATCH 04/21] Add mu, LF_TO_CMF tests --- tardis/montecarlo/tests/test_montecarlo.py | 38 ++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/tardis/montecarlo/tests/test_montecarlo.py b/tardis/montecarlo/tests/test_montecarlo.py index e0db0f6a752..f0158c0edfe 100644 --- a/tardis/montecarlo/tests/test_montecarlo.py +++ b/tardis/montecarlo/tests/test_montecarlo.py @@ -51,6 +51,9 @@ import tardis.montecarlo.formal_integral as formal_integral import tardis.montecarlo.montecarlo_numba.r_packet as r_packet import tardis.montecarlo.montecarlo_configuration as mc +from tardis import constants as const +C_SPEED_OF_LIGHT = const.c.to('cm/s').value + from ctypes import ( @@ -435,6 +438,41 @@ def test_get_inverse_doppler_factor_full_relativity(mu, # Perform required assertions assert_almost_equal(obtained, expected) +def test_get_random_mu_different_output(): + """ + Ensure that different calls results + """ + output1 = r_packet.get_random_mu() + output2 = r_packet.get_random_mu() + assert output1 != output2 + +def test_get_random_mu_different_output(): + """ + Ensure that different calls results + """ + output1 = r_packet.get_random_mu() + output2 = r_packet.get_random_mu() + assert output1 != output2 + + +@pytest.mark.parametrize( + ['mu', 'r', 'time_explosion'], + [(1, C_SPEED_OF_LIGHT, 1), + (1, 1, C_SPEED_OF_LIGHT)] +) +def test_angle_ab_LF_to_CMF_diverge(mu, r, time_explosion, packet): + """ + This equation should diverage as costheta --> 1 and beta --> 1 + """ + C_SPEED_OF_LIGHT + packet.r = r + packet.mu = mu + with pytest.raises(ZeroDivisionError): + obtained = r_packet.angle_aberration_LF_to_CMF( + packet, + time_explosion, + mu) + @pytest.mark.parametrize( ['packet_params', 'expected_params'], From 415dd26846b04d8762a33b224bd4ba6aec946de5 Mon Sep 17 00:00:00 2001 From: Arjun Savel Date: Mon, 13 Jul 2020 21:04:04 -0700 Subject: [PATCH 05/21] Change r_inner, r_outer names; reference plasma correctly in integral test --- tardis/montecarlo/formal_integral.py | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/tardis/montecarlo/formal_integral.py b/tardis/montecarlo/formal_integral.py index 2788f2309c0..3f793df58e5 100644 --- a/tardis/montecarlo/formal_integral.py +++ b/tardis/montecarlo/formal_integral.py @@ -9,6 +9,8 @@ from tardis.montecarlo.montecarlo_numba import njit_dict +from tardis.montecarlo.montecarlo_numba.numba_interface \ + import numba_plasma_initialize # from tardis.montecarlo.montecarlo import formal_integral from tardis.montecarlo.spectrum import TARDISSpectrum @@ -33,7 +35,7 @@ class FormalIntegrator(object): def __init__(self, model, plasma, runner, points=1000): self.model = model - self.plasma = plasma + self.plasma = numba_plasma_initialize(plasma) self.runner = runner self.points = points @@ -81,8 +83,7 @@ def calculate_spectrum(self, frequency, points=None, frequency = frequency.to('Hz', u.spectral()) luminosity = u.Quantity( - formal_integral( - self, + self.formal_integral( frequency, N), 'erg' @@ -122,7 +123,7 @@ def make_source_function(self): plasma = self.plasma runner = self.runner atomic_data = self.plasma.atomic_data - macro_ref = atomic_data.macro_atom_references + macro_ref = plasma.macro_block_references macro_data = atomic_data.macro_atom_data no_lvls = len(atomic_data.levels) @@ -242,10 +243,6 @@ def formal_integral(self, nu, N): res = self.make_source_function() - - - - att_S_ul = res[0].flatten(order='F') Jred_lu = res[1].flatten(order='F') Jblue_lu = res[2].flatten(order='F') @@ -267,22 +264,20 @@ def _formal_integral(self, iT, inu, inu_size, att_S_ul, Jred_lu, Jblue_lu, N): # Initialize the output which is shared among threads L = np.zeros(inu_size) # global read-only values - size_line = self.model.no_of_lines # check - size_shell = self.model.no_of_shells_i # check + size_line = len(self.plasma.line_list_nu) + size_shell = self.model.no_of_shells # check size_tau = size_line * size_shell finished_nus = 0 - R_ph = self.runner.r_inner_i[0] - R_max = self.runner.r_outer_i[size_shell - 1] + R_ph = self.runner.r_inner[0] + R_max = self.runner.r_outer[size_shell - 1] pp = np.zeros(N) # check exp_tau = np.zeros(size_tau) # TODO: multiprocessing offset = 0 - i = 0 size_z = 0 idx_nu_start = 0 direction = 0 - first = 0 I_nu = np.zeros(N) shell_id = np.zeros(2 * self.runner.no_of_shells_i) # check # instantiate more variables here, maybe? From ac0156955ead70e7ddb478f5b187f2fc9138a360 Mon Sep 17 00:00:00 2001 From: Arjun Savel Date: Mon, 13 Jul 2020 21:05:04 -0700 Subject: [PATCH 06/21] Clean up montecarlo testing --- tardis/montecarlo/tests/test_montecarlo.py | 630 ++++++--------------- 1 file changed, 185 insertions(+), 445 deletions(-) diff --git a/tardis/montecarlo/tests/test_montecarlo.py b/tardis/montecarlo/tests/test_montecarlo.py index f0158c0edfe..c8ca25dde24 100644 --- a/tardis/montecarlo/tests/test_montecarlo.py +++ b/tardis/montecarlo/tests/test_montecarlo.py @@ -52,6 +52,8 @@ import tardis.montecarlo.montecarlo_numba.r_packet as r_packet import tardis.montecarlo.montecarlo_configuration as mc from tardis import constants as const +from tardis.montecarlo.montecarlo_numba import numba_interface +from tardis.montecarlo.montecarlo_numba import macro_atom C_SPEED_OF_LIGHT = const.c.to('cm/s').value @@ -263,7 +265,7 @@ def d_boundary_setter(d_boundary, model, packet): {'result': 0, 'ret_val': TARDIS_ERROR_BOUNDS_ERROR})] ) def test_reverse_binary_search(x, x_insert, imin, imax, expected_params): - x = (c_double * (imax - imin + 1))(*x) + # x = (c_double * (imax - imin + 1))(*x) obtained_result = 0 obtained_tardis_error = formal_integral.reverse_binary_search( @@ -285,7 +287,7 @@ def test_reverse_binary_search(x, x_insert, imin, imax, expected_params): {'result': 0, 'ret_val': TARDIS_ERROR_OK})] ) def test_line_search(nu, nu_insert, number_of_lines, expected_params): - nu = (c_double * number_of_lines)(*nu) + # nu = (c_double * number_of_lines)(*nu) nu_insert = c_double(nu_insert) obtained_result = 0 @@ -342,7 +344,7 @@ def test_get_doppler_factor(mu, r, inv_t_exp, expected, packet, model): ['mu', 'r', 'inv_t_exp'], [(-0.3, 5, 1e10)] ) -def test_unphysical_doppler_factor(mu, r, inv_t_exp, expected, packet, model): +def test_unphysical_doppler_factor(mu, r, inv_t_exp, packet, model): # Set the params from test cases here # TODO: add relativity tests time_explosion = 1 / inv_t_exp @@ -457,14 +459,13 @@ def test_get_random_mu_different_output(): @pytest.mark.parametrize( ['mu', 'r', 'time_explosion'], - [(1, C_SPEED_OF_LIGHT, 1), - (1, 1, C_SPEED_OF_LIGHT)] + [(1, C_SPEED_OF_LIGHT, 1)] ) def test_angle_ab_LF_to_CMF_diverge(mu, r, time_explosion, packet): """ This equation should diverage as costheta --> 1 and beta --> 1 """ - C_SPEED_OF_LIGHT + packet = r_packet.RPacket(packet.r, packet.mu, packet.nu, packet.energy) packet.r = r packet.mu = mu with pytest.raises(ZeroDivisionError): @@ -472,6 +473,111 @@ def test_angle_ab_LF_to_CMF_diverge(mu, r, time_explosion, packet): packet, time_explosion, mu) +@pytest.mark.parametrize( + ['mu', 'r', 'time_explosion'], + [(0.3, 1e7, 1e10), + (-0.3, 1e7, 1e11)] +) +def test_both_angle_aberrations(mu, r, time_explosion, packet): + """ + The angle aberration functions should be the functional inverse of one + another. + """ + packet = r_packet.RPacket(r, packet.mu, packet.nu, packet.energy) + packet.r = r + obtained_mu = r_packet.angle_aberration_LF_to_CMF( + packet, + time_explosion, + mu) + inverse_obtained_mu = r_packet.angle_aberration_CMF_to_LF( + packet, + time_explosion, + obtained_mu) + assert_almost_equal(inverse_obtained_mu, mu) + + +@pytest.mark.parametrize( + ['mu', 'r', 'time_explosion'], + [(0.3, 7.5e14, 5.2e5), + (-0.3, 7.5e14, 5.2e5)] +) +def test_both_angle_aberrations_inverse(mu, r, time_explosion, packet): + """ + The angle aberration functions should be the functional inverse of one + another. + """ + packet = r_packet.RPacket(r, packet.mu, packet.nu, packet.energy) + packet.r = r + obtained_mu = r_packet.angle_aberration_CMF_to_LF( + packet, + time_explosion, + mu) + inverse_obtained_mu = r_packet.angle_aberration_LF_to_CMF( + packet, + time_explosion, + obtained_mu) + assert_almost_equal(inverse_obtained_mu, mu) + + +@pytest.mark.parametrize( + ['current_shell_id', 'delta_shell', 'no_of_shells'], + [(132, 11, 132), + (132, 1, 133), + (132, 2, 133)] +) +def test_move_packet_across_shell_boundary_emitted(packet, current_shell_id, + delta_shell, + no_of_shells): + packet = r_packet.RPacket(packet.r, packet.mu, packet.nu, packet.energy) + packet.current_shell_id = current_shell_id + r_packet.move_packet_across_shell_boundary(packet, delta_shell, + no_of_shells) + assert packet.status == r_packet.PacketStatus.EMITTED + +@pytest.mark.parametrize( + ['current_shell_id', 'delta_shell', 'no_of_shells'], + [(132, -133, 132), + (132, -133, 133), + (132, -1e9, 133)] +) +def test_move_packet_across_shell_boundary_reabsorbed(packet, current_shell_id, + delta_shell, + no_of_shells): + packet = r_packet.RPacket(packet.r, packet.mu, packet.nu, packet.energy) + packet.current_shell_id = current_shell_id + r_packet.move_packet_across_shell_boundary(packet, delta_shell, + no_of_shells) + assert packet.status == r_packet.PacketStatus.REABSORBED + + +@pytest.mark.parametrize( + ['current_shell_id', 'delta_shell', 'no_of_shells'], + [(132, -1, 199), + (132, 0, 132), + (132, 20, 154)] +) +def test_move_packet_across_shell_boundary_increment(packet, current_shell_id, + delta_shell, + no_of_shells): + packet = r_packet.RPacket(packet.r, packet.mu, packet.nu, packet.energy) + packet.current_shell_id = current_shell_id + r_packet.move_packet_across_shell_boundary(packet, delta_shell, + no_of_shells) + assert packet.current_shell_id == current_shell_id + delta_shell + + +@pytest.mark.parametrize( + ['distance_trace', 'time_explosion', 'mu', 'r'], + [(0, 1, 0, 0), + (0, 1, 1, 0), + (0, 1, 0, 1)] +) +def test_packet_energy_limit_one(packet, distance_trace, time_explosion, mu, r): + initial_energy = packet.energy + packet = r_packet.RPacket(r, mu, packet.nu, packet.energy) + new_energy = r_packet.calc_packet_energy(packet, distance_trace, time_explosion) + assert new_energy == initial_energy + @pytest.mark.parametrize( @@ -595,74 +701,56 @@ def test_move_packet(packet_params, expected_params, expected_nubar, rtol=5e-7) -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['packet_params', 'j_blue_idx', 'expected'], - [({'nu': 0.30, 'energy': 0.30}, 0, 1.0), - ({'nu': 0.20, 'energy': 1.e5}, 0, 5e5), - ({'nu': 2e15, 'energy': 0.50}, 1, 2.5e-16), - ({'nu': 0.40, 'energy': 1e-7}, 1, 2.5e-7)], -) -def test_increment_j_blue_estimator_full_relativity(clib, packet_params, - j_blue_idx, expected, - packet, model): - packet.nu = packet_params['nu'] - packet.energy = packet_params['energy'] - model.full_relativity = True - - clib.increment_j_blue_estimator(byref(packet), byref(model), - c_double(packet.d_line), - c_int64(j_blue_idx)) - - assert_almost_equal(model.line_lists_j_blues[j_blue_idx], expected) - - -@pytest.mark.parametrize( - ['packet_params', 'j_blue_idx', 'expected'], - [({'nu': 0.1, 'mu': 0.3, 'r': 7.5e14}, 0, 8.998643292289723), - ({'nu': 0.2, 'mu': -.3, 'r': 7.7e14}, 0, 4.499971133976377), - ({'nu': 0.5, 'mu': 0.5, 'r': 7.9e14}, 1, 0.719988453650551), - ({'nu': 0.6, 'mu': -.5, 'r': 8.1e14}, 1, 0.499990378058792)] -) -def test_increment_j_blue_estimator(clib, packet_params, j_blue_idx, expected, packet, model): - packet.nu = packet_params['nu'] - packet.mu = packet_params['mu'] - packet.r = packet_params['r'] - - clib.compute_distance2line(byref(packet), byref(model)) - clib.move_packet(byref(packet), byref(model), c_double(1.e13)) - clib.increment_j_blue_estimator(byref(packet), byref(model), - c_double(packet.d_line), c_int64(j_blue_idx)) - - assert_almost_equal(model.line_lists_j_blues[j_blue_idx], expected) - - -@pytest.mark.parametrize( - ['packet_params', 'expected_params'], - [({'virtual_packet': 0, 'current_shell_id': 0, 'next_shell_id': 1}, - {'status': TARDIS_PACKET_STATUS_IN_PROCESS, 'current_shell_id': 1}), - - ({'virtual_packet': 1, 'current_shell_id': 1, 'next_shell_id': 1}, - {'status': TARDIS_PACKET_STATUS_EMITTED, 'current_shell_id': 1, - 'tau_event': 29000000000000.008}), - - ({'virtual_packet': 1, 'current_shell_id': 0, 'next_shell_id': -1}, - {'status': TARDIS_PACKET_STATUS_REABSORBED, 'current_shell_id': 0, - 'tau_event': 29000000000000.008})] -) -def test_move_packet_across_shell_boundary(clib, packet_params, expected_params, - packet, model, mt_state): - packet.virtual_packet = packet_params['virtual_packet'] - packet.current_shell_id = packet_params['current_shell_id'] - packet.next_shell_id = packet_params['next_shell_id'] - - clib.move_packet_across_shell_boundary(byref(packet), byref(model), - c_double(1.e13), byref(mt_state)) +# @pytest.mark.continuumtest +# @pytest.mark.parametrize( +# ['packet_params', 'j_blue_idx', 'expected'], +# [({'nu': 0.30, 'energy': 0.30}, 0, 1.0), +# ({'nu': 0.20, 'energy': 1.e5}, 0, 5e5), +# ({'nu': 2e15, 'energy': 0.50}, 1, 2.5e-16), +# ({'nu': 0.40, 'energy': 1e-7}, 1, 2.5e-7)], +# ) +# def test_increment_j_blue_estimator_full_relativity(packet_params, +# j_blue_idx, expected, +# packet, model): +# packet.nu = packet_params['nu'] +# packet.energy = packet_params['energy'] +# model.full_relativity = True +# +# r_packet.increment_j_blue_estimator(byref(packet), byref(model), +# c_double(packet.d_line), +# c_int64(j_blue_idx)) +# +# assert_almost_equal(model.line_lists_j_blues[j_blue_idx], expected) +# +# +# @pytest.mark.parametrize( +# ['packet_params', 'cur_line_id', 'expected'], +# [({'nu': 0.1, 'mu': 0.3, 'r': 7.5e14}, 0, 8.998643292289723), +# ({'nu': 0.2, 'mu': -.3, 'r': 7.7e14}, 0, 4.499971133976377), +# ({'nu': 0.5, 'mu': 0.5, 'r': 7.9e14}, 1, 0.719988453650551), +# ({'nu': 0.6, 'mu': -.5, 'r': 8.1e14}, 1, 0.499990378058792)] +# ) +# def test_increment_j_blue_estimator(packet_params, cur_line_id, expected, packet): +# +# numba_interface +# packet = r_packet.RPacket(packet_params['r'], +# packet_params['mu'], +# packet_params['nu'], +# packet.energy) +# +# r_packet.compute_distance2line(byref(packet), byref(model)) +# r_packet.move_r_packet(packet, +# distance, +# model.time_explosion, +# numba_estimator) +# r_packet.move_packet(byref(packet), byref(model), c_double(1.e13)) +# r_packet.update_line_estimators(estimators, r_packet, +# cur_line_id, +# model.distance_trace, +# model.time_explosion) +# +# assert_almost_equal(model.line_lists_j_blues[j_blue_idx], expected) - if packet_params['virtual_packet'] == 1: - assert_almost_equal(packet.tau_event, expected_params['tau_event']) - assert packet.status == expected_params['status'] - assert packet.current_shell_id == expected_params['current_shell_id'] @pytest.mark.parametrize( @@ -712,34 +800,6 @@ def test_montecarlo_line_scatter(clib, packet_params, expected_params, packet, m assert_almost_equal(packet.next_line_id, expected_params['next_line_id']) -@pytest.mark.parametrize( - ['distances', 'expected'], - [({'boundary': 1.3e13, 'continuum': 1e14, 'line': 1e15}, - {'handler': 'move_packet_across_shell_boundary', 'distance': 1.3e13}), - - ({'boundary': 1.3e13, 'continuum': 1e14, 'line': 2.5e12}, - {'handler': 'montecarlo_line_scatter', 'distance': 2.5e12}), - - ({'boundary': 1.3e13, 'continuum': 1e11, 'line': 2.5e12}, - {'handler': 'montecarlo_thomson_scatter', 'distance': 1e11})] -) -def test_get_event_handler(clib, packet, model, mt_state, distances, expected): - d_cont_setter(distances['continuum'], model, packet) - d_line_setter(distances['line'], model, packet) - d_boundary_setter(distances['boundary'], model, packet) - obtained_distance = c_double() - - clib.get_event_handler.restype = c_void_p - obtained_handler = clib.get_event_handler(byref(packet), byref(model), - byref(obtained_distance), - byref(mt_state)) - - expected_handler = getattr(clib, expected['handler']) - expected_handler = cast(expected_handler, c_void_p).value - - assert_equal(obtained_handler, expected_handler) - assert_allclose(obtained_distance.value, expected['distance'], rtol=1e-10) - @pytest.mark.parametrize( ['z_random', 'packet_params', 'expected'], @@ -758,12 +818,12 @@ def test_get_event_handler(clib, packet, model, mt_state, distances, expected): (0.22443743797312765, {'activation_level': 2, 'shell_id': 1}, 0)] # Direct deactivation ) -def test_macro_atom(clib, model_3lvlatom, packet, z_random, packet_params, get_rkstate, expected): +def test_macro_atom(packet, z_random, packet_params, get_rkstate, expected): packet.macro_atom_activation_level = packet_params['activation_level'] packet.current_shell_id = packet_params['shell_id'] rkstate = get_rkstate(z_random) - clib.macro_atom(byref(packet), byref(model_3lvlatom), byref(rkstate)) + macro_atom.macro_atom(byref(packet), byref(model_3lvlatom), byref(rkstate)) obtained_line_id = model_3lvlatom.last_line_interaction_out_id[packet.id] assert_equal(obtained_line_id, expected) @@ -773,7 +833,7 @@ def test_macro_atom(clib, model_3lvlatom, packet, z_random, packet_params, get_r """ Simple Tests: ---------------- -These test check very simple pices of code still work. +These test check very simple pieces of code still work. """ @pytest.mark.parametrize( @@ -791,31 +851,6 @@ def test_increment_Edotlu_estimator(clib, packet_params, line_idx, expected, pac assert_almost_equal(model.line_lists_Edotlu[line_idx], expected) -""" -Difficult Tests: ----------------- -The tests written further are more complex than previous tests. They require -proper design procedure. They are not taken up yet and intended to be -completed together in future. -""" - - - -@pytest.mark.skipif(True, reason="Yet to be written.") -def test_montecarlo_one_packet(packet, model, mt_state): - pass - - -@pytest.mark.skipif(True, reason="Yet to be written.") -def test_montecarlo_one_packet_loop(packet, model, mt_state): - pass - - -@pytest.mark.skipif(True, reason="Yet to be written.") -def test_montecarlo_main_loop(packet, model, mt_state): - pass - - """ Continuum Tests: @@ -825,295 +860,6 @@ def test_montecarlo_main_loop(packet, model, mt_state): """ - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - 't_electron', [2500., 15000.] -) -def test_sample_nu_free_free(clib, t_electron, packet, model, mt_state_seeded, expected_ff_emissivity): - model.t_electrons[packet.current_shell_id] = t_electron - clib.sample_nu_free_free.restype = c_double - - nu_bins, expected_emissivity = expected_ff_emissivity(t_electron) - - nus = [] - for _ in range(int(1e5)): - nu = clib.sample_nu_free_free(byref(packet), byref(model), byref(mt_state_seeded)) - nus.append(nu) - - obtained_emissivity, _ = np.histogram(nus, density=True, bins=nu_bins) - - assert_allclose(obtained_emissivity, expected_emissivity, rtol=1e-10) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['packet_params', 't_electrons', 'chi_ff_factor', 'expected'], - [({'nu': 4.5e14, 'mu': 0.0, 'current_shell_id': 1}, 15000, 2.0, 1.6746639430359494e-44), - ({'nu': 3.0e15, 'mu': 0.0, 'current_shell_id': 0}, 5000, 3.0, 1.1111111111107644e-46), - ({'nu': 3.0e15, 'mu': 0.4, 'current_shell_id': 0}, 10000, 4.0, 1.5638286016098277e-46)] -) -def test_calculate_chi_ff(clib, packet, model, packet_params, t_electrons, chi_ff_factor, expected): - packet.mu = packet_params['mu'] - packet.nu = packet_params['nu'] - packet.current_shell_id = packet_params['current_shell_id'] - packet.r = 1.04e17 - - model.t_electrons[packet_params['current_shell_id']] = t_electrons - model.chi_ff_factor[packet_params['current_shell_id']] = chi_ff_factor - - clib.calculate_chi_ff(byref(packet), byref(model)) - obtained = packet.chi_ff - - assert_equal(obtained, expected) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['continuum_status', 'z_random', 'packet_params', 'expected'], - [(CONTINUUM_OFF, 0.94183547596539363, - {'chi_c': 1.0, 'chi_th': 0.4, 'chi_bf': 0.5}, - 'montecarlo_thomson_scatter'), - - (CONTINUUM_ON, 0.22443743797312765, - {'chi_c': 1.0, 'chi_th': 0.4, 'chi_bf': 0.5}, - 'montecarlo_thomson_scatter'), - - (CONTINUUM_ON, 0.54510721066252377, - {'chi_c': 1.0, 'chi_th': 0.4, 'chi_bf': 0.5}, - 'montecarlo_bound_free_scatter'), - - (CONTINUUM_ON, 0.94183547596539363, - {'chi_c': 1.0, 'chi_th': 0.4, 'chi_bf': 0.5}, - 'montecarlo_free_free_scatter'), - - (CONTINUUM_ON, 0.22443743797312765, - {'chi_c': 1e2, 'chi_th': 1e1, 'chi_bf': 2e1}, - 'montecarlo_bound_free_scatter')] -) -def test_montecarlo_continuum_event_handler(clib, continuum_status, expected, z_random, - packet_params, packet, model, get_rkstate): - packet.chi_cont = packet_params['chi_c'] - packet.chi_th = packet_params['chi_th'] - packet.chi_bf = packet_params['chi_bf'] - model.cont_status = continuum_status - - rkstate = get_rkstate(z_random) - - clib.montecarlo_continuum_event_handler.restype = c_void_p - obtained = clib.montecarlo_continuum_event_handler(byref(packet), - byref(model), byref(rkstate)) - expected = getattr(clib, expected) - expected = cast(expected, c_void_p).value - - assert_equal(obtained, expected) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['nu', 'continuum_id', 'expected', 'bf_treatment'], - [(4.40e14, 1, 0.00, BoundFreeTreatment.LIN_INTERPOLATION), - (3.25e14, 1, 0.75, BoundFreeTreatment.LIN_INTERPOLATION), - (4.03e14, 0, 0.97, BoundFreeTreatment.LIN_INTERPOLATION), - (4.10e14 + 1e-1, 0, 0.90, BoundFreeTreatment.LIN_INTERPOLATION), - pytest.param(4.1e14, 0, 0.90, BoundFreeTreatment.LIN_INTERPOLATION, - marks=pytest.mark.xfail), - (6.50e14, 0, 0.23304506144742834, BoundFreeTreatment.HYDROGENIC), - (3.40e14, 2, 1.1170364339507428, BoundFreeTreatment.HYDROGENIC)] -) -def test_bf_cross_section(clib, nu, continuum_id, model_w_edges, expected, bf_treatment): - model_w_edges.bf_treatment = bf_treatment.value - - clib.bf_cross_section.restype = c_double - obtained = clib.bf_cross_section(byref(model_w_edges), continuum_id, c_double(nu)) - - assert_almost_equal(obtained, expected) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['packet_params', 'expected'], - [({'nu': 4.13e14, 'mu': 0.0, 'current_shell_id': 1}, - [3.2882087455641473, 0.0, 0.0]), - - ({'nu': 3.27e14, 'mu': 0.0, 'current_shell_id': 0}, - [0.0, 1.3992114634681028, 5.702548202131454]), - - ({'nu': 3.27e14, 'mu': -0.4, 'current_shell_id': 0}, - [0.0, 1.2670858, 5.4446587])] -) -def test_calculate_chi_bf(clib, packet_params, expected, packet, model_w_edges): - model_w_edges.l_pop = (c_double * 6)(*range(1, 7)) - model_w_edges.l_pop_r = (c_double * 6)(*np.linspace(0.1, 0.6, 6)) - model_w_edges.t_electrons[packet_params['current_shell_id']] = 1e4 - - packet.mu = packet_params['mu'] - packet.nu = packet_params['nu'] - packet.r = 1.04e17 - packet.current_shell_id = packet_params['current_shell_id'] - packet.chi_bf_tmp_partial = (c_double * model_w_edges.no_of_edges)() - - clib.calculate_chi_bf(byref(packet), byref(model_w_edges)) - - obtained_chi_bf_tmp = np.ctypeslib.as_array(packet.chi_bf_tmp_partial, shape=(model_w_edges.no_of_edges,)) - expected_chi_bf_tmp = np.array(expected) - expected_chi_bf = expected_chi_bf_tmp[expected_chi_bf_tmp > 0][-1] - - assert_almost_equal(obtained_chi_bf_tmp, expected_chi_bf_tmp) - assert_almost_equal(packet.chi_bf, expected_chi_bf) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['comov_energy', 'distance', 'chi_ff', 'no_of_updates', 'expected'], - [(1.3, 1.3e14, 3e-12, 1, 507.), - (0.9, 0.7e15, 2e-12, 10, 1.260e4), - (0.8, 0.8e13, 1.5e-12, 35, 336.)] -) -def test_increment_continuum_estimators_ff_heating_estimator(clib, packet, model_w_edges, comov_energy, distance, - chi_ff, no_of_updates, expected): - packet.chi_ff = chi_ff - - for _ in range(no_of_updates): - clib.increment_continuum_estimators(byref(packet), byref(model_w_edges), c_double(distance), - c_double(0), c_double(comov_energy)) - obtained = model_w_edges.ff_heating_estimator[packet.current_shell_id] - - assert_almost_equal(obtained, expected) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['comov_nus', 'expected'], - [([4.05e14, 4.17e14, 3.3e14, 3.2e14, 2.9e14], [2, 2, 3]), - ([4.15e15, 3.25e14, 3.3e14, 2.85e14, 2.9e14], [0, 2, 4])] -) -def test_increment_continuum_estimators_photo_ion_estimator_statistics(clib, packet, model_w_edges, comov_nus, expected): - for comov_nu in comov_nus: - clib.increment_continuum_estimators(byref(packet), byref(model_w_edges), c_double(1e13), - c_double(comov_nu), c_double(1.0)) - - no_of_edges = model_w_edges.no_of_edges - no_of_shells = model_w_edges.no_of_shells - - obtained = np.ctypeslib.as_array(model_w_edges.photo_ion_estimator_statistics, - shape=(no_of_edges * no_of_shells,)) - obtained = np.reshape(obtained, newshape=(no_of_shells, no_of_edges), order='F') - obtained = obtained[packet.current_shell_id] - expected = np.array(expected) - - assert_array_equal(obtained, expected) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['comov_energy', 'distance', 'comov_nus', 'expected'], - [(1.3, 1.3e14, [4.05e14, 2.65e14], - {"photo_ion": [0.39641975308641975, 0., 0.], - "stim_recomb": [0.056757061269242064, 0., 0.], - "bf_heating": [1.9820987654321076e12, 0., 0.], - "stim_recomb_cooling": [283784812699.75476, 0., 0.]}), - - (0.9, 0.7e15, [3.25e14, 2.85e14], - {"photo_ion": [0., 1.4538461538461538, 7.315141700404858], - "stim_recomb": [0., 0.3055802, 1.7292954], - "bf_heating": [0., 36346153846153.82, 156760323886639.69], - "stim_recomb_cooling": [0., 7639505724285.9746, 33907776077426.875]})] -) -def test_increment_continuum_estimators_bf_estimators(clib, packet, model_w_edges, comov_energy, - distance, comov_nus, expected): - for comov_nu in comov_nus: - clib.increment_continuum_estimators(byref(packet), byref(model_w_edges), c_double(distance), - c_double(comov_nu), c_double(comov_energy)) - - no_of_edges = model_w_edges.no_of_edges - no_of_shells = model_w_edges.no_of_shells - - for estim_name, expected_value in expected.items(): - obtained = np.ctypeslib.as_array(getattr(model_w_edges, estim_name + "_estimator"), - shape=(no_of_edges * no_of_shells,)) - obtained = np.reshape(obtained, newshape=(no_of_shells, no_of_edges), order='F') - obtained = obtained[packet.current_shell_id] - - assert_almost_equal(obtained, np.array(expected_value)) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['packet_params', 'expected_params'], - [({'nu_comov': 1.1e16, 'mu': 0.0, 'r': 1.4e14}, - {'next_line_id': 5, 'last_line': True, 'nu': 1.1e16, 'type_id': 3}), - - ({'nu_comov': 1.3e16, 'mu': 0.3, 'r': 7.5e14}, - {'next_line_id': 0, 'last_line': False, 'nu': 1.30018766e+16, 'type_id': 3}), - - ({'nu_comov': 1.24e16, 'mu': -0.3, 'r': 7.5e14}, - {'next_line_id': 2, 'last_line': False, 'nu': 1.23982106e+16, 'type_id': 4})] -) -def test_continuum_emission(clib, packet, model, mock_sample_nu, packet_params, expected_params, mt_state): - packet.nu = packet_params['nu_comov'] # Is returned by mock function mock_sample_nu - packet.mu = packet_params['mu'] - packet.r = packet_params['r'] - expected_interaction_out_type = expected_params['type_id'] - - clib.continuum_emission(byref(packet), byref(model), byref(mt_state), - mock_sample_nu, expected_interaction_out_type) - - obtained_next_line_id = packet.next_line_id - obtained_last_interaction_out_type = model.last_interaction_out_type[0] - - assert_equal(obtained_next_line_id, expected_params['next_line_id']) - assert_equal(packet.last_line, expected_params['last_line']) - assert_equal(expected_interaction_out_type, obtained_last_interaction_out_type) - assert_allclose(packet.nu, expected_params['nu'], rtol=1e-7) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['packet_params', 'expected'], - [({'next_line_id': 3, 'last_line': 0}, 1), - ({'next_line_id': 5, 'last_line': 1}, 0), - ({'next_line_id': 2, 'last_line': 0}, 0), - ({'next_line_id': 1, 'last_line': 0}, 1)] -) -def test_test_for_close_line(clib, packet, model, packet_params, expected): - packet.nu_line = model.line_list_nu[packet_params['next_line_id'] - 1] - packet.next_line_id = packet_params['next_line_id'] - - clib.test_for_close_line(byref(packet), byref(model)) - - assert_equal(expected, packet.close_line) - - -@pytest.mark.continuumtest -@pytest.mark.parametrize( - ['packet_params', 'z_random', 'expected'], - [({'current_continuum_id': 1, 'chi_bf_tmp_partial': [0.0, 0.23e13, 1.0e13]}, - 0.22443743797312765, 1), - - ({'current_continuum_id': 1, 'chi_bf_tmp_partial': [0.0, 0.23e10, 1.0e10]}, - 0.78961460371187597, 2), - - ({'current_continuum_id': 0, 'chi_bf_tmp_partial': [0.2e5, 0.5e5, 0.6e5, 1.0e5, 1.0e5]}, - 0.78961460371187597, 3)] -) -def test_montecarlo_bound_free_scatter_continuum_selection(clib, packet, model_3lvlatom, packet_params, - get_rkstate, z_random, expected): - rkstate = get_rkstate(z_random) - packet.current_continuum_id = packet_params['current_continuum_id'] - - chi_bf_tmp = packet_params['chi_bf_tmp_partial'] - packet.chi_bf_tmp_partial = (c_double * len(chi_bf_tmp))(*chi_bf_tmp) - packet.chi_bf = chi_bf_tmp[-1] - model_3lvlatom.no_of_edges = len(chi_bf_tmp) - - clib.montecarlo_bound_free_scatter(byref(packet), byref(model_3lvlatom), - c_double(1.e13), byref(rkstate)) - - assert_equal(packet.current_continuum_id, expected) - assert_equal(model_3lvlatom.last_line_interaction_in_id[packet.id], expected) - - @pytest.mark.continuumtest @pytest.mark.parametrize( ['mu', 'r', 'inv_t_exp', 'full_relativity'], @@ -1123,19 +869,19 @@ def test_montecarlo_bound_free_scatter_continuum_selection(clib, packet, model_3 (0.0, 7.5e14, 1 / 2.2e5, 1), (-0.7, 7.5e14, 1 / 5.2e5, 0)] ) -def test_frame_transformations(clib, packet, model, mu, r, +def test_frame_transformations(packet, model, mu, r, inv_t_exp, full_relativity): packet.r = r packet.mu = mu + mc.full_relativity = bool(full_relativity) model.inverse_time_explosion = inv_t_exp model.full_relativity = full_relativity - clib.rpacket_doppler_factor.restype = c_double - clib.rpacket_inverse_doppler_factor.restype = c_double - inverse_doppler_factor = clib.rpacket_inverse_doppler_factor(byref(packet), byref(model)) - clib.angle_aberration_CMF_to_LF(byref(packet), byref(model)) + inverse_doppler_factor = r_packet.get_inverse_doppler_factor(r, mu, 1/inv_t_exp) + r_packet.angle_aberration_CMF_to_LF(byref(packet), byref(model)) - doppler_factor = clib.rpacket_doppler_factor(byref(packet), byref(model)) + doppler_factor = r_packet.get_doppler_factor(r, mu, 1/inv_t_exp) + mc.full_relativity = False assert_almost_equal(doppler_factor * inverse_doppler_factor, 1.0) @@ -1149,17 +895,14 @@ def test_frame_transformations(clib, packet, model, mu, r, (0.0, 7.5e14, 1 / 2.2e5), (-0.7, 7.5e14, 1 / 5.2e5)] ) -def test_angle_transformation_invariance(clib, packet, model, +def test_angle_transformation_invariance(packet, model, mu, r, inv_t_exp): - packet.r = r - packet.mu = mu - model.inverse_time_explosion = inv_t_exp + packet = r_packet.RPacket(r, mu, packet.nu, packet.energy) model.full_relativity = 1 - clib.angle_aberration_LF_to_CMF.restype = c_double - clib.angle_aberration_CMF_to_LF(byref(packet), byref(model)) - mu_obtained = clib.angle_aberration_LF_to_CMF( - byref(packet), byref(model), c_double(packet.mu)) + mu1 = r_packet.angle_aberration_CMF_to_LF(packet, 1/inv_t_exp, mu) + mu_obtained = r_packet.angle_aberration_LF_to_CMF( + packet, 1/inv_t_exp, mu1) assert_almost_equal(mu_obtained, mu) @@ -1178,7 +921,7 @@ def test_angle_transformation_invariance(clib, packet, model, (-0.7, 7.5e14, 5.2e5, 1.0e15, 9.8e14), (-1.0, 6.3e14, 2.2e5, 6.0e12, 6.55e12)] ) -def test_compute_distance2line_relativistic(clib, mu, r, t_exp, nu, nu_line, +def test_compute_distance2line_relativistic(mu, r, t_exp, nu, nu_line, full_relativity, packet, model): packet.r = r packet.mu = mu @@ -1186,20 +929,17 @@ def test_compute_distance2line_relativistic(clib, mu, r, t_exp, nu, nu_line, packet.nu_line = nu_line model.inverse_time_explosion = 1 / t_exp model.time_explosion = t_exp - model.full_relativity = full_relativity + mc.full_relativity = bool(full_relativity) - clib.rpacket_doppler_factor.restype = c_double - - clib.compute_distance2line(byref(packet), byref(model)) - clib.move_packet(byref(packet), byref(model), c_double(packet.d_line)) + doppler_factor = r_packet.get_doppler_factor(r, mu, t_exp) + comov_nu = packet.nu * doppler_factor + distance = r_packet.calculate_distance_line( + packet, + comov_nu, nu_line, t_exp) + r_packet.move_r_packet(packet, distance, t_exp, model) - doppler_factor = clib.rpacket_doppler_factor(byref(packet), byref(model)) + doppler_factor = r_packet.get_doppler_factor(r, mu, t_exp) comov_nu = packet.nu * doppler_factor + mc.full_relativity = False assert_allclose(comov_nu, nu_line, rtol=1e-14) - - -@pytest.mark.continuumtest -@pytest.mark.skipif(True, reason="Yet to be written.") -def test_montecarlo_free_free_scatter(packet, model, mt_state): - pass From 35136c063d1fba642c5a52eb0baa35ac8ff2ed15 Mon Sep 17 00:00:00 2001 From: Arjun Savel Date: Mon, 13 Jul 2020 21:06:29 -0700 Subject: [PATCH 07/21] Move location of coverage file --- tardis/tests/coveragerc => coveragerc | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename tardis/tests/coveragerc => coveragerc (100%) diff --git a/tardis/tests/coveragerc b/coveragerc similarity index 100% rename from tardis/tests/coveragerc rename to coveragerc From 49b39d4fd317786c8d587e1031b8ffd99cfe6b09 Mon Sep 17 00:00:00 2001 From: Arjun Savel Date: Mon, 13 Jul 2020 21:07:03 -0700 Subject: [PATCH 08/21] Remove extraneous imports --- tardis/tests/test_tardis_full_formal_integral.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tardis/tests/test_tardis_full_formal_integral.py b/tardis/tests/test_tardis_full_formal_integral.py index d3de68501c4..4704491871b 100644 --- a/tardis/tests/test_tardis_full_formal_integral.py +++ b/tardis/tests/test_tardis_full_formal_integral.py @@ -1,13 +1,11 @@ import os import pytest -import numpy as np import numpy.testing as npt from astropy import units as u from astropy.tests.helper import assert_quantity_allclose from tardis.simulation.base import Simulation from tardis.io.config_reader import Configuration -import astropy config_line_modes = ['downbranch', 'macroatom'] interpolate_shells = [-1, 30] From 9f246105e5668ef439736d825d5b43192221c3ee Mon Sep 17 00:00:00 2001 From: Arjun Savel Date: Thu, 16 Jul 2020 14:15:54 -0700 Subject: [PATCH 09/21] Refer to previous plasma object for atomic data --- tardis/montecarlo/formal_integral.py | 31 ++++++++++++++-------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/tardis/montecarlo/formal_integral.py b/tardis/montecarlo/formal_integral.py index 3f793df58e5..9064df54e1d 100644 --- a/tardis/montecarlo/formal_integral.py +++ b/tardis/montecarlo/formal_integral.py @@ -35,7 +35,9 @@ class FormalIntegrator(object): def __init__(self, model, plasma, runner, points=1000): self.model = model - self.plasma = numba_plasma_initialize(plasma) + if plasma: + self.plasma = numba_plasma_initialize(plasma) + self.atomic_data = plasma.atomic_data self.runner = runner self.points = points @@ -120,25 +122,24 @@ def make_source_function(self): """ model = self.model - plasma = self.plasma runner = self.runner - atomic_data = self.plasma.atomic_data - macro_ref = plasma.macro_block_references - macro_data = atomic_data.macro_atom_data - no_lvls = len(atomic_data.levels) + macro_ref = self.plasma.macro_block_references + macro_data = self.atomic_data.macro_atom_data + + no_lvls = len(self.tomic_data.levels) no_shells = len(model.w) if runner.line_interaction_type == 'macroatom': internal_jump_mask = (macro_data.transition_type >= 0).values ma_int_data = macro_data[internal_jump_mask] - internal = plasma.transition_probabilities[internal_jump_mask] + internal = self.plasma.transition_probabilities[internal_jump_mask] source_level_idx = ma_int_data.source_level_idx.values destination_level_idx = ma_int_data.destination_level_idx.values Edotlu_norm_factor = (1 / (runner.time_of_simulation * model.volume)) - exptau = 1 - np.exp(- plasma.tau_sobolevs) + exptau = 1 - np.exp(-self.plasma.tau_sobolev) Edotlu = Edotlu_norm_factor * exptau * runner.Edotlu_estimator # The following may be achieved by calling the appropriate plasma @@ -150,7 +151,7 @@ def make_source_function(self): # the transition l->u Jbluelu = runner.j_blue_estimator * Jbluelu_norm_factor - upper_level_index = atomic_data.lines.index.droplevel('level_number_lower') + upper_level_index = self.atomic_data.lines.index.droplevel('level_number_lower') e_dot_lu = pd.DataFrame(Edotlu, index=upper_level_index) e_dot_u = e_dot_lu.groupby(level=[0, 1, 2]).sum() e_dot_u_src_idx = macro_ref.loc[e_dot_u.index].references_idx.values @@ -170,12 +171,12 @@ def make_source_function(self): C_frame[shell] = sp.linalg.spsolve(inv_N.T, e_dot_u_vec) e_dot_u.index.names = ['atomic_number', 'ion_number', 'source_level_number'] # To make the q_ul e_dot_u product work, could be cleaner - transitions = atomic_data.macro_atom_data[atomic_data.macro_atom_data.transition_type == -1].copy() + transitions = self.atomic_data.macro_atom_data[self.atomic_data.macro_atom_data.transition_type == -1].copy() transitions_index = transitions.set_index(['atomic_number', 'ion_number', 'source_level_number']).index.copy() - tmp = plasma.transition_probabilities[(atomic_data.macro_atom_data.transition_type == -1).values] + tmp = self.plasma.transition_probabilities[(self.atomic_data.macro_atom_data.transition_type == -1).values] q_ul = tmp.set_index(transitions_index) t = model.time_explosion.value - lines = atomic_data.lines.set_index('line_id') + lines = self.atomic_data.lines.set_index('line_id') wave = lines.wavelength_cm.loc[transitions.transition_line_id].values.reshape(-1,1) if runner.line_interaction_type == 'macroatom': e_dot_u = C_frame.loc[e_dot_u.index] @@ -186,15 +187,15 @@ def make_source_function(self): # Jredlu should already by in the correct order, i.e. by wavelength of # the transition l->u (similar to Jbluelu) - Jredlu = Jbluelu * np.exp(-plasma.tau_sobolevs.values) + att_S_ul + Jredlu = Jbluelu * np.exp(-self.plasma.tau_sobolev.values) + att_S_ul if self.interpolate_shells > 0: att_S_ul, Jredlu, Jbluelu, e_dot_u = self.interpolate_integrator_quantities( att_S_ul, Jredlu, Jbluelu, e_dot_u) else: runner.r_inner_i = runner.r_inner_cgs runner.r_outer_i = runner.r_outer_cgs - runner.tau_sobolevs_integ = plasma.tau_sobolevs.values - runner.electron_densities_integ = plasma.electron_densities.values + runner.tau_sobolevs_integ = self.plasma.tau_sobolevs.values + runner.electron_densities_integ = self.plasma.electron_densities.values return att_S_ul, Jredlu, Jbluelu, e_dot_u From 416c4bc5c3bf79ac0b7e07d4c24eafb64440b789 Mon Sep 17 00:00:00 2001 From: Arjun Savel Date: Mon, 20 Jul 2020 10:11:53 -0400 Subject: [PATCH 10/21] Debug formal integral to ensure correct plasma is referenced --- tardis/montecarlo/formal_integral.py | 75 +++++++++++++++------------- 1 file changed, 40 insertions(+), 35 deletions(-) diff --git a/tardis/montecarlo/formal_integral.py b/tardis/montecarlo/formal_integral.py index 9064df54e1d..75249793422 100644 --- a/tardis/montecarlo/formal_integral.py +++ b/tardis/montecarlo/formal_integral.py @@ -6,6 +6,7 @@ from astropy import units as u from tardis import constants as const from numba import jitclass, njit +import pdb from tardis.montecarlo.montecarlo_numba import njit_dict @@ -19,6 +20,7 @@ M_PI = np.arccos(-1) KB_CGS = 1.3806488e-16 H_CGS = 6.62606957e-27 +SIGMA_THOMSON = 6.652486e-25 class IntegrationError(Exception): pass @@ -38,6 +40,7 @@ def __init__(self, model, plasma, runner, points=1000): if plasma: self.plasma = numba_plasma_initialize(plasma) self.atomic_data = plasma.atomic_data + self.original_plasma = plasma self.runner = runner self.points = points @@ -124,16 +127,16 @@ def make_source_function(self): model = self.model runner = self.runner - macro_ref = self.plasma.macro_block_references + macro_ref = self.atomic_data.macro_atom_references macro_data = self.atomic_data.macro_atom_data - no_lvls = len(self.tomic_data.levels) + no_lvls = len(self.atomic_data.levels) no_shells = len(model.w) if runner.line_interaction_type == 'macroatom': internal_jump_mask = (macro_data.transition_type >= 0).values ma_int_data = macro_data[internal_jump_mask] - internal = self.plasma.transition_probabilities[internal_jump_mask] + internal = self.original_plasma.transition_probabilities[internal_jump_mask] source_level_idx = ma_int_data.source_level_idx.values destination_level_idx = ma_int_data.destination_level_idx.values @@ -170,10 +173,11 @@ def make_source_function(self): e_dot_u_vec[e_dot_u_src_idx] = e_dot_u[shell].values C_frame[shell] = sp.linalg.spsolve(inv_N.T, e_dot_u_vec) + e_dot_u.index.names = ['atomic_number', 'ion_number', 'source_level_number'] # To make the q_ul e_dot_u product work, could be cleaner - transitions = self.atomic_data.macro_atom_data[self.atomic_data.macro_atom_data.transition_type == -1].copy() + transitions = self.original_plasma.atomic_data.macro_atom_data[self.original_plasma.atomic_data.macro_atom_data.transition_type == -1].copy() transitions_index = transitions.set_index(['atomic_number', 'ion_number', 'source_level_number']).index.copy() - tmp = self.plasma.transition_probabilities[(self.atomic_data.macro_atom_data.transition_type == -1).values] + tmp = self.original_plasma.transition_probabilities[(self.atomic_data.macro_atom_data.transition_type == -1).values] q_ul = tmp.set_index(transitions_index) t = model.time_explosion.value lines = self.atomic_data.lines.set_index('line_id') @@ -187,15 +191,15 @@ def make_source_function(self): # Jredlu should already by in the correct order, i.e. by wavelength of # the transition l->u (similar to Jbluelu) - Jredlu = Jbluelu * np.exp(-self.plasma.tau_sobolev.values) + att_S_ul + Jredlu = Jbluelu * np.exp(-self.plasma.tau_sobolev) + att_S_ul if self.interpolate_shells > 0: att_S_ul, Jredlu, Jbluelu, e_dot_u = self.interpolate_integrator_quantities( att_S_ul, Jredlu, Jbluelu, e_dot_u) else: runner.r_inner_i = runner.r_inner_cgs runner.r_outer_i = runner.r_outer_cgs - runner.tau_sobolevs_integ = self.plasma.tau_sobolevs.values - runner.electron_densities_integ = self.plasma.electron_densities.values + runner.tau_sobolevs_integ = self.plasma.tau_sobolev + runner.electron_densities_integ = self.plasma.electron_density return att_S_ul, Jredlu, Jbluelu, e_dot_u @@ -215,12 +219,12 @@ def interpolate_integrator_quantities(self, att_S_ul, Jredlu, r_middle_integ = (r_integ[:-1] + r_integ[1:]) / 2. runner.electron_densities_integ = interp1d( - r_middle, plasma.electron_densities, + r_middle, plasma.electron_density, fill_value='extrapolate', kind='nearest')(r_middle_integ) # Assume tau_sobolevs to be constant within a shell # (as in the MC simulation) runner.tau_sobolevs_integ = interp1d( - r_middle, plasma.tau_sobolevs, + r_middle, plasma.tau_sobolev, fill_value='extrapolate', kind='nearest')(r_middle_integ) att_S_ul = interp1d( r_middle, att_S_ul, fill_value='extrapolate')(r_middle_integ) @@ -270,23 +274,23 @@ def _formal_integral(self, iT, inu, inu_size, att_S_ul, Jred_lu, Jblue_lu, N): size_tau = size_line * size_shell finished_nus = 0 - R_ph = self.runner.r_inner[0] - R_max = self.runner.r_outer[size_shell - 1] + R_ph = self.runner.r_inner_i[0] + R_max = self.runner.r_outer_i[size_shell - 1] pp = np.zeros(N) # check exp_tau = np.zeros(size_tau) # TODO: multiprocessing offset = 0 size_z = 0 + z = np.zeros(2 * self.model.no_of_shells) idx_nu_start = 0 direction = 0 I_nu = np.zeros(N) - shell_id = np.zeros(2 * self.runner.no_of_shells_i) # check + shell_id = np.zeros(2 * self.model.no_of_shells) # instantiate more variables here, maybe? # prepare exp_tau - for i in range(size_tau): - exp_tau[i] = np.exp(-self.model.line_lists_tau_sobolevs_i[i]) # check - pp = self.calculate_p_values(R_max, N, pp) + exp_tau = np.exp(-self.plasma.tau_sobolev) # check + pp = calculate_p_values(R_max, N, pp) # done with instantiation # now loop over wavelength in spectrum @@ -309,7 +313,7 @@ def _formal_integral(self, iT, inu, inu_size, att_S_ul, Jred_lu, Jblue_lu, N): # find first contributing lines nu_start = nu * z[0] nu_end = nu * z[1] - idx_nu_start = line_search(self.runner.line_list_nu, + idx_nu_start = line_search(self.plasma.line_list_nu, nu_start, size_line, idx_nu_start) offset = shell_id[0] * size_line @@ -317,31 +321,31 @@ def _formal_integral(self, iT, inu, inu_size, att_S_ul, Jred_lu, Jblue_lu, N): zstart = self.model.time_explosion / C_INV * (1. - z[0]) # Initialize "pointers" - pline = self.runner.line_list_nu + idx_nu_start; - pexp_tau = exp_tau + offset + idx_nu_start; - patt_S_ul = att_S_ul + offset + idx_nu_start; - pJred_lu = Jred_lu + offset + idx_nu_start; - pJblue_lu = Jblue_lu + offset + idx_nu_start; + pline = self.plasma.line_list_nu + idx_nu_start + pexp_tau = exp_tau + offset + idx_nu_start + patt_S_ul = att_S_ul + offset + idx_nu_start + pJred_lu = Jred_lu + offset + idx_nu_start + pJblue_lu = Jblue_lu + offset + idx_nu_start # flag for first contribution to integration on current p-ray first = 1 # loop over all interactions for i in range(size_z - 1): - escat_op = self.model.electron_densities_i[shell_id[i]] * sigma_thomson # change sigma_thomson + escat_op = self.plasma.electron_density[int(shell_id[i])] * SIGMA_THOMSON nu_end = nu * z[i + 1] - while pline < self.model.ine_list_nu + size_line: # check ;pline + while np.all(pline < self.plasma.line_list_nu + size_line): # check all condition # increment all pointers simulatenously pline += 1 pexp_tau += 1 patt_S_ul += 1 pJblue_lu += 1 - if (pline[0] < nu_end): + if (pline[0] < nu_end.value): break # calculate e-scattering optical depth to next resonance point - zend = self.model.time_explosion / C_INV * (1. - pline / nu) # check + zend = self.model.time_explosion / C_INV * (1. - pline[0] / nu.value) # check if first == 1: # first contribution to integration @@ -349,26 +353,27 @@ def _formal_integral(self, iT, inu, inu_size, att_S_ul, Jred_lu, Jblue_lu, N): # by boundary conditions) is not in Lucy 1999; # should be re-examined carefully escat_contrib += (zend - zstart) * escat_op * ( - pJblue_lu - I_nu[p_idx]); + pJblue_lu[0] - I_nu[p_idx]); first = 0; else: # Account for e-scattering, c.f. Eqs 27, 28 in Lucy 1999 - Jkkp = 0.5 * (pJred_lu + pJblue_lu); + Jkkp = 0.5 * (pJred_lu[0] + pJblue_lu[0]); escat_contrib += (zend - zstart) * escat_op * ( Jkkp - I_nu[p_idx]) # this introduces the necessary ffset of one element between # pJblue_lu and pJred_lu pJred_lu += 1 - I_nu[p_idx] = I_nu[p_idx] + escat_contrib; + # pdb.set_trace() + I_nu[p_idx] = I_nu[p_idx] + escat_contrib.value # // Lucy 1999, Eq 26 - I_nu[p_idx] = I_nu[p_idx] * (pexp_tau) + patt_S_ul # check about taking about asterisks beforehand elsewhere + I_nu[p_idx] = I_nu[p_idx] * (pexp_tau[0][0]) + patt_S_ul[0] # check about taking about asterisks beforehand elsewhere # // reset e-scattering opacity escat_contrib = 0 zstart = zend # calculate e-scattering optical depth to grid cell boundary - Jkkp = 0.5 * (pJred_lu + pJblue_lu) + Jkkp = 0.5 * (pJred_lu[0] + pJblue_lu[0]) zend = self.model.time_explosion / C_INV * (1. - nu_end / nu) # check escat_contrib += (zend - zstart) * escat_op * ( Jkkp - I_nu[p_idx]) @@ -402,14 +407,14 @@ def populate_z(self, p, oz, oshell_id): :oshell_id: (int64) will be set with the corresponding shell_ids """ # abbreviations - r = self.model.r_outer_i - N = self.model.no_of_shells_i # check + r = self.runner.r_outer_i + N = self.model.no_of_shells # check print(N) inv_t = 1/self.model.time_explosion z = 0 offset = N - if p <= self.model.r_inner_i[0]: + if p <= self.runner.r_inner_i[0]: # intersect the photosphere for i in range(N): oz[i] = 1 - self.calculate_z(r[i], p, inv_t) @@ -449,7 +454,7 @@ def calculate_z(self, r, p, inv_t): :inv_t: (double) inverse time_explosio is needed to norm to unit-length """ if r > p: - return np.sqrt(r * r - p * p) * C_INV * inv_t + return np.sqrt(r * r - p * p) * C_INV * inv_t.value else: return 0 From 35fa02d8395a6c54541afd809071423dfd16f166 Mon Sep 17 00:00:00 2001 From: Arjun Savel Date: Mon, 20 Jul 2020 10:25:17 -0400 Subject: [PATCH 11/21] Adding runner to test fixture --- tardis/montecarlo/tests/conftest.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/tardis/montecarlo/tests/conftest.py b/tardis/montecarlo/tests/conftest.py index 054f6b455c4..30b4cd121a6 100644 --- a/tardis/montecarlo/tests/conftest.py +++ b/tardis/montecarlo/tests/conftest.py @@ -1,5 +1,7 @@ import os import pytest +from tardis.io import config_reader + from ctypes import ( CDLL, @@ -15,6 +17,7 @@ CONTINUUM_OFF, BoundFreeTreatment ) +from tardis.montecarlo.base import MontecarloRunner @pytest.fixture(scope="function") @@ -138,3 +141,10 @@ def mt_state_seeded(clib, mt_state): seed = 23111963 clib.rk_seed(seed, byref(mt_state)) return mt_state + +@pytest.fixture(scope="function") +def runner(): + config_fname = 'tardis/io/tests/data/tardis_configv1_verysimply.yml' + config = config_reader.Configuration.from_yaml(config_fname) + runner = MontecarloRunner.from_config(config) + return runner From 7274a14a1902cba19fe7ef6f47c73661fe7f48c3 Mon Sep 17 00:00:00 2001 From: Arjun Savel Date: Tue, 21 Jul 2020 10:20:37 -0400 Subject: [PATCH 12/21] Add capacity to use Estimators --- tardis/montecarlo/tests/test_montecarlo.py | 61 +++++++++++++--------- 1 file changed, 35 insertions(+), 26 deletions(-) diff --git a/tardis/montecarlo/tests/test_montecarlo.py b/tardis/montecarlo/tests/test_montecarlo.py index c8ca25dde24..2ca4c6397c5 100644 --- a/tardis/montecarlo/tests/test_montecarlo.py +++ b/tardis/montecarlo/tests/test_montecarlo.py @@ -52,7 +52,7 @@ import tardis.montecarlo.montecarlo_numba.r_packet as r_packet import tardis.montecarlo.montecarlo_configuration as mc from tardis import constants as const -from tardis.montecarlo.montecarlo_numba import numba_interface +from tardis.montecarlo.montecarlo_numba.numba_interface import Estimators from tardis.montecarlo.montecarlo_numba import macro_atom C_SPEED_OF_LIGHT = const.c.to('cm/s').value @@ -288,7 +288,6 @@ def test_reverse_binary_search(x, x_insert, imin, imax, expected_params): ) def test_line_search(nu, nu_insert, number_of_lines, expected_params): # nu = (c_double * number_of_lines)(*nu) - nu_insert = c_double(nu_insert) obtained_result = 0 obtained_tardis_error = formal_integral.line_search( @@ -660,7 +659,7 @@ def test_compute_distance2line(packet_params, expected_params, packet, model): # assert_almost_equal(packet.d_cont, expected_params['d_cont']) -# @pytest.mark.parametrize('full_relativity', [1, 0]) +@pytest.mark.parametrize('full_relativity', [1, 0]) @pytest.mark.parametrize( ['packet_params', 'expected_params'], [({'nu': 0.4, 'mu': 0.3, 'energy': 0.9, 'r': 7.5e14}, @@ -672,32 +671,41 @@ def test_compute_distance2line(packet_params, expected_params, packet, model): 'j': 5001298975563.031, 'nubar': 3001558973156.1387})] ) def test_move_packet(packet_params, expected_params, - packet, model): - # TODO: test relativity + packet, model, full_relativity): + distance = 1e13 packet.nu = packet_params['nu'] packet.mu = packet_params['mu'] packet.energy = packet_params['energy'] packet.r = packet_params['r'] # model.full_relativity = full_relativity + mc.full_relativity = full_relativity doppler_factor = r_packet.get_doppler_factor(packet.r, packet.mu, model.time_explosion) + numba_estimator = Estimators( + packet_params['j'], + packet_params['nu_bar'], + 0, + 0) r_packet.move_r_packet( - packet, distance, time_explosion, numba_estimator) + packet, distance, model.time_explosion, numba_estimator) assert_almost_equal(packet.mu, expected_params['mu']) assert_almost_equal(packet.r, expected_params['r']) + expected_j = expected_params['j'] expected_nubar = expected_params['nubar'] - # if full_relativity: - # expected_j *= doppler_factor - # expected_nubar *= doppler_factor + if full_relativity: + expected_j *= doppler_factor + expected_nubar *= doppler_factor + + mc.full_relativity = False - assert_allclose(model.js[packet.current_shell_id], + assert_allclose(numba_estimator.j_estimator[packet.current_shell_id], expected_j, rtol=5e-7) - assert_allclose(model.nubars[packet.current_shell_id], + assert_allclose(numba_estimator.nu_bar_estimator[packet.current_shell_id], expected_nubar, rtol=5e-7) @@ -820,10 +828,12 @@ def test_montecarlo_line_scatter(clib, packet_params, expected_params, packet, m ) def test_macro_atom(packet, z_random, packet_params, get_rkstate, expected): packet.macro_atom_activation_level = packet_params['activation_level'] + packet = r_packet.RPacket(r=packet.r, mu=packet.mu, nu=packet.nu, + energy=packet.energy) packet.current_shell_id = packet_params['shell_id'] rkstate = get_rkstate(z_random) - macro_atom.macro_atom(byref(packet), byref(model_3lvlatom), byref(rkstate)) + macro_atom.macro_atom(packet, numba_plasma) obtained_line_id = model_3lvlatom.last_line_interaction_out_id[packet.id] assert_equal(obtained_line_id, expected) @@ -869,16 +879,14 @@ def test_increment_Edotlu_estimator(clib, packet_params, line_idx, expected, pac (0.0, 7.5e14, 1 / 2.2e5, 1), (-0.7, 7.5e14, 1 / 5.2e5, 0)] ) -def test_frame_transformations(packet, model, mu, r, +def test_frame_transformations(packet, mu, r, inv_t_exp, full_relativity): - packet.r = r - packet.mu = mu + packet = r_packet.RPacket(r=r, mu=mu, energy=packet.energy, nu=packet.nu) mc.full_relativity = bool(full_relativity) - model.inverse_time_explosion = inv_t_exp - model.full_relativity = full_relativity + mc.full_relativity = full_relativity inverse_doppler_factor = r_packet.get_inverse_doppler_factor(r, mu, 1/inv_t_exp) - r_packet.angle_aberration_CMF_to_LF(byref(packet), byref(model)) + r_packet.angle_aberration_CMF_to_LF(packet, 1/inv_t_exp, packet.mu) doppler_factor = r_packet.get_doppler_factor(r, mu, 1/inv_t_exp) mc.full_relativity = False @@ -922,13 +930,14 @@ def test_angle_transformation_invariance(packet, model, (-1.0, 6.3e14, 2.2e5, 6.0e12, 6.55e12)] ) def test_compute_distance2line_relativistic(mu, r, t_exp, nu, nu_line, - full_relativity, packet, model): - packet.r = r - packet.mu = mu - packet.nu = nu - packet.nu_line = nu_line - model.inverse_time_explosion = 1 / t_exp - model.time_explosion = t_exp + full_relativity, packet, runner): + packet = r_packet.RPacket(r=r, nu=nu, mu=mu, energy=packet.energy) + # packet.nu_line = nu_line + numba_estimator = Estimators( + runner.j_estimator, + runner.nu_bar_estimator, + runner.j_blue_estimator, + runner.Edotlu_estimator) mc.full_relativity = bool(full_relativity) doppler_factor = r_packet.get_doppler_factor(r, mu, t_exp) @@ -936,7 +945,7 @@ def test_compute_distance2line_relativistic(mu, r, t_exp, nu, nu_line, distance = r_packet.calculate_distance_line( packet, comov_nu, nu_line, t_exp) - r_packet.move_r_packet(packet, distance, t_exp, model) + r_packet.move_r_packet(packet, distance, t_exp, numba_estimator) doppler_factor = r_packet.get_doppler_factor(r, mu, t_exp) comov_nu = packet.nu * doppler_factor From a5db1fd313ab2d49cda2d3f4ad716a8eaeda44f7 Mon Sep 17 00:00:00 2001 From: Jack O'Brien Date: Wed, 14 Oct 2020 11:22:59 -0400 Subject: [PATCH 13/21] Update r_packet.py --- .../montecarlo/montecarlo_numba/r_packet.py | 26 ++----------------- 1 file changed, 2 insertions(+), 24 deletions(-) diff --git a/tardis/montecarlo/montecarlo_numba/r_packet.py b/tardis/montecarlo/montecarlo_numba/r_packet.py index 0a0eb61e6b0..e1e90bcadf8 100644 --- a/tardis/montecarlo/montecarlo_numba/r_packet.py +++ b/tardis/montecarlo/montecarlo_numba/r_packet.py @@ -154,17 +154,10 @@ def calculate_tau_electron(electron_density, distance): @njit(**njit_dict) def get_doppler_factor(r, mu, time_explosion): -<<<<<<< HEAD inv_c = 1 / C_SPEED_OF_LIGHT inv_t = 1/time_explosion beta = r * inv_t * inv_c if not numba_config.ENABLE_FULL_RELATIVITY: -======= - beta = r / (time_explosion * C_SPEED_OF_LIGHT) - if beta > 1.: # TODO: check speed penalty - raise SuperluminalError - if not montecarlo_configuration.full_relativity: ->>>>>>> upstream/pr/1217 return get_doppler_factor_partial_relativity(mu, beta) else: return get_doppler_factor_full_relativity(mu, beta) @@ -175,27 +168,16 @@ def get_doppler_factor_partial_relativity(mu, beta): @njit(**njit_dict) def get_doppler_factor_full_relativity(mu, beta): -<<<<<<< HEAD return (1.0 - mu * beta) / math.sqrt(1 - beta * beta) -======= - # todo: check numpy vs other sqrt - return (1.0 - mu * beta) / np.sqrt(1 - beta * beta) ->>>>>>> upstream/pr/1217 + @njit(**njit_dict) def get_inverse_doppler_factor(r, mu, time_explosion): -<<<<<<< HEAD inv_c = 1 / C_SPEED_OF_LIGHT inv_t = 1 / time_explosion beta = r * inv_t * inv_c if not numba_config.ENABLE_FULL_RELATIVITY: -======= - beta = (r / time_explosion) / C_SPEED_OF_LIGHT - if beta > 1.: - raise SuperluminalError - if not montecarlo_configuration.full_relativity: ->>>>>>> upstream/pr/1217 return get_inverse_doppler_factor_partial_relativity(mu, beta) else: return get_inverse_doppler_factor_full_relativity(mu, beta) @@ -206,12 +188,8 @@ def get_inverse_doppler_factor_partial_relativity(mu, beta): @njit(**njit_dict) def get_inverse_doppler_factor_full_relativity(mu, beta): -<<<<<<< HEAD return (1.0 + mu * beta) / math.sqrt(1 - beta * beta) -======= - # TODO: check numpy vs. other sqrt - return (1.0 + mu * beta) / np.sqrt(1 - beta * beta) ->>>>>>> upstream/pr/1217 + @njit(**njit_dict) def get_random_mu(): From 5492d4b8da64b3bd4047cdaa5bc16af72e8bc93b Mon Sep 17 00:00:00 2001 From: Andrew Fullard Date: Wed, 14 Oct 2020 11:29:18 -0400 Subject: [PATCH 14/21] Added empty test files --- tardis/montecarlo/montecarlo_numba/tests/test_base.py | 0 tardis/montecarlo/montecarlo_numba/tests/test_interaction.py | 0 tardis/montecarlo/montecarlo_numba/tests/test_macro_atom.py | 0 .../montecarlo/montecarlo_numba/tests/test_montecarlo_logger.py | 0 tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py | 0 tardis/montecarlo/montecarlo_numba/tests/test_r_packet.py | 0 .../montecarlo/montecarlo_numba/tests/test_single_packet_loop.py | 0 tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py | 0 8 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 tardis/montecarlo/montecarlo_numba/tests/test_base.py create mode 100644 tardis/montecarlo/montecarlo_numba/tests/test_interaction.py create mode 100644 tardis/montecarlo/montecarlo_numba/tests/test_macro_atom.py create mode 100644 tardis/montecarlo/montecarlo_numba/tests/test_montecarlo_logger.py create mode 100644 tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py create mode 100644 tardis/montecarlo/montecarlo_numba/tests/test_r_packet.py create mode 100644 tardis/montecarlo/montecarlo_numba/tests/test_single_packet_loop.py create mode 100644 tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_base.py b/tardis/montecarlo/montecarlo_numba/tests/test_base.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_interaction.py b/tardis/montecarlo/montecarlo_numba/tests/test_interaction.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_macro_atom.py b/tardis/montecarlo/montecarlo_numba/tests/test_macro_atom.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_montecarlo_logger.py b/tardis/montecarlo/montecarlo_numba/tests/test_montecarlo_logger.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py b/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_r_packet.py b/tardis/montecarlo/montecarlo_numba/tests/test_r_packet.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_single_packet_loop.py b/tardis/montecarlo/montecarlo_numba/tests/test_single_packet_loop.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py b/tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py new file mode 100644 index 00000000000..e69de29bb2d From af140a7a72c1d15429efcab6a6b2501d96197a72 Mon Sep 17 00:00:00 2001 From: Andrew Fullard Date: Wed, 14 Oct 2020 12:30:47 -0400 Subject: [PATCH 15/21] Packet tests --- .../montecarlo/montecarlo_numba/r_packet.py | 4 +- .../montecarlo_numba/tests/test_packet.py | 208 ++++++++++++++++++ .../montecarlo_numba/tests/test_r_packet.py | 0 3 files changed, 210 insertions(+), 2 deletions(-) create mode 100644 tardis/montecarlo/montecarlo_numba/tests/test_packet.py delete mode 100644 tardis/montecarlo/montecarlo_numba/tests/test_r_packet.py diff --git a/tardis/montecarlo/montecarlo_numba/r_packet.py b/tardis/montecarlo/montecarlo_numba/r_packet.py index e1e90bcadf8..d487967de10 100644 --- a/tardis/montecarlo/montecarlo_numba/r_packet.py +++ b/tardis/montecarlo/montecarlo_numba/r_packet.py @@ -120,8 +120,8 @@ def calculate_distance_line( distance = (nu_diff / nu) * C_SPEED_OF_LIGHT * time_explosion else: print('WARNING: nu difference is less than 0.0') - #raise MonteCarloException('nu difference is less than 0.0; for more' - # ' information, see print statement beforehand') + raise MonteCarloException('nu difference is less than 0.0; for more' + ' information, see print statement beforehand') if numba_config.ENABLE_FULL_RELATIVITY: return calculate_distance_line_full_relativity(nu_line, nu, diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py new file mode 100644 index 00000000000..d37702e0509 --- /dev/null +++ b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py @@ -0,0 +1,208 @@ +import os +import pytest +import numpy as np +import pandas as pd +import tardis.montecarlo.formal_integral as formal_integral +import tardis.montecarlo.montecarlo_numba.r_packet as r_packet +import tardis.montecarlo.montecarlo_configuration as mc +import tardis.montecarlo.montecarlo_numba.numba_interface as numba_interface +from tardis import constants as const +from tardis.montecarlo.montecarlo_numba.numba_interface import Estimators +from tardis.montecarlo.montecarlo_numba import macro_atom +C_SPEED_OF_LIGHT = const.c.to('cm/s').value + +from numpy.testing import ( + assert_equal, + assert_almost_equal, + assert_array_equal, + assert_allclose + ) + +@pytest.fixture(scope="function") +def packet(): + return r_packet.RPacket( + r = 7.5e14, + nu = 0.4, + mu = 0.3, + energy = 0.9, + seed = 1963, + index = 0, + is_close_line = 0 + ) + +@pytest.fixture(scope="function") +def model(): + return numba_interface.NumbaModel( + r_inner = np.array([6.912e14, 8.64e14], dtype=np.float64), + r_outer = np.array([8.64e14, 1.0368e15], dtype=np.float64), + time_explosion = 5.2e7 + ) + + +@pytest.mark.parametrize( + ['mu', 'r', 'inv_t_exp', 'expected'], + [(0.3, 7.5e14, 1 / 5.2e7, 0.9998556693818854), + (-.3, 0, 1 / 2.6e7, 1.0), + (0, 1, 1 / 2.6e7, 1.0)] +) +def test_get_doppler_factor(mu, r, inv_t_exp, expected): + # Set the params from test cases here + # TODO: add relativity tests + time_explosion = 1/inv_t_exp + + # Perform any other setups just before this, they can be additional calls + # to other methods or introduction of some temporary variables + + obtained = r_packet.get_doppler_factor(r, mu, time_explosion) + + + # Perform required assertions + assert_almost_equal(obtained, expected) + + +@pytest.mark.parametrize( + ['mu', 'r', 'inv_t_exp', 'expected'], + [(0.3, 7.5e14, 1 / 5.2e7, 1/0.9998556693818854), + (-.3, 0, 1 / 2.6e7, 1.0), + (0, 1, 1 / 2.6e7, 1.0)] +) +def test_get_inverse_doppler_factor(mu, r, inv_t_exp, expected): + # Set the params from test cases here + # TODO: add relativity tests + time_explosion = 1/inv_t_exp + + # Perform any other setups just before this, they can be additional calls + # to other methods or introduction of some temporary variables + + obtained = r_packet.get_inverse_doppler_factor(r, mu, time_explosion) + + # Perform required assertions + assert_almost_equal(obtained, expected) + + +def test_get_random_mu(): + """ + Ensure that different calls results + """ + output1 = r_packet.get_random_mu() + output2 = r_packet.get_random_mu() + assert output1 != output2 + +@pytest.mark.parametrize( + ['current_shell_id', 'delta_shell', 'no_of_shells'], + [(132, 11, 132), + (132, 1, 133), + (132, 2, 133)] +) +def test_move_packet_across_shell_boundary_emitted(packet, current_shell_id, + delta_shell, + no_of_shells): + packet.current_shell_id = current_shell_id + r_packet.move_packet_across_shell_boundary(packet, delta_shell, + no_of_shells) + assert packet.status == r_packet.PacketStatus.EMITTED + +@pytest.mark.parametrize( + ['current_shell_id', 'delta_shell', 'no_of_shells'], + [(132, -133, 132), + (132, -133, 133), + (132, -1e9, 133)] +) +def test_move_packet_across_shell_boundary_reabsorbed(packet, current_shell_id, + delta_shell, + no_of_shells): + packet.current_shell_id = current_shell_id + r_packet.move_packet_across_shell_boundary(packet, delta_shell, + no_of_shells) + assert packet.status == r_packet.PacketStatus.REABSORBED + + +@pytest.mark.parametrize( + ['current_shell_id', 'delta_shell', 'no_of_shells'], + [(132, -1, 199), + (132, 0, 132), + (132, 20, 154)] +) +def test_move_packet_across_shell_boundary_increment(packet, current_shell_id, + delta_shell, + no_of_shells): + packet.current_shell_id = current_shell_id + r_packet.move_packet_across_shell_boundary(packet, delta_shell, + no_of_shells) + assert packet.current_shell_id == current_shell_id + delta_shell + +#SAYS WE NEED TO FIX/REMOVE PACKET CALC ENERGY BY MOVING DOPPLER FACTOR TO FUNCTION +""" +@pytest.mark.parametrize( + ['distance_trace', 'time_explosion'], + [(0, 1), + (1, 1), + (1, 1e9)] +) +def test_packet_energy_limit_one(packet, distance_trace, time_explosion): + initial_energy = packet.energy + new_energy = r_packet.calc_packet_energy(packet, distance_trace, time_explosion) + assert_almost_equal(new_energy, initial_energy) +""" + +@pytest.mark.parametrize( + ['packet_params', 'expected_params'], + [({'mu': 0.3, 'r': 7.5e14}, + {'d_boundary': 259376919351035.88}), + + ({'mu': -.3, 'r': 7.5e13}, + {'d_boundary': -664987228972291.5}), + + ({'mu': -.3, 'r': 7.5e14}, + {'d_boundary': 709376919351035.9})] +) +def test_calculate_distance_boundary(packet_params, expected_params, model): + mu = packet_params['mu'] + r = packet_params['r'] + + d_boundary = r_packet.calculate_distance_boundary( + r, mu, model.r_inner[0], model.r_outer[0]) + + #Accuracy to within 0.1cm + assert_almost_equal(d_boundary[0], expected_params['d_boundary'], decimal=1) +# +# +# TODO: split this into two tests - one to assert errors and other for d_line +@pytest.mark.parametrize( + ['packet_params', 'expected_params'], + [({'nu_line': 0.1, 'next_line_id': 0, 'is_last_line': True}, + {'tardis_error': None, 'd_line': 1e+99}), + + ({'nu_line': 0.2, 'next_line_id': 1, 'is_last_line': False}, + {'tardis_error': None, 'd_line': 7.792353908000001e+17}), + + ({'nu_line': 0.5, 'next_line_id': 1, 'is_last_line': False}, + {'tardis_error': r_packet.MonteCarloException, 'd_line': 0.0}), + + ({'nu_line': 0.6, 'next_line_id': 0, 'is_last_line': False}, + {'tardis_error': r_packet.MonteCarloException, 'd_line': 0.0})] +) +def test_calculate_distance_line(packet_params, expected_params, packet, model): + nu_line = packet_params['nu_line'] + is_last_line = packet_params['is_last_line'] + + time_explosion = model.time_explosion + + doppler_factor = r_packet.get_doppler_factor(packet.r, + packet.mu, + time_explosion) + comov_nu = packet.nu * doppler_factor + + d_line = 0 + obtained_tardis_error = None + try: + d_line = r_packet.calculate_distance_line(packet, + comov_nu, + is_last_line, + nu_line, + time_explosion) + except r_packet.MonteCarloException: + obtained_tardis_error = r_packet.MonteCarloException + + assert_almost_equal(d_line, expected_params['d_line']) + assert obtained_tardis_error == expected_params['tardis_error'] \ No newline at end of file diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_r_packet.py b/tardis/montecarlo/montecarlo_numba/tests/test_r_packet.py deleted file mode 100644 index e69de29bb2d..00000000000 From fda42500955e3fb093b455f172ab650bfeb1728f Mon Sep 17 00:00:00 2001 From: Andrew Fullard Date: Wed, 14 Oct 2020 13:58:21 -0400 Subject: [PATCH 16/21] Added empty tests for packet-related functions --- .../montecarlo_numba/tests/test_packet.py | 151 ++++++++++-------- 1 file changed, 88 insertions(+), 63 deletions(-) diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py index d37702e0509..47fbdc56ddc 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py @@ -39,6 +39,74 @@ def model(): ) +@pytest.mark.parametrize( + ['packet_params', 'expected_params'], + [({'mu': 0.3, 'r': 7.5e14}, + {'d_boundary': 259376919351035.88}), + + ({'mu': -.3, 'r': 7.5e13}, + {'d_boundary': -664987228972291.5}), + + ({'mu': -.3, 'r': 7.5e14}, + {'d_boundary': 709376919351035.9})] +) +def test_calculate_distance_boundary(packet_params, expected_params, model): + mu = packet_params['mu'] + r = packet_params['r'] + + d_boundary = r_packet.calculate_distance_boundary( + r, mu, model.r_inner[0], model.r_outer[0]) + + #Accuracy to within 0.1cm + assert_almost_equal(d_boundary[0], expected_params['d_boundary'], decimal=1) +# +# +# TODO: split this into two tests - one to assert errors and other for d_line +@pytest.mark.parametrize( + ['packet_params', 'expected_params'], + [({'nu_line': 0.1, 'next_line_id': 0, 'is_last_line': True}, + {'tardis_error': None, 'd_line': 1e+99}), + + ({'nu_line': 0.2, 'next_line_id': 1, 'is_last_line': False}, + {'tardis_error': None, 'd_line': 7.792353908000001e+17}), + + ({'nu_line': 0.5, 'next_line_id': 1, 'is_last_line': False}, + {'tardis_error': r_packet.MonteCarloException, 'd_line': 0.0}), + + ({'nu_line': 0.6, 'next_line_id': 0, 'is_last_line': False}, + {'tardis_error': r_packet.MonteCarloException, 'd_line': 0.0})] +) +def test_calculate_distance_line(packet_params, expected_params, packet, model): + nu_line = packet_params['nu_line'] + is_last_line = packet_params['is_last_line'] + + time_explosion = model.time_explosion + + doppler_factor = r_packet.get_doppler_factor(packet.r, + packet.mu, + time_explosion) + comov_nu = packet.nu * doppler_factor + + d_line = 0 + obtained_tardis_error = None + try: + d_line = r_packet.calculate_distance_line(packet, + comov_nu, + is_last_line, + nu_line, + time_explosion) + except r_packet.MonteCarloException: + obtained_tardis_error = r_packet.MonteCarloException + + assert_almost_equal(d_line, expected_params['d_line']) + assert obtained_tardis_error == expected_params['tardis_error'] + +def test_calculate_distance_electron(): + pass + +def test_calculate_tau_electron(): + pass + @pytest.mark.parametrize( ['mu', 'r', 'inv_t_exp', 'expected'], [(0.3, 7.5e14, 1 / 5.2e7, 0.9998556693818854), @@ -59,7 +127,6 @@ def test_get_doppler_factor(mu, r, inv_t_exp, expected): # Perform required assertions assert_almost_equal(obtained, expected) - @pytest.mark.parametrize( ['mu', 'r', 'inv_t_exp', 'expected'], [(0.3, 7.5e14, 1 / 5.2e7, 1/0.9998556693818854), @@ -88,6 +155,24 @@ def test_get_random_mu(): output2 = r_packet.get_random_mu() assert output1 != output2 +def test_update_line_estimators(): + pass + +def test_trace_packet(): + pass + +def test_move_r_packet(): + pass + +def test_set_estimators(): + pass + +def test_set_estimators_full_relativity(): + pass + +def test_line_emission(): + pass + @pytest.mark.parametrize( ['current_shell_id', 'delta_shell', 'no_of_shells'], [(132, 11, 132), @@ -144,65 +229,5 @@ def test_packet_energy_limit_one(packet, distance_trace, time_explosion): new_energy = r_packet.calc_packet_energy(packet, distance_trace, time_explosion) assert_almost_equal(new_energy, initial_energy) """ - -@pytest.mark.parametrize( - ['packet_params', 'expected_params'], - [({'mu': 0.3, 'r': 7.5e14}, - {'d_boundary': 259376919351035.88}), - - ({'mu': -.3, 'r': 7.5e13}, - {'d_boundary': -664987228972291.5}), - - ({'mu': -.3, 'r': 7.5e14}, - {'d_boundary': 709376919351035.9})] -) -def test_calculate_distance_boundary(packet_params, expected_params, model): - mu = packet_params['mu'] - r = packet_params['r'] - - d_boundary = r_packet.calculate_distance_boundary( - r, mu, model.r_inner[0], model.r_outer[0]) - - #Accuracy to within 0.1cm - assert_almost_equal(d_boundary[0], expected_params['d_boundary'], decimal=1) -# -# -# TODO: split this into two tests - one to assert errors and other for d_line -@pytest.mark.parametrize( - ['packet_params', 'expected_params'], - [({'nu_line': 0.1, 'next_line_id': 0, 'is_last_line': True}, - {'tardis_error': None, 'd_line': 1e+99}), - - ({'nu_line': 0.2, 'next_line_id': 1, 'is_last_line': False}, - {'tardis_error': None, 'd_line': 7.792353908000001e+17}), - - ({'nu_line': 0.5, 'next_line_id': 1, 'is_last_line': False}, - {'tardis_error': r_packet.MonteCarloException, 'd_line': 0.0}), - - ({'nu_line': 0.6, 'next_line_id': 0, 'is_last_line': False}, - {'tardis_error': r_packet.MonteCarloException, 'd_line': 0.0})] -) -def test_calculate_distance_line(packet_params, expected_params, packet, model): - nu_line = packet_params['nu_line'] - is_last_line = packet_params['is_last_line'] - - time_explosion = model.time_explosion - - doppler_factor = r_packet.get_doppler_factor(packet.r, - packet.mu, - time_explosion) - comov_nu = packet.nu * doppler_factor - - d_line = 0 - obtained_tardis_error = None - try: - d_line = r_packet.calculate_distance_line(packet, - comov_nu, - is_last_line, - nu_line, - time_explosion) - except r_packet.MonteCarloException: - obtained_tardis_error = r_packet.MonteCarloException - - assert_almost_equal(d_line, expected_params['d_line']) - assert obtained_tardis_error == expected_params['tardis_error'] \ No newline at end of file +def test_test_for_close_line(): + pass \ No newline at end of file From f23529f6c48ec2a8864e719e92719fe58f02f219 Mon Sep 17 00:00:00 2001 From: Andrew Fullard Date: Wed, 14 Oct 2020 14:57:09 -0400 Subject: [PATCH 17/21] Added update line estimator test --- .../montecarlo_numba/tests/test_packet.py | 28 +++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py index 47fbdc56ddc..0e9113b71d0 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py @@ -38,6 +38,14 @@ def model(): time_explosion = 5.2e7 ) +@pytest.fixture(scope="function") +def estimators(): + return numba_interface.Estimators( + j_estimator = np.array([0.0, 0.0], dtype=np.float64), + nu_bar_estimator = np.array([0.0, 0.0], dtype=np.float64), + j_blue_estimator = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], dtype=np.float64), + Edotlu_estimator = np.array([[0.0, 0.0, 1.0], [0.0, 0.0, 1.0]], dtype=np.float64) + ) @pytest.mark.parametrize( ['packet_params', 'expected_params'], @@ -155,8 +163,24 @@ def test_get_random_mu(): output2 = r_packet.get_random_mu() assert output1 != output2 -def test_update_line_estimators(): - pass + +@pytest.mark.parametrize( + ['cur_line_id', 'distance_trace', 'time_explosion', + 'expected_j_blue', 'expected_Edotlu'], + [(0, 1e12, 5.2e7, + [[2.249673812803061, 0.0, 0.0], [0.0, 0.0, 0.0]], [[2.249673812803061*0.4, 0.0, 1.0], [0.0, 0.0, 1.0]]), + (0, 0, 5.2e7, + [[2.249675256109242, 0.0, 0.0], [0.0, 0.0, 0.0]], [[2.249675256109242*0.4, 0.0, 1.0], [0.0, 0.0, 1.0],]), + (1, 1e5, 1e10, + [[0.0, 0.0, 0.0], [2.249998311331767, 0.0, 0.0]], [[0.0, 0.0, 1.0], [2.249998311331767*0.4, 0.0, 1.0]])] +) +def test_update_line_estimators(estimators, packet, cur_line_id, distance_trace, + time_explosion, expected_j_blue, expected_Edotlu): + r_packet.update_line_estimators(estimators, packet, cur_line_id, distance_trace, + time_explosion) + + assert_allclose(estimators.j_blue_estimator, expected_j_blue) + assert_allclose(estimators.Edotlu_estimator, expected_Edotlu) def test_trace_packet(): pass From ef60d7835b2158a848e7514a29c5a81217eb3249 Mon Sep 17 00:00:00 2001 From: Andrew Fullard Date: Wed, 14 Oct 2020 14:57:35 -0400 Subject: [PATCH 18/21] Nicer formatting --- .../montecarlo/montecarlo_numba/tests/test_packet.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py index 0e9113b71d0..bc5862578a9 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py @@ -166,13 +166,17 @@ def test_get_random_mu(): @pytest.mark.parametrize( ['cur_line_id', 'distance_trace', 'time_explosion', - 'expected_j_blue', 'expected_Edotlu'], + 'expected_j_blue', + 'expected_Edotlu'], [(0, 1e12, 5.2e7, - [[2.249673812803061, 0.0, 0.0], [0.0, 0.0, 0.0]], [[2.249673812803061*0.4, 0.0, 1.0], [0.0, 0.0, 1.0]]), + [[2.249673812803061, 0.0, 0.0], [0.0, 0.0, 0.0]], + [[2.249673812803061*0.4, 0.0, 1.0], [0.0, 0.0, 1.0]]), (0, 0, 5.2e7, - [[2.249675256109242, 0.0, 0.0], [0.0, 0.0, 0.0]], [[2.249675256109242*0.4, 0.0, 1.0], [0.0, 0.0, 1.0],]), + [[2.249675256109242, 0.0, 0.0], [0.0, 0.0, 0.0]], + [[2.249675256109242*0.4, 0.0, 1.0], [0.0, 0.0, 1.0],]), (1, 1e5, 1e10, - [[0.0, 0.0, 0.0], [2.249998311331767, 0.0, 0.0]], [[0.0, 0.0, 1.0], [2.249998311331767*0.4, 0.0, 1.0]])] + [[0.0, 0.0, 0.0], [2.249998311331767, 0.0, 0.0]], + [[0.0, 0.0, 1.0], [2.249998311331767*0.4, 0.0, 1.0]])] ) def test_update_line_estimators(estimators, packet, cur_line_id, distance_trace, time_explosion, expected_j_blue, expected_Edotlu): From d8a1b2e6470e894a71aa3c180a14d71b5449acc7 Mon Sep 17 00:00:00 2001 From: Andrew Fullard Date: Wed, 14 Oct 2020 15:17:19 -0400 Subject: [PATCH 19/21] Added empty tests that assert False --- tardis/montecarlo/montecarlo_numba/tests/test_base.py | 5 +++++ .../montecarlo_numba/tests/test_interaction.py | 8 ++++++++ .../montecarlo/montecarlo_numba/tests/test_macro_atom.py | 2 ++ .../montecarlo_numba/tests/test_numba_interface.py | 5 +++++ .../montecarlo_numba/tests/test_single_packet_loop.py | 8 ++++++++ tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py | 9 +++++++++ 6 files changed, 37 insertions(+) diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_base.py b/tardis/montecarlo/montecarlo_numba/tests/test_base.py index e69de29bb2d..a86308c18c3 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_base.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_base.py @@ -0,0 +1,5 @@ +def test_montecarlo_radial1d(): + assert False + +def test_montecarlo_main_loop(): + assert False \ No newline at end of file diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_interaction.py b/tardis/montecarlo/montecarlo_numba/tests/test_interaction.py index e69de29bb2d..2f977553187 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_interaction.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_interaction.py @@ -0,0 +1,8 @@ +def test_thomson_scatter(): + assert False + +def test_line_scatter(): + assert False + +def test_line_emission(): + assert False \ No newline at end of file diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_macro_atom.py b/tardis/montecarlo/montecarlo_numba/tests/test_macro_atom.py index e69de29bb2d..02d56da7aaa 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_macro_atom.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_macro_atom.py @@ -0,0 +1,2 @@ +def test_macro_atom(): + assert False \ No newline at end of file diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py b/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py index e69de29bb2d..5741b585559 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py @@ -0,0 +1,5 @@ +def test_numba_plasma_initialize(): + assert False + +def test_configuration_initialize(): + assert False \ No newline at end of file diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_single_packet_loop.py b/tardis/montecarlo/montecarlo_numba/tests/test_single_packet_loop.py index e69de29bb2d..71cb081c801 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_single_packet_loop.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_single_packet_loop.py @@ -0,0 +1,8 @@ +def test_single_packet_loop(): + assert False + +def test_set_packet_props_partial_relativity(): + assert False + +def test_set_packet_props_full_relativity(): + assert False \ No newline at end of file diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py b/tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py index e69de29bb2d..3c902083e07 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py @@ -0,0 +1,9 @@ +def test_trace_vpacket_within_shell(): + assert False + +def test_trace_vpacket(): + assert False + +def test_trace_vpacket_volley(): + assert False + From 8cf526c8d2c44c5285c860f948bcdc28be22943d Mon Sep 17 00:00:00 2001 From: rodot- Date: Wed, 14 Oct 2020 16:52:17 -0400 Subject: [PATCH 20/21] Wrote unit tests for r_packet.move_r_packet. Should be noted, we found that when updating global variables use in numba functions, the numba functions must be recompiled in order for the changes to take effect. Keep this in mind when testing. The cases for full relativity being enabled SHOULD FAIL, there is a discrepancy between how the code is written and what the tests check. Someone who is more confident in their knowledge of the physics than me should check it out. (i.e. should we be off by a factor of doppler factor, or doppler factor squared?) --- .../montecarlo_numba/tests/test_packet.py | 48 +++++++++++++++++-- 1 file changed, 45 insertions(+), 3 deletions(-) diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py index bc5862578a9..9ac37fc3755 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py @@ -8,6 +8,7 @@ import tardis.montecarlo.montecarlo_numba.numba_interface as numba_interface from tardis import constants as const from tardis.montecarlo.montecarlo_numba.numba_interface import Estimators +import tardis.montecarlo.montecarlo_numba.numba_config as numba_config from tardis.montecarlo.montecarlo_numba import macro_atom C_SPEED_OF_LIGHT = const.c.to('cm/s').value @@ -189,8 +190,49 @@ def test_update_line_estimators(estimators, packet, cur_line_id, distance_trace, def test_trace_packet(): pass -def test_move_r_packet(): - pass +@pytest.mark.parametrize('ENABLE_FULL_RELATIVITY', [True, False]) +@pytest.mark.parametrize( + ['packet_params', 'expected_params'], + [({'nu': 0.4, 'mu': 0.3, 'energy': 0.9, 'r': 7.5e14}, + {'mu': 0.3120599529139568, 'r': 753060422542573.9, + 'j': 8998701024436.969, 'nubar': 3598960894542.354}), + + ({'nu': 0.6, 'mu': -.5, 'energy': 0.5, 'r': 8.1e14}, + {'mu': -.4906548373534084, 'r': 805046582503149.2, + 'j': 5001298975563.031, 'nubar': 3001558973156.1387})] +) +def test_move_r_packet(packet_params, expected_params, packet, model, estimators, ENABLE_FULL_RELATIVITY): + distance = 1.e13 + packet.nu = packet_params['nu'] + packet.mu = packet_params['mu'] + packet.energy = packet_params['energy'] + packet.r = packet_params['r'] + + numba_config.ENABLE_FULL_RELATIVITY = ENABLE_FULL_RELATIVITY + r_packet.move_r_packet.recompile() # This must be done as move_r_packet was jitted with ENABLE_FULL_RELATIVITY + doppler_factor = r_packet.get_doppler_factor( + packet.r, + packet.mu, + model.time_explosion) + + r_packet.move_r_packet( + packet, distance, model.time_explosion, estimators) + + assert_almost_equal(packet.mu, expected_params['mu']) + assert_almost_equal(packet.r, expected_params['r']) + + expected_j = expected_params['j'] + expected_nubar = expected_params['nubar'] + + if ENABLE_FULL_RELATIVITY: + expected_j *= doppler_factor + expected_nubar *= doppler_factor + + numba_config.ENABLE_FULL_RELATIVITY = False + assert_allclose(estimators.j_estimator[packet.current_shell_id], + expected_j, rtol=5e-7) + assert_allclose(estimators.nu_bar_estimator[packet.current_shell_id], + expected_nubar, rtol=5e-7) def test_set_estimators(): pass @@ -258,4 +300,4 @@ def test_packet_energy_limit_one(packet, distance_trace, time_explosion): assert_almost_equal(new_energy, initial_energy) """ def test_test_for_close_line(): - pass \ No newline at end of file + pass From 89ced056af2005b9975b98b404bcaa24523bf161 Mon Sep 17 00:00:00 2001 From: Andrew Fullard Date: Wed, 14 Oct 2020 16:55:38 -0400 Subject: [PATCH 21/21] Bad vpacket test --- .../montecarlo_numba/tests/test_vpacket.py | 68 ++++++++++++++++++- 1 file changed, 67 insertions(+), 1 deletion(-) diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py b/tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py index 3c902083e07..dd7cdb42a9d 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py @@ -1,4 +1,70 @@ -def test_trace_vpacket_within_shell(): +import os +import pytest +import numpy as np +import pandas as pd +import tardis.montecarlo.formal_integral as formal_integral +import tardis.montecarlo.montecarlo_numba.r_packet as r_packet +import tardis.montecarlo.montecarlo_numba.vpacket as vpacket +import tardis.montecarlo.montecarlo_configuration as mc +import tardis.montecarlo.montecarlo_numba.numba_interface as numba_interface +from tardis import constants as const +from tardis.montecarlo.montecarlo_numba.numba_interface import Estimators +from tardis.montecarlo.montecarlo_numba import macro_atom +C_SPEED_OF_LIGHT = const.c.to('cm/s').value + +from numpy.testing import ( + assert_equal, + assert_almost_equal, + assert_array_equal, + assert_allclose + ) + + +@pytest.fixture(scope="function") +def v_packet(): + return vpacket.VPacket( + r = 7.5e14, + nu = 0.4, + mu = 0.3, + energy = 0.9, + current_shell_id = 0, + next_line_id = 0, + index = 0, + is_close_line = 0 + ) + +@pytest.fixture(scope="function") +def model(): + return numba_interface.NumbaModel( + r_inner = np.array([6.912e14, 8.64e14], dtype=np.float64), + r_outer = np.array([8.64e14, 1.0368e15], dtype=np.float64), + time_explosion = 5.2e7 + ) + +@pytest.fixture(scope="function") +def plasma(): + return numba_interface.NumbaPlasma( + electron_density=1.0e9*np.ones(2, dtype=np.float64), + line_list_nu=np.array([ + 1.26318289e+16, + 1.26318289e+16, + 1.23357675e+16, + 1.23357675e+16, + 1.16961598e+16], dtype=np.float64), + tau_sobolev=np.ones((2, 1000), dtype=np.float64), + transition_probabilities=np.zeros((2, 2), dtype=np.float64), + line2macro_level_upper=np.zeros(2, dtype=np.int64), + macro_block_references=np.zeros(2, dtype=np.int64), + transition_type=np.zeros(2, dtype=np.int64), + destination_level_id=np.zeros(2, dtype=np.int64), + transition_line_id=np.zeros(2, dtype=np.int64) + ) + +def test_trace_vpacket_within_shell(v_packet, model, plasma): + tau_trace_combined, distance_boundary, delta_shell = vpacket.trace_vpacket_within_shell( + v_packet, + model, + plasma) assert False def test_trace_vpacket():