diff --git a/tardis/montecarlo/setup_package.py b/tardis/montecarlo/setup_package.py index 6de507d283e..2d462f3199a 100644 --- a/tardis/montecarlo/setup_package.py +++ b/tardis/montecarlo/setup_package.py @@ -22,27 +22,13 @@ def get_extensions(): sources = ['tardis/montecarlo/montecarlo.pyx'] sources += [os.path.relpath(fname) for fname in glob( - os.path.join(os.path.dirname(__file__), 'src', '*.c')) - if not os.path.basename(fname).startswith('test')] + os.path.join(os.path.dirname(__file__), 'src', '*.c'))] sources += [os.path.relpath(fname) for fname in glob( os.path.join(os.path.dirname(__file__), 'src/randomkit', '*.c'))] deps = [os.path.relpath(fname) for fname in glob( os.path.join(os.path.dirname(__file__), 'src', '*.h'))] deps += [os.path.relpath(fname) for fname in glob( os.path.join(os.path.dirname(__file__), 'src/randomkit', '*.h'))] - #tests = [os.path.relpath(fname) for fname in os.wakk( - # os.path.join(os.path.dirname(__file__), 'src', '*.c')) - # if os.path.basename(fname).startswith('test')] - test_sources = [ - os.path.relpath(os.path.join(r,fname)) for r,p,fnames in os.walk(os.path.dirname(__file__)) - for fname in fnames - if os.path.basename(fname).endswith('.c') - ] - test_deps = [ - os.path.relpath(os.path.join(r,fname)) for r,p,fnames in os.walk(os.path.dirname(__file__)) - for fname in fnames - if os.path.basename(fname).endswith('.h') - ] return cythonize( Extension('tardis.montecarlo.montecarlo', sources, @@ -53,19 +39,8 @@ def get_extensions(): extra_compile_args=compile_args, extra_link_args=link_args, define_macros=define_macros) - ) + [ -# Test extension - Extension('tardis.montecarlo.test_montecarlo', test_sources, - include_dirs=['tardis/montecarlo/src', - 'tardis/montecarlo/src/randomkit', - 'numpy'], - depends=test_deps, - extra_compile_args=compile_args, - extra_link_args=link_args, -# library_dirs=['tardis/montecarlo'], -# libraries= ['montecarlo'], - define_macros=define_macros), - ] + ) + def get_package_data(): return {'tardis.montecarlo.tests':['data/*.npy']} diff --git a/tardis/montecarlo/src/test_cmontecarlo.c b/tardis/montecarlo/src/test_cmontecarlo.c deleted file mode 100644 index 079e13eb927..00000000000 --- a/tardis/montecarlo/src/test_cmontecarlo.c +++ /dev/null @@ -1,429 +0,0 @@ -#include -#include -#include -#include -#include - -#include "cmontecarlo.h" - - -rpacket_t * rp; -storage_model_t * sm; - -double TIME_EXPLOSION = 5.2e7; /* 10 days(in seconds) ~ 51840000.0 */ -double R_INNER_VALUE = 6.2e11; /* 12,000xTIME_EXPLOSION ~ 622080000000.0 */ -double R_OUTER_VALUE = 7.8e12; /* 15,000xTIME_EXPLOSION ~ 7776600000000.0 */ -double SIGMA_THOMSON = 6.652486e-25; - -double test_compute_distance2boundary(void); -double test_compute_distance2line(void); -double test_compute_distance2continuum(void); -double test_rpacket_doppler_factor(void); -//double test_move_packet(void); -bool test_move_packet_across_shell_boundary(void); -int64_t test_montecarlo_one_packet(void); -int64_t test_montecarlo_one_packet_loop(void); -bool test_montecarlo_line_scatter(void); -double test_increment_j_blue_estimator(void); -bool test_montecarlo_thomson_scatter(void); -//bool test_macro_atom(void); -double test_calculate_chi_bf(void); -bool test_montecarlo_bound_free_scatter(void); -double test_bf_cross_section(void); -int64_t test_montecarlo_free_free_scatter(void); - -/* initialise RPacket */ -static void init_rpacket(rpacket_t *rp){ - double MU = 0.3; - double R = 7.5e14; - double ENERGY = 0.9; - int NEXT_LINE_ID = 1; - double NU = 0.4; - double NU_LINE = 0.2; - int CURRENT_SHELL_ID = 0; - - double TAU_EVENT = 2.9e13; - - rpacket_set_current_shell_id(rp, CURRENT_SHELL_ID); - rpacket_set_next_shell_id(rp, CURRENT_SHELL_ID+1); - rpacket_set_mu(rp, MU); - rpacket_set_nu(rp, NU); - rpacket_set_r(rp, R); - rpacket_set_last_line(rp, false); - rpacket_set_close_line(rp, false); - rpacket_set_nu_line(rp, NU_LINE); - - rpacket_set_next_line_id(rp, NEXT_LINE_ID); - - rpacket_set_tau_event(rp, TAU_EVENT); - rpacket_set_virtual_packet(rp, 0); - rpacket_set_energy(rp, ENERGY); - rpacket_set_virtual_packet_flag(rp, true); - rpacket_set_status(rp, TARDIS_PACKET_STATUS_IN_PROCESS); - rpacket_set_id(rp, 0); - - rpacket_set_current_continuum_id(rp, 1); -} - -/* initialise storage model */ -static void init_storage_model(storage_model_t *sm){ - int NUMBER_OF_SHELLS = 2; - - sm->time_explosion = TIME_EXPLOSION; - sm->inverse_time_explosion = 1.0/TIME_EXPLOSION; - sm->inverse_sigma_thomson = 1.0/SIGMA_THOMSON; - - /* R_OUTER = {8.64e14, 1.0368e15} */ - sm->r_outer = (double *) malloc(sizeof(double)*NUMBER_OF_SHELLS); - sm->r_outer[0] = 8.64e14; - sm->r_outer[1] = 1.0368e15; - - /* R_INNER = {6.912e14, 8.64e14} */ - sm->r_inner = (double *) malloc(sizeof(double)*NUMBER_OF_SHELLS); - sm->r_inner[0] = 6.912e14; - sm->r_inner[1] = 8.64e14; - sm->no_of_lines = 2; - - /* LINE_LIST_NU = { 1.26318289e+16, 1.26318289e+16, - 1.23357675e+16, 1.23357675e+16, - 1.16961598e+16 }; - */ - sm->line_list_nu = (double *) malloc(sizeof(double)*5); - sm->line_list_nu[0] = 1.26318289e+16; - sm->line_list_nu[1] = 1.26318289e+16; - sm->line_list_nu[2] = 1.23357675e+16; - sm->line_list_nu[3] = 1.23357675e+16; - sm->line_list_nu[4] = 1.16961598e+16; - - /* INVERSE_ELECTRON_DENSITIES = {} */ - sm->inverse_electron_densities = (double *) malloc(sizeof(double)*NUMBER_OF_SHELLS); - sm->inverse_electron_densities[0] = 1e-9; - sm->inverse_electron_densities[1] = 1e-9; - sm->electron_densities = (double *) malloc(sizeof(double )*2); - sm->electron_densities[0] = 1e9; - sm->electron_densities[1] = 1e9; - - - sm->js = (double *) malloc(sizeof(double)*NUMBER_OF_SHELLS); - sm->js[0] = 0; - sm->js[1] = 0; - - sm->nubars = (double *) malloc(sizeof(double)*NUMBER_OF_SHELLS); - sm->nubars[0] = 0; - sm->nubars[1] = 0; - - sm->last_line_interaction_in_id = (int64_t *) malloc(sizeof(int64_t)*NUMBER_OF_SHELLS); - sm->last_line_interaction_in_id[0] = 0; - sm->last_line_interaction_in_id[1] = 0; - - sm->last_line_interaction_shell_id = (int64_t *) malloc(sizeof(int64_t)*NUMBER_OF_SHELLS); - sm->last_line_interaction_shell_id[0] = 0; - sm->last_line_interaction_shell_id[1] = 0; - - sm->last_interaction_type = (int64_t *) malloc(sizeof(int64_t)*1); - sm->last_interaction_type[0] = 2; - - sm->line_interaction_id = 0; - - sm->line_lists_j_blues_nd = 0; - - sm->line_lists_j_blues = (double *) malloc(sizeof(double )*2); - sm->line_lists_j_blues[0] = 1e-10; - sm->line_lists_j_blues[1] = 1e-10; - - sm->line_lists_tau_sobolevs = (double *) malloc(sizeof(double )*2); - sm->line_lists_tau_sobolevs[0] = 1e-5; - sm->line_lists_tau_sobolevs[1] = 1000; - - sm->line2macro_level_upper = (int64_t *) malloc(sizeof(int64_t)*2); - sm->line2macro_level_upper[0] = 0; - sm->line2macro_level_upper[1] = 0; - - sm->reflective_inner_boundary = false; - sm->inner_boundary_albedo = 0.0; - sm->no_of_shells = NUMBER_OF_SHELLS; - - sm->spectrum_start_nu = 1.e14; - sm->spectrum_delta_nu = 293796608840.0; - sm->spectrum_end_nu = 6.e15; -/* FIXME: copy paste error below? - sm->spectrum_virt_start_nu = 1.e14; - sm->spectrum_virt_end_nu = 293796608840.0; - sm->spectrum_delta_nu = 6.e15; - */ - - sm->spectrum_virt_start_nu = 1.e14; - sm->spectrum_virt_end_nu = 6.e15; - - sm->spectrum_virt_nu = (double *) calloc(20000,sizeof(double )); - - /* - * Initialising the below values to 0 untill - * I get the real data ! - */ - - sm->t_electrons = (double *) malloc(sizeof(double )*2); - sm->t_electrons[0] = 2; - sm->t_electrons[1] = 2; - - sm->l_pop = (double *) malloc(sizeof(double )*20000); - for (size_t i=0; i<20000; ++i) - sm->l_pop[i]=2; - - sm->l_pop_r = (double *) malloc(sizeof(double )* 20000); - for (size_t i=0; i<20000; ++i) - sm->l_pop_r[i]=3; - - sm->continuum_list_nu = (double *) malloc(sizeof(double )* 20000); - for (size_t i=0; i<20000; ++i) - sm->continuum_list_nu[i]=1e13; - - sm->no_of_edges = 100; - - sm->chi_bf_tmp_partial = (double *) malloc(sizeof(double )* 20000); - for (size_t i=0; i<20000; ++i) - sm->chi_bf_tmp_partial[i]=160; - -} -static void dealloc_storage_model(storage_model_t *sm){ - free(sm->r_outer); - free(sm->r_inner); - free(sm->line_list_nu); - free(sm->inverse_electron_densities); - free(sm->electron_densities); - free(sm->js); - free(sm->nubars); - free(sm->last_line_interaction_in_id); - free(sm->last_line_interaction_shell_id); - free(sm->last_interaction_type); - free(sm->line_lists_j_blues); - free(sm->line_lists_tau_sobolevs); - free(sm->line2macro_level_upper); - free(sm->spectrum_virt_nu); - free(sm->t_electrons); - free(sm->l_pop); - free(sm->l_pop_r); - free(sm->continuum_list_nu); - free(sm->chi_bf_tmp_partial); -} - -static void irandom(rk_state *mt_state) -{ -memset(mt_state,0,sizeof(rk_state)); -} - -double -test_compute_distance2boundary(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - compute_distance2boundary(&rp, &sm); - double D_BOUNDARY = rpacket_get_d_boundary(&rp); - dealloc_storage_model(&sm); - return D_BOUNDARY; -} -double -test_compute_distance2line(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - // FIXME MR: return status of compute_distance2line() is ignored - compute_distance2line(&rp, &sm); - double D_LINE = rpacket_get_d_line(&rp); - dealloc_storage_model(&sm); - return D_LINE; -} - -double -test_compute_distance2continuum(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - compute_distance2continuum(&rp, &sm); - dealloc_storage_model(&sm); - return rp.d_cont; -} - -double -test_rpacket_doppler_factor(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - double res= rpacket_doppler_factor(&rp, &sm); - dealloc_storage_model(&sm); - return res; -} - -//double -//test_move_packet(void){ -// double DISTANCE = 1e13; -// rpacket_t rp; -// storage_model_t sm; -// init_rpacket(&rp); -// init_storage_model(&sm); -// double res = move_packet(&rp, &sm, DISTANCE); -// dealloc_storage_model(&sm); -// return res; -//} - -double -test_increment_j_blue_estimator(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - int64_t j_blue_idx = 0; - compute_distance2line(&rp, &sm); - move_packet(&rp, &sm, 1e13); - double d_line = rpacket_get_d_line(&rp); - increment_j_blue_estimator(&rp, &sm, d_line, j_blue_idx); - double res = sm.line_lists_j_blues[j_blue_idx]; - dealloc_storage_model(&sm); - return res; -} - -bool -test_montecarlo_line_scatter(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - rk_state mt_state; - irandom(&mt_state); - double DISTANCE = 1e13; - montecarlo_line_scatter(&rp, &sm, DISTANCE, &mt_state); - dealloc_storage_model(&sm); - return true; -} - -bool -test_montecarlo_thomson_scatter(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - rk_state mt_state; - irandom(&mt_state); - double DISTANCE = 1e13; - montecarlo_thomson_scatter(&rp, &sm, DISTANCE, &mt_state); - dealloc_storage_model(&sm); - return true; -} - -bool -test_move_packet_across_shell_boundary(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - rk_state mt_state; - irandom(&mt_state); - double DISTANCE = 0.95e13; -// MR: wrong: move_packet_across_shell_boundary() returns void - move_packet_across_shell_boundary(&rp, &sm, DISTANCE, &mt_state); - dealloc_storage_model(&sm); - return true; -} - - -int64_t -test_montecarlo_one_packet(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - rk_state mt_state; - irandom(&mt_state); - int64_t res = montecarlo_one_packet(&sm, &rp, 1, &mt_state); - dealloc_storage_model(&sm); - return res; -} - -int64_t -test_montecarlo_one_packet_loop(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - rk_state mt_state; - irandom(&mt_state); - int64_t res= montecarlo_one_packet_loop(&sm, &rp, 1, &mt_state); - dealloc_storage_model(&sm); - return res; -} - -//bool -//test_macro_atom(void){ -// macro_atom(rp, sm); -// return true; -//} - -double -test_calculate_chi_bf(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - rk_state mt_state; - irandom(&mt_state); - int64_t j_blue_idx = 0; - compute_distance2line(&rp, &sm); - move_packet(&rp, &sm, 1e13); - double d_line = rpacket_get_d_line(&rp); - increment_j_blue_estimator(&rp, &sm, d_line, j_blue_idx); - double DISTANCE = 1e13; - montecarlo_line_scatter(&rp, &sm, DISTANCE, &mt_state); - DISTANCE = 0.95e13; -// MR: wrong: move_packet_across_shell_boundary() returns void - move_packet_across_shell_boundary(&rp, &sm, DISTANCE, &mt_state); - montecarlo_one_packet(&sm, &rp, 1, &mt_state); - montecarlo_one_packet_loop(&sm, &rp, 1, &mt_state); - DISTANCE = 1e13; - montecarlo_thomson_scatter(&rp, &sm, DISTANCE, &mt_state); - calculate_chi_bf(&rp, &sm); - double res = rpacket_doppler_factor (&rp, &sm); - dealloc_storage_model(&sm); - return res; -} - -bool -test_montecarlo_bound_free_scatter(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - rk_state mt_state; - irandom(&mt_state); - double DISTANCE = 1e13; - montecarlo_bound_free_scatter(&rp, &sm, DISTANCE, &mt_state); - dealloc_storage_model(&sm); - return (rpacket_get_status(&rp)!=0); -} - -double -test_bf_cross_section(void){ - storage_model_t sm; - init_storage_model(&sm); - double CONV_MU = 0.4; - double res = bf_cross_section(&sm, 1, CONV_MU); - dealloc_storage_model(&sm); - return res; -} - -int64_t -test_montecarlo_free_free_scatter(void){ - rpacket_t rp; - storage_model_t sm; - init_rpacket(&rp); - init_storage_model(&sm); - rk_state mt_state; - irandom(&mt_state); - double DISTANCE = 1e13; - montecarlo_free_free_scatter(&rp, &sm, DISTANCE,&mt_state); - dealloc_storage_model(&sm); - return rpacket_get_status(&rp); -} diff --git a/tardis/montecarlo/src/test_rpacket.c b/tardis/montecarlo/src/test_rpacket.c deleted file mode 100644 index 60d0fa10499..00000000000 --- a/tardis/montecarlo/src/test_rpacket.c +++ /dev/null @@ -1,166 +0,0 @@ -#include -#include -#include -#include -#include - -#include "rpacket.h" -#include "status.h" - -bool test_rpacket_get_nu(double); -bool test_rpacket_get_mu(double); -bool test_rpacket_get_energy(double); -bool test_rpacket_get_r(double); -bool test_rpacket_get_tau_event(double); -bool test_rpacket_get_nu_line(double); -bool test_rpacket_get_d_boundary(double); -bool test_rpacket_get_d_electron(double); -bool test_rpacket_get_d_line(double); - - -bool test_rpacket_get_current_shell_id(unsigned int); -bool test_rpacket_get_next_line_id(unsigned int); - -bool test_rpacket_get_virtual_packet_flag(int); -bool test_rpacket_get_virtual_packet(int); -bool test_rpacket_get_next_shell_id(int); - -bool test_rpacket_get_last_line(void); -bool test_rpacket_get_close_line(void); -bool test_rpacket_get_status(void); -bool test_rpacket_get_id(void); - -bool -test_rpacket_get_nu(double value){ - rpacket_t rp; - rpacket_set_nu(&rp, value); - return value==rpacket_get_nu(&rp); -} - -bool -test_rpacket_get_mu(double value){ - rpacket_t rp; - rpacket_set_mu(&rp, value); - return value==rpacket_get_mu(&rp); -} - -bool -test_rpacket_get_energy(double value){ - rpacket_t rp; - rpacket_set_energy(&rp, value); - return value==rpacket_get_energy(&rp); -} - -bool -test_rpacket_get_r(double value){ - rpacket_t rp; - rpacket_set_r(&rp, value); - return value==rpacket_get_r(&rp); -} - -bool -test_rpacket_get_tau_event(double value){ - rpacket_t rp; - rpacket_set_tau_event(&rp, value); - return value==rpacket_get_tau_event(&rp); -} - -bool -test_rpacket_get_nu_line(double value){ - rpacket_t rp; - rpacket_set_nu_line(&rp, value); - return value==rpacket_get_nu_line(&rp); -} - -bool -test_rpacket_get_current_shell_id(unsigned int value){ - rpacket_t rp; - rpacket_set_current_shell_id(&rp, value); - return value==rpacket_get_current_shell_id(&rp); -} - -bool -test_rpacket_get_next_line_id(unsigned int value){ - rpacket_t rp; - rpacket_set_next_line_id(&rp, value); - return value==rpacket_get_next_line_id(&rp); -} - -bool -test_rpacket_get_last_line(void){ - rpacket_t rp; - rpacket_set_last_line(&rp, true); - return rpacket_get_last_line(&rp); -} - -bool -test_rpacket_get_close_line(void){ - rpacket_t rp; - rpacket_set_last_line(&rp, true); - return rpacket_get_last_line(&rp); -} - -bool -test_rpacket_get_virtual_packet_flag(int value){ - rpacket_t rp; - rpacket_set_virtual_packet_flag(&rp, value); - return value==rpacket_get_virtual_packet_flag(&rp); -} - -bool -test_rpacket_get_virtual_packet(int value){ - rpacket_t rp; - rpacket_set_virtual_packet(&rp, value); - return value==rpacket_get_virtual_packet(&rp); -} - -bool -test_rpacket_get_d_boundary(double value){ - rpacket_t rp; - rpacket_set_d_boundary(&rp, value); - return value==rpacket_get_d_boundary(&rp); -} - -bool -test_rpacket_get_d_electron(double value){ - rpacket_t rp; - rpacket_set_d_electron(&rp, value); - return value==rpacket_get_d_electron(&rp); -} - -bool -test_rpacket_get_d_line(double value){ - rpacket_t rp; - rpacket_set_d_line(&rp, value); - return value==rpacket_get_d_line(&rp); -} - -bool -test_rpacket_get_next_shell_id(int value){ - rpacket_t rp; - rpacket_set_next_shell_id(&rp, value); - return value==rpacket_get_next_shell_id(&rp); -} - -bool -test_rpacket_get_status(void){ - rpacket_status_t inProcess = TARDIS_PACKET_STATUS_IN_PROCESS; - rpacket_status_t emitted = TARDIS_PACKET_STATUS_EMITTED; - rpacket_status_t reabsorbed = TARDIS_PACKET_STATUS_REABSORBED; - - rpacket_t rp; - rpacket_set_status(&rp, inProcess); - bool res= inProcess==rpacket_get_status(&rp); - rpacket_set_status(&rp, emitted); - res &= emitted==rpacket_get_status(&rp); - rpacket_set_status(&rp, reabsorbed); - res &= reabsorbed==rpacket_get_status(&rp); - return res; -} - -bool -test_rpacket_get_id(void){ - rpacket_t rp; - rpacket_set_id(&rp, 2); - return rpacket_get_id(&rp) == 2; -} diff --git a/tardis/montecarlo/struct.py b/tardis/montecarlo/struct.py index 976d84be834..a68bbfde450 100644 --- a/tardis/montecarlo/struct.py +++ b/tardis/montecarlo/struct.py @@ -99,15 +99,6 @@ class StorageModel(Structure): ('virt_array_size', c_int64) ] - -class RKState(Structure): - _fields_ = [ - ('key', POINTER(c_ulong)), - ('pos', c_int), - ('has_gauss', c_int), - ('gauss', c_double) - ] - # Variables corresponding to `tardis_error_t` enum. TARDIS_ERROR_OK = 0 TARDIS_ERROR_BOUNDS_ERROR = 1 @@ -121,3 +112,15 @@ class RKState(Structure): # Variables corresponding to `ContinuumProcessesStatus` enum. CONTINUUM_OFF = 0 CONTINUUM_ON = 1 + +# `rk_state` specific macros. +RK_STATE_LEN = 624 + + +class RKState(Structure): + _fields_ = [ + ('key', c_ulong * RK_STATE_LEN), + ('pos', c_int), + ('has_gauss', c_int), + ('gauss', c_double) + ] diff --git a/tardis/montecarlo/tests/test_cmontecarlo.py b/tardis/montecarlo/tests/test_cmontecarlo.py index 1c6e4c8a5ac..ed1eacf59bd 100644 --- a/tardis/montecarlo/tests/test_cmontecarlo.py +++ b/tardis/montecarlo/tests/test_cmontecarlo.py @@ -45,8 +45,8 @@ import os import pytest -from ctypes import CDLL, byref, c_uint, c_int64, c_double, c_ulong -from numpy.testing import assert_almost_equal +from ctypes import CDLL, byref, c_uint, c_int64, c_double, c_ulong, POINTER +from numpy.testing import assert_equal, assert_almost_equal from tardis import __path__ as path from tardis.montecarlo.struct import ( @@ -61,11 +61,6 @@ CONTINUUM_ON ) -# Wrap the shared object containing tests for C methods, written in C. -# TODO: Shift all tests here in Python and completely remove this test design. -test_path = os.path.join(path[0], 'montecarlo', 'test_montecarlo.so') -cmontecarlo_tests = CDLL(test_path) - # Wrap the shared object containing C methods, which are tested here. cmontecarlo_filepath = os.path.join(path[0], 'montecarlo', 'montecarlo.so') cmontecarlo_methods = CDLL(cmontecarlo_filepath) @@ -88,8 +83,10 @@ def packet(): current_continuum_id=1, virtual_packet_flag=1, virtual_packet=0, + next_shell_id=1, status=TARDIS_PACKET_STATUS_IN_PROCESS, - id=0 + id=0, + chi_cont=6.652486e-16 ) @@ -165,6 +162,63 @@ def mt_state(): ) +""" +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) + x_insert = c_double(x_insert) + imin = c_int64(imin) + imax = c_int64(imax) + obtained_result = c_int64(0) + + cmontecarlo_methods.reverse_binary_search.restype = c_uint + obtained_tardis_error = cmontecarlo_methods.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(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) + + cmontecarlo_methods.line_search.restype = c_uint + obtained_tardis_error = cmontecarlo_methods.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( ['mu', 'r', 'inv_t_exp', 'expected'], [(0.3, 7.5e14, 1 / 5.2e7, 0.9998556693818854), @@ -253,6 +307,31 @@ def test_compute_distance2continuum(packet_params, expected_params, packet, mode assert_almost_equal(packet.d_cont, expected_params['d_cont']) +@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): + packet.nu = packet_params['nu'] + packet.mu = packet_params['mu'] + packet.energy = packet_params['energy'] + packet.r = packet_params['r'] + + cmontecarlo_methods.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']) + + assert_almost_equal(model.js[packet.current_shell_id], expected_params['j']) + assert_almost_equal(model.nubars[packet.current_shell_id], expected_params['nubar']) + + @pytest.mark.parametrize( ['packet_params', 'j_blue_idx', 'expected'], [({'nu': 0.1, 'mu': 0.3, 'r': 7.5e14}, 0, 8.998643292289723), @@ -273,6 +352,56 @@ def test_increment_j_blue_estimator(packet_params, j_blue_idx, expected, packet, 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(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'] + + cmontecarlo_methods.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(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'] + + cmontecarlo_methods.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. @@ -298,50 +427,110 @@ def test_montecarlo_line_scatter(packet_params, expected_params, packet, model, assert_almost_equal(packet.next_line_id, expected_params['next_line_id']) -# TODO: redesign subsequent tests according to tests written above. -@pytest.mark.skipif(True, reason='Bad test design') -def test_move_packet(): - doppler_factor = 0.9998556693818854 - cmontecarlo_tests.test_move_packet.restype = c_double - assert_almost_equal(cmontecarlo_tests.test_move_packet(), - doppler_factor) + +""" +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 -def test_move_packet_across_shell_boundary(): - assert cmontecarlo_tests.test_move_packet_across_shell_boundary() +@pytest.mark.skipif(True, reason="Yet to be written.") +def test_montecarlo_main_loop(packet, model, mt_state): + pass -def test_montecarlo_one_packet(): - assert cmontecarlo_tests.test_montecarlo_one_packet() +@pytest.mark.skipif(True, reason="Yet to be written.") +def test_montecarlo_event_handler(packet, model, mt_state): + pass -def test_montecarlo_one_packet_loop(): - assert cmontecarlo_tests.test_montecarlo_one_packet_loop() == 0 +""" +Not Yet Relevant Tests: +----------------------- +The tests written further (till next block comment is encountered) are for the +methods related to Continuum interactions. These are not required to be tested +on current master and can be skipped for now. +""" + + +@pytest.mark.skipif(True, reason="Not yet relevant") +@pytest.mark.parametrize( + ['packet_params', 'expected'], + [({'nu': 0.1, 'mu': 0.3, 'r': 7.5e14}, 2.5010827921809502e+26), + ({'nu': 0.2, 'mu': -.3, 'r': 7.7e14}, 3.123611229395459e+25)] +) +def test_bf_cross_section(packet_params, expected, packet, model): + packet.nu = packet_params['nu'] + packet.mu = packet_params['mu'] + packet.r = packet_params['r'] + + cmontecarlo_methods.rpacket_doppler_factor.restype = c_double + doppler_factor = cmontecarlo_methods.rpacket_doppler_factor(byref(packet), byref(model)) + comov_nu = packet.nu * doppler_factor + + cmontecarlo_methods.bf_cross_section.restype = c_double + obtained = cmontecarlo_methods.bf_cross_section(byref(model), c_int64(0), + c_double(comov_nu)) + + assert_almost_equal(obtained, expected) + + +# TODO: fix underlying method and update expected values in testcases. +# For loop is not being executed in original method, and hence bf_helper +# always remains zero. Reason for for loop not executed: +# "current_continuum_id = no_of_continuum edges" +@pytest.mark.skipif(True, reason="Not yet relevant") +@pytest.mark.parametrize( + ['packet_params', 'expected'], + [({'nu': 0.1, 'mu': 0.3, 'r': 7.5e14}, 0.0), + ({'nu': 0.2, 'mu': -.3, 'r': 7.7e14}, 0.0)] +) +def test_calculate_chi_bf(packet_params, expected, packet, model): + packet.nu = packet_params['nu'] + packet.mu = packet_params['mu'] + packet.r = packet_params['r'] + + cmontecarlo_methods.calculate_chi_bf(byref(packet), byref(model)) + assert_almost_equal(packet.chi_bf, expected) -def test_montecarlo_thomson_scatter(): - assert cmontecarlo_tests.test_montecarlo_thomson_scatter() +@pytest.mark.skipif(True, reason="Not yet relevant") +def test_montecarlo_continuum_event_handler(packet_params, continuum_status, expected, + packet, model, mt_state): + packet.chi_cont = packet_params['chi_cont'] + packet.chi_th = packet_params['chi_th'] + packet.chi_bf = packet.chi_cont - packet.chi_th + model.cont_status = continuum_status -def test_calculate_chi_bf(): - chi_bf = 1.0006697327643788 - cmontecarlo_tests.test_calculate_chi_bf.restype = c_double - assert_almost_equal(cmontecarlo_tests.test_calculate_chi_bf(), - chi_bf) + obtained = cmontecarlo_methods.montecarlo_continuum_event_handler(byref(packet), + byref(model), byref(mt_state)) -@pytest.mark.xfail -def test_montecarlo_bound_free_scatter(): - assert cmontecarlo_tests.test_montecarlo_bound_free_scatter() == 1 +@pytest.mark.skipif(True, reason="Not yet relevant") +def test_montecarlo_free_free_scatter(packet, model, mt_state): + cmontecarlo_methods.montecarlo_free_free_scatter(byref(packet), byref(model), + c_double(1.e13), byref(mt_state)) + assert_equal(packet.status, TARDIS_PACKET_STATUS_REABSORBED) -@pytest.mark.xfail -def test_bf_cross_section(): - bf_cross_section = 0.0 - cmontecarlo_tests.test_bf_cross_section.restype = c_double - assert_almost_equal(cmontecarlo_tests.test_bf_cross_section(), - bf_cross_section) +@pytest.mark.skipif(True, reason="Not yet relevant") +def test_montecarlo_bound_free_scatter(packet, model, mt_state): + cmontecarlo_methods.montecarlo_bound_free_scatter(byref(packet), byref(model), + c_double(1.e13), byref(mt_state)) -def test_montecarlo_free_free_scatter(): - assert cmontecarlo_tests.test_montecarlo_free_free_scatter() == 2 + assert_equal(packet.status, TARDIS_PACKET_STATUS_REABSORBED) diff --git a/tardis/montecarlo/tests/test_rpacket.py b/tardis/montecarlo/tests/test_rpacket.py deleted file mode 100644 index b167c42757a..00000000000 --- a/tardis/montecarlo/tests/test_rpacket.py +++ /dev/null @@ -1,122 +0,0 @@ -import os -import random -from ctypes import CDLL, c_double - -import pytest -import numpy as np - -from tardis import __path__ as path - -test_path = os.path.join(path[0], 'montecarlo', 'test_montecarlo.so') - -#pytestmark = pytest.mark.skipif(True, reason='problem with the files') - -tests = CDLL(test_path) - -np.random.seed(1) - - -def get_doubles(size=10): - from sys import float_info as floats - MAX_FLOAT, MIN_FLOAT = floats[0], floats[3] - return map(c_double, np.random.uniform( - MIN_FLOAT, MAX_FLOAT, size=10)) - -@pytest.fixture(params=get_doubles()) -def double_value(request): - return request.param - - -def get_integers(size=10): - MAX_INT, MIN_INT = 10000, -10000 - return np.random.randint( - MIN_INT, MAX_INT, size=10) - -@pytest.fixture(params=get_integers()) -def int_value(request): - return request.param - - -def get_unsigned_integers(size=10): - MAX_INT = 10000 - return np.random.randint( - 0, MAX_INT, size=10) - -@pytest.fixture(params=get_unsigned_integers()) -def unsigned_int_value(request): - return request.param - - -# Testing functions with float(C double) valued parameters - -def test_rpacket_get_nu(double_value): - assert tests.test_rpacket_get_nu(double_value) - - -def test_rpacket_get_mu(double_value): - assert tests.test_rpacket_get_mu(double_value) - - -def test_rpacket_get_energy(double_value): - assert tests.test_rpacket_get_energy(double_value) - - -def test_rpacket_get_r(double_value): - assert tests.test_rpacket_get_r(double_value) - - -def test_rpacket_get_tau_event(double_value): - assert tests.test_rpacket_get_tau_event(double_value) - - -def test_rpacket_get_nu_line(double_value): - assert tests.test_rpacket_get_nu_line(double_value) - - -def test_rpacket_get_d_boundary(double_value): - assert tests.test_rpacket_get_d_boundary(double_value) - - -def test_rpacket_get_d_electron(double_value): - assert tests.test_rpacket_get_d_electron(double_value) - - -def test_rpacket_get_d_line(double_value): - assert tests.test_rpacket_get_d_line(double_value) - - -# Testing functions with Unsigned Integer valued parameters - -def test_rpacket_get_current_shell_id(unsigned_int_value): - assert tests.test_rpacket_get_current_shell_id(unsigned_int_value) - - -def test_rpacket_get_next_line_id(unsigned_int_value): - assert tests.test_rpacket_get_next_line_id(unsigned_int_value) - - -# Testing functions with Integer valued parameters - -def test_rpacket_get_virtual_packet_flag(int_value): - assert tests.test_rpacket_get_virtual_packet_flag(int_value) - - -def test_rpacket_get_virtual_packet(int_value): - assert tests.test_rpacket_get_virtual_packet(int_value) - - -def test_rpacket_get_next_shell_id(int_value): - assert tests.test_rpacket_get_next_shell_id(int_value) - - -# Testing functions without any parameters -def test_rpacket_get_last_line(): - assert tests.test_rpacket_get_last_line() - - -def test_rpacket_get_close_line(): - assert tests.test_rpacket_get_close_line() - - -def test_rpacket_get_status(): - assert tests.test_rpacket_get_status()