diff --git a/tardis/montecarlo/formal_integral.py b/tardis/montecarlo/formal_integral.py index fe9f82a90de..087d1b9ddf8 100644 --- a/tardis/montecarlo/formal_integral.py +++ b/tardis/montecarlo/formal_integral.py @@ -6,9 +6,12 @@ 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 +from tardis.montecarlo.montecarlo_numba.numba_interface \ + import numba_plasma_initialize from tardis.montecarlo.montecarlo import formal_integral from tardis.montecarlo.spectrum import TARDISSpectrum @@ -17,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 @@ -33,7 +37,10 @@ class FormalIntegrator(object): def __init__(self, model, plasma, runner, points=1000): self.model = model - self.plasma = plasma + if plasma: + self.plasma = numba_plasma_initialize(plasma) + self.atomic_data = plasma.atomic_data + self.original_plasma = plasma self.runner = runner self.points = points @@ -81,8 +88,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' @@ -119,25 +125,24 @@ def make_source_function(self): """ model = self.model - plasma = self.plasma runner = self.runner - atomic_data = self.plasma.atomic_data - macro_ref = atomic_data.macro_atom_references - macro_data = atomic_data.macro_atom_data - no_lvls = len(atomic_data.levels) + macro_ref = self.atomic_data.macro_atom_references + macro_data = self.atomic_data.macro_atom_data + + 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 = 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 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 @@ -149,7 +154,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 @@ -168,13 +173,14 @@ 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 = atomic_data.macro_atom_data[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 = plasma.transition_probabilities[(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 = 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] @@ -185,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(-plasma.tau_sobolevs.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 = plasma.tau_sobolevs.values - runner.electron_densities_integ = 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 @@ -213,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) @@ -242,10 +248,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,8 +269,8 @@ 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 @@ -278,19 +280,17 @@ def _formal_integral(self, iT, inu, inu_size, att_S_ul, Jred_lu, Jblue_lu, N): exp_tau = np.zeros(size_tau) # TODO: multiprocessing offset = 0 - i = 0 size_z = 0 + z = np.zeros(2 * self.model.no_of_shells) 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 + 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 @@ -313,38 +313,39 @@ 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.plasma.line_list_nu, + nu_start, size_line, idx_nu_start) offset = shell_id[0] * size_line # start tracking accumulated e-scattering optical depth 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 @@ -352,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]) @@ -405,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) @@ -452,67 +454,93 @@ 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 - @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): diff --git a/tardis/montecarlo/montecarlo_numba/r_packet.py b/tardis/montecarlo/montecarlo_numba/r_packet.py index 1d9bfba5c6d..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, @@ -171,6 +171,7 @@ def get_doppler_factor_full_relativity(mu, beta): return (1.0 - mu * beta) / math.sqrt(1 - beta * beta) + @njit(**njit_dict) def get_inverse_doppler_factor(r, mu, time_explosion): inv_c = 1 / C_SPEED_OF_LIGHT @@ -189,6 +190,7 @@ def get_inverse_doppler_factor_partial_relativity(mu, beta): def get_inverse_doppler_factor_full_relativity(mu, beta): return (1.0 + mu * beta) / math.sqrt(1 - beta * beta) + @njit(**njit_dict) def get_random_mu(): return 2.0 * np.random.random() - 1.0 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..a86308c18c3 --- /dev/null +++ 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 new file mode 100644 index 00000000000..2f977553187 --- /dev/null +++ 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 new file mode 100644 index 00000000000..02d56da7aaa --- /dev/null +++ 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_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..5741b585559 --- /dev/null +++ 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_packet.py b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py new file mode 100644 index 00000000000..9ac37fc3755 --- /dev/null +++ b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py @@ -0,0 +1,303 @@ +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 +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 + +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.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'], + [({'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), + (-.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( + ['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 + +@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 + +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), + (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) +""" +def test_test_for_close_line(): + pass 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..71cb081c801 --- /dev/null +++ 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 new file mode 100644 index 00000000000..dd7cdb42a9d --- /dev/null +++ b/tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py @@ -0,0 +1,75 @@ +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(): + assert False + +def test_trace_vpacket_volley(): + assert False + 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 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..2ca4c6397c5 100644 --- a/tardis/montecarlo/tests/test_montecarlo.py +++ b/tardis/montecarlo/tests/test_montecarlo.py @@ -0,0 +1,954 @@ +""" +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 +import tardis.montecarlo.montecarlo_configuration as mc +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 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) + 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, 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, 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_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) + +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)] +) +def test_angle_ab_LF_to_CMF_diverge(mu, r, time_explosion, packet): + """ + This equation should diverage as costheta --> 1 and beta --> 1 + """ + packet = r_packet.RPacket(packet.r, packet.mu, packet.nu, packet.energy) + 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( + ['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( + ['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, 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, 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 + + mc.full_relativity = False + + assert_allclose(numba_estimator.j_estimator[packet.current_shell_id], + expected_j, rtol=5e-7) + assert_allclose(numba_estimator.nu_bar_estimator[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(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) + + + +@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( + ['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(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(packet, numba_plasma) + 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 pieces 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) + + + +""" +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( + ['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(packet, mu, r, + inv_t_exp, full_relativity): + packet = r_packet.RPacket(r=r, mu=mu, energy=packet.energy, nu=packet.nu) + mc.full_relativity = bool(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(packet, 1/inv_t_exp, packet.mu) + + 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) + + +@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(packet, model, + mu, r, inv_t_exp): + packet = r_packet.RPacket(r, mu, packet.nu, packet.energy) + model.full_relativity = 1 + + 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) + + +@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(mu, r, t_exp, nu, nu_line, + 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) + 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, numba_estimator) + + 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) diff --git a/tardis/tests/test_tardis_full_formal_integral.py b/tardis/tests/test_tardis_full_formal_integral.py index 53a7343e7ea..67a40819763 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]