From feee98472414c1a33485d7a492b0dc702c10658a Mon Sep 17 00:00:00 2001 From: Arjun Savel <35353555+arjunsavel@users.noreply.github.com> Date: Thu, 2 Jul 2020 07:50:28 -0700 Subject: [PATCH] Debugging relativity on numba_montecarlo (#1136) * removing kwarg from MonteCarloConfiguration, as jitclass does not support them * need to pass enable_full_relativity to configuration_initialize, as removed kwarg in MonteCarloConfiguration * move_r_packet and trace_packet both require montecarlo_configuration arg * r_packet object used in distance calculation in place of r_packet.r value * the full_relativity bool was being passed as a different kwarg * added montecarlo_configuration to args of trace_vpacket_volley and trace_vpacket, so that relativity attribute continues to get passed down * making sure the vpacket funcs that now require montecarlo_configuration get called with this arg * Added random seed; does not work with numba * Changed to numpy random seed * Commented out MonteCarloException block; will test how this impacts output * Need the MonteCarloException block to be run, otherwise large discrepancies between numba branch and master branch * Adding custom seeds to vpackets; few line adjustments * changed from njit to jit to ease debugging * More descriptive montecarloexception, with optional printing beforehand * Changed njit to jit in single_loop to ease debugging * Added single-packet debugging capabilities to MonteCarloRunner. * Removing plotting funcs from montecarlo * Added basic logging capabilities to montecarlo_numba * Generalized logger to *args * Tweaking configuration of logging in module * Logger can now handle kwargs * Renamed logger * Make logging more in line with other TARDIS loggers. * Share global variable DEBUG_MODE across modules with montecarlo init. * No longer calling removed plot_single_packet; also referring to __init__ DEBUG_MODE * Now referring to montecarlo __init__ DEBUG_MODE global * Remove dangling else statement * no longer need plotting in montecarlo_numba base.py * Remove import of function that no longer exists. * No longer holding global variables in montecarlo init * Importing mc_logger for global variables now * Moving log_decorator from base.py * Adding logger to new montecarlo_logger file * Add new logger to, remove print statements from r_packet.py * Ensure that DEBUG_MODE is referenced within the montecarlo_logger file * Allow the function being decorated to take its kwargs. * Apply wraps correctly * Make else block logic more apparent * Get rid of extra blank line * Got rid of a few todos * Add note about why decorator behavior will not change according to config * removing reference to logger in base * Added crude buffer * Increase buffer; check DEBUG_MODE checked during call, not init * Set default buffer in montecarlo logger to stated default * Allow specification of buffer in montecarlo yml * Make sure logger_buffer is passed through to montecarlo logger * Remove todo, question from log_decorator * Adding better docstrings to log_decorator * Add newline at end of montecarlo_logger file * Make logger config happen outside of decorator * Moved print statement to exception; * Added logging to file * No longer profiling calculate_distance_line * Add functionality to profile incomplete packet runs * Rename previous debug mode, as it was not for single packets * Added single packet seed to yml * Propagate single_packet_seed throughout configurations * Reference single packet seed from montecarloconfiguration * Add bool type to single packet seed type * Catch exception better, catch extra random seed * Set whole index, including energies, not just random seed * Change where loop is broken. * Move close_line_threshold check past montecarloexception * Add single-packet debug documentation and script * Add jitclass arg for single packet seed in montecarloconfiguration * Now the default for single_packet_seed cannot be a bool * Added new debug page to developer part of documentation * Now the single_packet_seed has to be not 1, not not False * Pause using log_decorator, remove try/except for nopython * Vpacket jitclass now allows for an int64 index * Exceptions args must be compile-time constants in nopython mode * Allow 0 nu_diff to not throw exception. * v_packets should have the same random seed as their r_packet * Fix MonteCarloException throw * Fix estimator typos * Make sure return is not made before relativity block * Rename esitmators to keep them consistent * Also rename the estimators in montecarlo.pyx * Rename estimators in numba_interface * Clean up estimator typos; change back to njit * Include Doppler factor in energy calculation for full_relativity * Add back not statement * Added relativistic Doppler factor * Add in Doppler factor relativity for interaction * Add angle aberration * add angle aberration to interactions * Update numba version * Remove relativity branching * Logger no longer prints to console if printing to file * Initialize configuration with external module, not jitclass * Remove extra repr of logger.handlers * Ensure the angle aberration is called on vpacket, remove jitclass reference * Convert montecarloconfiguration refs from jitclass to module * Add more references to global module, complete relativistic branching * Pass r_packet.mu to angle aberration calculation * Clarify angle aberration calc, include doppler factor to distance calc * Remove reference to MonteCarloConfiguration * Trying to delete C tests * Remove C montecarlo * Remove references to C module * Alter formal integral tests for python version * Rewrite formal integral in python --- docs/development/debug_numba.rst | 14 + docs/development/index.rst | 1 + tardis/io/schemas/montecarlo.yml | 13 + tardis/montecarlo/__init__.py | 2 +- tardis/montecarlo/base.py | 31 +- tardis/montecarlo/formal_integral.py | 322 +++- tardis/montecarlo/montecarlo.pyx | 378 ----- tardis/montecarlo/montecarlo_configuration.py | 5 + tardis/montecarlo/montecarlo_numba/base.py | 39 +- .../montecarlo_numba/interaction.py | 29 +- .../montecarlo_numba/montecarlo_logger.py | 67 + .../montecarlo_numba/numba_interface.py | 45 +- .../montecarlo/montecarlo_numba/r_packet.py | 223 ++- .../montecarlo_numba/single_packet_loop.py | 74 +- tardis/montecarlo/montecarlo_numba/vpacket.py | 33 +- tardis/montecarlo/setup_package.py | 44 +- tardis/montecarlo/src/abbrev.h | 24 - tardis/montecarlo/src/cmontecarlo.c | 1291 ----------------- tardis/montecarlo/src/cmontecarlo.h | 159 -- tardis/montecarlo/src/integrator.c | 343 ----- tardis/montecarlo/src/integrator.h | 19 - tardis/montecarlo/src/io.h | 32 - tardis/montecarlo/src/omp_helper.h | 7 - tardis/montecarlo/src/randomkit/LICENSE | 20 - tardis/montecarlo/src/randomkit/randomkit.h | 42 - tardis/montecarlo/src/randomkit/rk_isaac.c | 258 ---- tardis/montecarlo/src/randomkit/rk_isaac.h | 142 -- tardis/montecarlo/src/randomkit/rk_mt.c | 342 ----- tardis/montecarlo/src/randomkit/rk_mt.h | 194 --- .../montecarlo/src/randomkit/rk_primitive.c | 520 ------- .../montecarlo/src/randomkit/rk_primitive.h | 54 - tardis/montecarlo/src/randomkit/rk_sobol.c | 991 ------------- tardis/montecarlo/src/randomkit/rk_sobol.h | 173 --- tardis/montecarlo/src/rpacket.c | 55 - tardis/montecarlo/src/rpacket.h | 330 ----- tardis/montecarlo/src/status.h | 38 - tardis/montecarlo/src/storage.h | 102 -- tardis/montecarlo/tests/conftest.py | 7 - .../montecarlo/tests/test_formal_integral.py | 73 +- tardis/montecarlo/tests/test_montecarlo.py | 0 tardis/scripts/debug/run_numba_single.py | 22 + tardis/scripts/debug/run_numba_single.run.xml | 25 + .../scripts/debug/tardis_example_single.yml | 45 + tardis/simulation/base.py | 4 +- tardis/tests/test_tardis_full.py | 2 +- .../tests/test_tardis_full_formal_integral.py | 2 +- tardis_env3.yml | 2 +- 47 files changed, 889 insertions(+), 5749 deletions(-) create mode 100644 docs/development/debug_numba.rst delete mode 100644 tardis/montecarlo/montecarlo.pyx create mode 100644 tardis/montecarlo/montecarlo_configuration.py create mode 100644 tardis/montecarlo/montecarlo_numba/montecarlo_logger.py delete mode 100644 tardis/montecarlo/src/abbrev.h delete mode 100644 tardis/montecarlo/src/cmontecarlo.c delete mode 100644 tardis/montecarlo/src/cmontecarlo.h delete mode 100644 tardis/montecarlo/src/integrator.c delete mode 100644 tardis/montecarlo/src/integrator.h delete mode 100644 tardis/montecarlo/src/io.h delete mode 100644 tardis/montecarlo/src/omp_helper.h delete mode 100644 tardis/montecarlo/src/randomkit/LICENSE delete mode 100644 tardis/montecarlo/src/randomkit/randomkit.h delete mode 100644 tardis/montecarlo/src/randomkit/rk_isaac.c delete mode 100644 tardis/montecarlo/src/randomkit/rk_isaac.h delete mode 100644 tardis/montecarlo/src/randomkit/rk_mt.c delete mode 100644 tardis/montecarlo/src/randomkit/rk_mt.h delete mode 100644 tardis/montecarlo/src/randomkit/rk_primitive.c delete mode 100644 tardis/montecarlo/src/randomkit/rk_primitive.h delete mode 100644 tardis/montecarlo/src/randomkit/rk_sobol.c delete mode 100644 tardis/montecarlo/src/randomkit/rk_sobol.h delete mode 100644 tardis/montecarlo/src/rpacket.c delete mode 100644 tardis/montecarlo/src/rpacket.h delete mode 100644 tardis/montecarlo/src/status.h delete mode 100644 tardis/montecarlo/src/storage.h create mode 100644 tardis/montecarlo/tests/test_montecarlo.py create mode 100644 tardis/scripts/debug/run_numba_single.py create mode 100644 tardis/scripts/debug/run_numba_single.run.xml create mode 100644 tardis/scripts/debug/tardis_example_single.yml diff --git a/docs/development/debug_numba.rst b/docs/development/debug_numba.rst new file mode 100644 index 00000000000..5cd8c31d9e5 --- /dev/null +++ b/docs/development/debug_numba.rst @@ -0,0 +1,14 @@ +************************** +Debugging numba_montecarlo +************************** +To facilitate more in-depth debugging when interfacing with the `montecarlo_numba` +module, we provide a set of debugging configurations. PyCharm debugging +configurations, in addition to related scripts and .yml files, are contained in +`tardis.scripts.debug`. Currently, these include the ability to run TARDIS +in asingle-packet mode, with the packet seed identified at debug time. +`tardis_example_single.yml` is the configuration filethat is used to set up the +single-packet TARDIS run; `run_numba_single.py` is thePython script that runs +this .yml file; `run_numba_single.xml` is the PyCharmdebug configuration file +that can be used in conjunction with the above files. + + diff --git a/docs/development/index.rst b/docs/development/index.rst index e4fe42e68e9..03c4a9b2e2a 100644 --- a/docs/development/index.rst +++ b/docs/development/index.rst @@ -13,6 +13,7 @@ to the Astropy team for designing it. running_tests issues + debug_numba diff --git a/tardis/io/schemas/montecarlo.yml b/tardis/io/schemas/montecarlo.yml index b0638270e6c..2c32e9598cb 100644 --- a/tardis/io/schemas/montecarlo.yml +++ b/tardis/io/schemas/montecarlo.yml @@ -68,6 +68,19 @@ properties: default: false description: Enables a more complete treatment of relativitic effects. This includes angle aberration as well as use of the fully general Doppler formula. + debug_packets: + type: boolean + default: false + description: Decide whether to go into debugging mode. + logger_buffer: + type: number + default: 1 + description: Provides option to not log every line. + single_packet_seed: + type: + - number + default: -1 + description: If debug_packets is true, this is the seed for the only packet. required: - no_of_packets diff --git a/tardis/montecarlo/__init__.py b/tardis/montecarlo/__init__.py index 222ee95b0b8..d8d3283c243 100644 --- a/tardis/montecarlo/__init__.py +++ b/tardis/montecarlo/__init__.py @@ -1 +1 @@ -from tardis.montecarlo.base import * +from tardis.montecarlo.base import * \ No newline at end of file diff --git a/tardis/montecarlo/base.py b/tardis/montecarlo/base.py index 5665778ae05..507f61aeebb 100644 --- a/tardis/montecarlo/base.py +++ b/tardis/montecarlo/base.py @@ -13,10 +13,13 @@ from tardis.util.base import quantity_linspace from tardis.io.util import HDFWriterMixin -from tardis.montecarlo import montecarlo, packet_source as source +from tardis.montecarlo import packet_source as source from tardis.montecarlo.formal_integral import FormalIntegrator +from tardis.montecarlo import montecarlo_configuration as mc_config_module + from tardis.montecarlo.montecarlo_numba import montecarlo_radial1d +from tardis.montecarlo.montecarlo_numba import montecarlo_logger as mc_logger from tardis.montecarlo.montecarlo_numba.numba_interface import ( configuration_initialize) @@ -62,7 +65,8 @@ def __init__(self, seed, spectrum_frequency, virtual_spectrum_range, enable_full_relativity, inner_boundary_albedo, line_interaction_type, integrator_settings, v_packet_settings, spectrum_method, - packet_source=None): + packet_source=None, debug_packets=False, + logger_buffer=1, single_packet_seed=None): self.seed = seed if packet_source is None: @@ -77,11 +81,17 @@ def __init__(self, seed, spectrum_frequency, virtual_spectrum_range, self.inner_boundary_albedo = inner_boundary_albedo self.enable_full_relativity = enable_full_relativity self.line_interaction_type = line_interaction_type + self.single_packet_seed = single_packet_seed self.integrator_settings = integrator_settings self.v_packet_settings = v_packet_settings self.spectrum_method = spectrum_method self._integrator = None self._spectrum_integrated = None + + # set up logger based on config + mc_logger.DEBUG_MODE = debug_packets + mc_logger.BUFFER = logger_buffer + if self.spectrum_method == 'integrated': self.optional_hdf_properties.append('spectrum_integrated') @@ -100,8 +110,9 @@ def _initialize_estimator_arrays(self, tau_sobolev_shape): # Estimators self.j_estimator = np.zeros(tau_sobolev_shape[1], dtype=np.float64) self.nu_bar_estimator = np.zeros(tau_sobolev_shape[1], dtype=np.float64) - self.j_b_lu_estimator = np.zeros(tau_sobolev_shape) - self.edot_lu_estimator = np.zeros(tau_sobolev_shape) + self.j_blue_estimator = np.zeros(tau_sobolev_shape) + self.Edotlu_estimator = np.zeros(tau_sobolev_shape) + # TODO: this is the wrong attribute naming style. def _initialize_geometry_arrays(self, model): """ @@ -227,10 +238,9 @@ def run(self, model, plasma, no_of_packets, self._initialize_packets(model.t_inner.value, no_of_packets) - montecarlo_configuration = configuration_initialize( - self, no_of_virtual_packets) - - montecarlo_radial1d(model, plasma, self, montecarlo_configuration) + montecarlo_configuration = configuration_initialize(self, + no_of_virtual_packets) + montecarlo_radial1d(model, plasma, self) #montecarlo.montecarlo_radial1d( # model, plasma, self, # virtual_packet_flag=no_of_virtual_packets, @@ -450,4 +460,7 @@ def from_config(cls, config, packet_source=None): integrator_settings=config.spectrum.integrated, v_packet_settings=config.spectrum.virtual, spectrum_method=config.spectrum.method, - packet_source=packet_source) + packet_source=packet_source, + debug_packets=config.montecarlo.debug_packets, + logger_buffer=config.montecarlo.logger_buffer, + single_packet_seed=config.montecarlo.single_packet_seed) diff --git a/tardis/montecarlo/formal_integral.py b/tardis/montecarlo/formal_integral.py index 2d973c44603..872e0410e2c 100644 --- a/tardis/montecarlo/formal_integral.py +++ b/tardis/montecarlo/formal_integral.py @@ -5,14 +5,30 @@ from scipy.interpolate import interp1d from astropy import units as u from tardis import constants as const +from numba import jitclass, njit -from tardis.montecarlo.montecarlo import formal_integral + +from tardis.montecarlo.montecarlo_numba import njit_dict + +# from tardis.montecarlo.montecarlo import formal_integral from tardis.montecarlo.spectrum import TARDISSpectrum +C_INV = 3.33564e-11 +M_PI = np.arccos(-1) +KB_CGS = 1.3806488e-16 +H_CGS = 6.62606957e-27 class IntegrationError(Exception): pass +# integrator_spec = [ +# ('model', float64), +# ('plasma', float64), +# ('runner', float64), +# ('points', int64) +# ] +# +# @jitclass(integrator_spec) class FormalIntegrator(object): def __init__(self, model, plasma, runner, points=1000): @@ -122,7 +138,7 @@ def make_source_function(self): Edotlu_norm_factor = (1 / (runner.time_of_simulation * model.volume)) exptau = 1 - np.exp(- plasma.tau_sobolevs) - Edotlu = Edotlu_norm_factor * exptau * runner.edot_lu_estimator + Edotlu = Edotlu_norm_factor * exptau * runner.Edotlu_estimator # The following may be achieved by calling the appropriate plasma # functions @@ -131,7 +147,7 @@ def make_source_function(self): model.volume)).to("1/(cm^2 s)").value # Jbluelu should already by in the correct order, i.e. by wavelength of # the transition l->u - Jbluelu = runner.j_b_lu_estimator * Jbluelu_norm_factor + Jbluelu = runner.j_blue_estimator * Jbluelu_norm_factor upper_level_index = atomic_data.lines.index.droplevel('level_number_lower') e_dot_lu = pd.DataFrame(Edotlu, index=upper_level_index) @@ -219,3 +235,303 @@ def interpolate_integrator_quantities(self, att_S_ul, Jredlu, Jredlu = Jredlu.clip(0.) e_dot_u = e_dot_u.clip(0.) return att_S_ul, Jredlu, Jbluelu, e_dot_u + + def formal_integral(self, nu, N): + # TODO: get rid of storage later on + + + 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') + L = self._formal_integral( + self.model.t_inner.value, + nu, + nu.shape[0], + att_S_ul, + Jred_lu, + Jblue_lu, + N + ) + return np.array(L, np.NPY_DOUBLE, nu.shape[0]) + + + + def _formal_integral(self, iT, inu, inu_size, att_S_ul, Jred_lu, Jblue_lu, N): + # todo: add all the original todos + # 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_tau = size_line * size_shell + finished_nus = 0 + + R_ph = self.runner.r_inner_i[0] + R_max = self.runner.r_outer_i[size_shell - 1] + pp = np.zeros(N) # check + exp_tau = np.zeros(size_tau) + # TODO: multiprocessing + offset = 0 + i = 0 + size_z = 0 + idx_nu_start = 0 + direction = 0 + first = 0 + I_nu = np.zeros(N) + shell_id = np.zeros(2 * self.runner.no_of_shells_i) # check + # instantiate more variables here, maybe? + + # 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) + + # done with instantiation + # now loop over wavelength in spectrum + for nu_idx in range(inu_size): + nu = inu[nu_idx] + # now loop over discrete values along line + for p_idx in range(1, N): + escat_contrib = 0 + p = pp[p_idx] + + # initialize z intersections for p values + size_z = self.populate_z(p, z, shell_id) # check returns + + # initialize I_nu + if p <= R_ph: + I_nu[p_idx] = intensity_black_body(nu * z[0], iT) + else: + I_nu[p_idx] = 0 + + # find first contributing lines + nu_start = nu * z[0] + nu_end = nu * z[1] + self.line_search(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; + + # 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 + nu_end = nu * z[i + 1] + while pline < self.model.ine_list_nu + size_line: # check ;pline + # increment all pointers simulatenously + pline += 1 + pexp_tau += 1 + patt_S_ul += 1 + pJblue_lu += 1 + + if (pline[0] < nu_end): + break + + # calculate e-scattering optical depth to next resonance point + zend = self.model.time_explosion / C_INV * (1. - pline / nu) # check + + if first == 1: + # first contribution to integration + # NOTE: this treatment of I_nu_b (given + # 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]); + first = 0; + else: + # Account for e-scattering, c.f. Eqs 27, 28 in Lucy 1999 + Jkkp = 0.5 * (pJred_lu + pJblue_lu); + 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; + # // Lucy 1999, Eq 26 + I_nu[p_idx] = I_nu[p_idx] * (pexp_tau) + patt_S_ul # 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) + zend = self.model.time_explosion / C_INV * (1. - nu_end / nu) # check + escat_contrib += (zend - zstart) * escat_op * ( + Jkkp - I_nu[p_idx]) + zstart = zend + + if i < size_z - 1: + # advance pointers + direction = shell_id[i+1] - shell_id[i] + pexp_tau += direction * size_line + patt_S_ul += direction * size_line + pJred_lu += direction * size_line + pJblue_lu += direction * size_line + I_nu[p_idx] *= p; + L[nu_idx] = 8 * M_PI * M_PI * trapezoid_integration(I_nu, R_max / N, + N) + # something pragma op atomic + return L + + + # @njit(**njit_dict) + def populate_z(self, p, oz, oshell_id): + """Calculate p line intersections + + This function calculates the intersection points of the p-line with + each shell + + Inputs: + :p: (double) distance of the integration line to the center + :oz: (array of doubles) will be set with z values. the array is truncated + by the value `1`. + :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 + print(N) + inv_t = 1/self.model.time_explosion + z = 0 + offset = N + + if p <= self.model.r_inner_i[0]: + # intersect the photosphere + for i in range(N): + oz[i] = 1 - self.calculate_z(r[i], p, inv_t) + oshell_id[i] = i + return N + else: + # no intersection with photosphere + # that means we intersect each shell twice + for i in range(N): + z = self.calculate_z(r[i], p, inv_t) + if z == 0: + continue + if offset == N: + offset = i + # calculate the index in the resulting array + i_low = N - i - 1 # the far intersection with the shell + i_up = N + i - 2 * offset # the nearer intersection with the shell + + # setting the arrays; check return them? + oz[i_low] = 1 + z + oshell_id[i_low] = i + oz[i_up] = 1 - z + oshell_id[i_up] = i + return 2 * (N - offset) + + # @njit(**njit_dict) + def calculate_z(self, r, p, inv_t): + """ Calculate distance to p line + + Calculate half of the length of the p-line inside a shell + of radius r in terms of unit length (c * t_exp). + If shell and p-line do not intersect, return 0. + + Inputs: + :r: (double) radius of the shell + :p: (double) distance of the p-line to the center of the supernova + :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 + 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 + + 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: + ret_val = self.reverse_binary_search(nu, nu_insert, imin, imax, result) + result = result + 1 + return ret_val + + @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. + + 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): + # TODO: replace with np.trapz? + result = (array[0] + array[N - 1]) / 2; + for idx in range(1, N-1): + result += array[idx] + return result * h + +@njit +def intensity_black_body(nu, T): + if nu == 0: + return np.nan # to avoid ZeroDivisionError + beta_rad = 1 / (KB_CGS * T) + coefficient = 2 * H_CGS * C_INV * C_INV + return coefficient * nu * nu * nu / (np.exp(H_CGS * nu * beta_rad) - 1) + +@njit(**njit_dict) +def calculate_p_values(R_max, N, opp): + for i in range(N): + opp[i] = R_max / (N - 1) * (i) + return opp \ No newline at end of file diff --git a/tardis/montecarlo/montecarlo.pyx b/tardis/montecarlo/montecarlo.pyx deleted file mode 100644 index d2d2def18bb..00000000000 --- a/tardis/montecarlo/montecarlo.pyx +++ /dev/null @@ -1,378 +0,0 @@ -# cython: profile=False -# cython: boundscheck=False -# cython: wraparound=False -# cython: cdivision=True -# cython: cdivision=True -# cython: language_level=3 - - - -import numpy as np -cimport numpy as np -from numpy cimport PyArray_DATA -from tardis import constants -from astropy import units -from libc.stdlib cimport free - -np.import_array() - - - -ctypedef np.int64_t int_type_t - -cdef extern from "numpy/arrayobject.h": - void PyArray_ENABLEFLAGS(np.ndarray arr, int flags) - -cdef c_array_to_numpy(void *ptr, int dtype, np.npy_intp N): - cdef np.ndarray arr = np.PyArray_SimpleNewFromData(1, &N, dtype, ptr) - PyArray_ENABLEFLAGS(arr, np.NPY_OWNDATA) - return arr - -cdef extern from "src/cmontecarlo.h": - ctypedef enum ContinuumProcessesStatus: - CONTINUUM_OFF = 0 - CONTINUUM_ON = 1 - - cdef int LOG_VPACKETS - - ctypedef struct photo_xsect_1level: - double *nu - double *x_sect - int_type_t no_of_points - - ctypedef struct storage_model_t: - double *packet_nus - double *packet_mus - double *packet_energies - double *output_nus - double *output_energies - double *last_interaction_in_nu - int_type_t *last_line_interaction_in_id - int_type_t *last_line_interaction_out_id - int_type_t *last_line_interaction_shell_id - int_type_t *last_interaction_type - int_type_t *last_interaction_out_type - int_type_t no_of_packets - int_type_t no_of_shells - int_type_t no_of_shells_i - double *r_inner - double *r_outer - double *r_inner_i - double *r_outer_i - double *v_inner - double time_explosion - double inverse_time_explosion - double *electron_densities - double *electron_densities_i - double *inverse_electron_densities - double *line_list_nu - double *line_lists_tau_sobolevs - double *line_lists_tau_sobolevs_i - double *continuum_list_nu - int_type_t line_lists_tau_sobolevs_nd - double *line_lists_j_blues - double *line_lists_Edotlu - int_type_t line_lists_j_blues_nd - int_type_t no_of_lines - int_type_t no_of_edges - int_type_t line_interaction_id - double *transition_probabilities - int_type_t transition_probabilities_nd - int_type_t *line2macro_level_upper - int_type_t *macro_block_references - int_type_t *transition_type - int_type_t *destination_level_id - int_type_t *transition_line_id - double *js - double *nubars - double spectrum_virt_start_nu - double spectrum_virt_end_nu - double spectrum_start_nu - double spectrum_delta_nu - double spectrum_end_nu - double *spectrum_virt_nu - double sigma_thomson - double inverse_sigma_thomson - double inner_boundary_albedo - int_type_t reflective_inner_boundary - photo_xsect_1level ** photo_xsect - double *chi_ff_factor - double *t_electrons - double *l_pop - double *l_pop_r - ContinuumProcessesStatus cont_status - double *virt_packet_nus - double *virt_packet_energies - double *virt_packet_last_interaction_in_nu - int_type_t *virt_packet_last_interaction_type - int_type_t *virt_packet_last_line_interaction_in_id - int_type_t *virt_packet_last_line_interaction_out_id - int_type_t virt_packet_count - int_type_t virt_array_size - int_type_t kpacket2macro_level - int_type_t *cont_edge2macro_level - double *photo_ion_estimator - double *stim_recomb_estimator - int_type_t *photo_ion_estimator_statistics - double *bf_heating_estimator - double *ff_heating_estimator - double *stim_recomb_cooling_estimator - int full_relativity - double survival_probability - double tau_russian - double *tau_bias - int enable_biasing - - void montecarlo_main_loop(storage_model_t * storage, int_type_t virtual_packet_flag, int nthreads, unsigned long seed) - -cdef extern from "src/integrator.h": - double *_formal_integral( - const storage_model_t *storage, - double T, - double *nu, - int_type_t nu_size, - double *att_S_ul, - double *Jred_lu, - double *Jblue_lu, - int N) - - - -cdef initialize_storage_model(model, plasma, runner, storage_model_t *storage): - """ - - Parameters - ---------- - model - plasma - runner: tardis.montecarlo.base.MontecarloRunner - storage - - Returns - ------- - - """ - storage.no_of_packets = runner.input_nu.size - storage.packet_nus = PyArray_DATA(runner.input_nu) - storage.packet_mus = PyArray_DATA(runner.input_mu) - storage.packet_energies = PyArray_DATA(runner.input_energy) - - # Setup of structure - storage.no_of_shells = model.no_of_shells - - - storage.r_inner = PyArray_DATA(runner.r_inner_cgs) - storage.r_outer = PyArray_DATA(runner.r_outer_cgs) - storage.v_inner = PyArray_DATA(runner.v_inner_cgs) - - # Setup the rest - # times - storage.time_explosion = model.time_explosion.to('s').value - storage.inverse_time_explosion = 1.0 / storage.time_explosion - #electron density - storage.electron_densities = PyArray_DATA( - plasma.electron_densities.values) - - runner.inverse_electron_densities = ( - 1.0 / plasma.electron_densities.values) - storage.inverse_electron_densities = PyArray_DATA( - runner.inverse_electron_densities) - # Switch for continuum processes - storage.cont_status = CONTINUUM_OFF - # Continuum data - cdef np.ndarray[double, ndim=1] continuum_list_nu - cdef np.ndarray[double, ndim=1] l_pop - cdef np.ndarray[double, ndim=1] l_pop_r - - if storage.cont_status == CONTINUUM_ON: - continuum_list_nu = np.array([9.0e14, 8.223e14, 6.0e14, 3.5e14, 3.0e14]) # sorted list of threshold frequencies - storage.continuum_list_nu = continuum_list_nu.data - storage.no_of_edges = continuum_list_nu.size - l_pop = np.ones(storage.no_of_shells * continuum_list_nu.size, dtype=np.float64) - storage.l_pop = l_pop.data - l_pop_r = np.ones(storage.no_of_shells * continuum_list_nu.size, dtype=np.float64) - storage.l_pop_r = l_pop_r.data - - # Line lists - storage.no_of_lines = plasma.atomic_data.lines.nu.values.size - storage.line_list_nu = PyArray_DATA(plasma.atomic_data.lines.nu.values) - runner.line_lists_tau_sobolevs = ( - plasma.tau_sobolevs.values.flatten(order='F') - ) - storage.line_lists_tau_sobolevs = PyArray_DATA( - runner.line_lists_tau_sobolevs - ) - storage.line_lists_j_blues = PyArray_DATA( - runner.j_b_lu_estimator) - - storage.line_lists_Edotlu = PyArray_DATA( - runner.edot_lu_estimator) - - storage.line_interaction_id = runner.get_line_interaction_id( - runner.line_interaction_type) - - # macro atom & downbranch - if storage.line_interaction_id >= 1: - runner.transition_probabilities = ( - plasma.transition_probabilities.values.flatten(order='F') - ) - storage.transition_probabilities = PyArray_DATA( - runner.transition_probabilities - ) - storage.transition_probabilities_nd = ( - plasma.transition_probabilities.values.shape[0]) - storage.line2macro_level_upper = PyArray_DATA( - plasma.atomic_data.lines_upper2macro_reference_idx) - storage.macro_block_references = PyArray_DATA( - plasma.atomic_data.macro_atom_references['block_references'].values) - storage.transition_type = PyArray_DATA( - plasma.atomic_data.macro_atom_data['transition_type'].values) - - # Destination level is not needed and/or generated for downbranch - storage.destination_level_id = PyArray_DATA( - plasma.atomic_data.macro_atom_data['destination_level_idx'].values) - storage.transition_line_id = PyArray_DATA( - plasma.atomic_data.macro_atom_data['lines_idx'].values) - - storage.output_nus = PyArray_DATA(runner._output_nu) - storage.output_energies = PyArray_DATA(runner._output_energy) - - storage.last_line_interaction_in_id = PyArray_DATA( - runner.last_line_interaction_in_id) - storage.last_line_interaction_out_id = PyArray_DATA( - runner.last_line_interaction_out_id) - storage.last_line_interaction_shell_id = PyArray_DATA( - runner.last_line_interaction_shell_id) - storage.last_interaction_type = PyArray_DATA( - runner.last_interaction_type) - storage.last_interaction_in_nu = PyArray_DATA( - runner.last_interaction_in_nu) - - storage.js = PyArray_DATA(runner.j_estimator) - storage.nubars = PyArray_DATA(runner.nu_bar_estimator) - - storage.spectrum_start_nu = runner.spectrum_frequency.to('Hz').value.min() - storage.spectrum_end_nu = runner.spectrum_frequency.to('Hz').value.max() - # TODO: Linspace handling for virtual_spectrum_range - storage.spectrum_virt_start_nu = runner.virtual_spectrum_range.stop.to('Hz', units.spectral()).value - storage.spectrum_virt_end_nu = runner.virtual_spectrum_range.start.to('Hz', units.spectral()).value - storage.spectrum_delta_nu = runner.spectrum_frequency.to('Hz').value[1] - runner.spectrum_frequency.to('Hz').value[0] - - storage.spectrum_virt_nu = PyArray_DATA( - runner._montecarlo_virtual_luminosity.value) - - storage.sigma_thomson = runner.sigma_thomson.cgs.value - storage.inverse_sigma_thomson = 1.0 / storage.sigma_thomson - storage.reflective_inner_boundary = runner.enable_reflective_inner_boundary - storage.inner_boundary_albedo = runner.inner_boundary_albedo - storage.full_relativity = runner.enable_full_relativity - - storage.tau_russian = runner.v_packet_settings['tau_russian'] - storage.survival_probability = runner.v_packet_settings['survival_probability'] - storage.enable_biasing = runner.v_packet_settings['enable_biasing'] - - if runner.v_packet_settings['enable_biasing']: - # Calculate the integrated electron scattering optical depth - # at all cell interfaces. - runner.tau_bias = np.zeros(len(runner.r_inner_cgs) + 1) - runner.tau_bias[:-1] = ( - ((runner.r_outer_cgs - runner.r_inner_cgs) * - plasma.electron_densities.values * - runner.sigma_thomson.cgs.value)[::-1].cumsum()[::-1] - ) - storage.tau_bias = PyArray_DATA(runner.tau_bias) - - # Data for continuum implementation - cdef np.ndarray[double, ndim=1] t_electrons = plasma.t_electrons - storage.t_electrons = t_electrons.data - -def montecarlo_radial1d(model, plasma, runner, int_type_t virtual_packet_flag=0, - int nthreads=4,last_run=False): - """ - Parameters - ---------- - model : `tardis.model_radial_oned.ModelRadial1D` - complete model - param photon_packets : PacketSource object - photon packets - - Returns - ------- - output_nus : `numpy.ndarray` - output_energies : `numpy.ndarray` - - TODO - np.ndarray[double, ndim=1] line_list_nu, - np.ndarray[double, ndim=2] tau_lines, - np.ndarray[double, ndim=1] ne, - double packet_energy, - np.ndarray[double, ndim=2] p_transition, - np.ndarray[int_type_t, ndim=1] type_transition, - np.ndarray[int_type_t, ndim=1] target_level_id, - np.ndarray[int_type_t, ndim=1] target_line_id, - np.ndarray[int_type_t, ndim=1] unroll_reference, - np.ndarray[int_type_t, ndim=1] line2level, - int_type_t log_packets, - int_type_t do_scatter - """ - - cdef storage_model_t storage - - initialize_storage_model(model, plasma, runner, &storage) - - montecarlo_main_loop(&storage, virtual_packet_flag, nthreads, runner.seed) - runner.virt_logging = LOG_VPACKETS - if LOG_VPACKETS != 0: - runner.virt_packet_nus = c_array_to_numpy(storage.virt_packet_nus, np.NPY_DOUBLE, storage.virt_packet_count) - runner.virt_packet_energies = c_array_to_numpy(storage.virt_packet_energies, np.NPY_DOUBLE, storage.virt_packet_count) - runner.virt_packet_last_interaction_in_nu = c_array_to_numpy(storage.virt_packet_last_interaction_in_nu, np.NPY_DOUBLE, storage.virt_packet_count) - runner.virt_packet_last_interaction_type = c_array_to_numpy(storage.virt_packet_last_interaction_type, np.NPY_INT64, storage.virt_packet_count) - runner.virt_packet_last_line_interaction_in_id = c_array_to_numpy(storage.virt_packet_last_line_interaction_in_id, np.NPY_INT64, - storage.virt_packet_count) - runner.virt_packet_last_line_interaction_out_id = c_array_to_numpy(storage.virt_packet_last_line_interaction_out_id, np.NPY_INT64, - storage.virt_packet_count) - else: - runner.virt_packet_nus = np.zeros(0) - runner.virt_packet_energies = np.zeros(0) - runner.virt_packet_last_interaction_in_nu = np.zeros(0) - runner.virt_packet_last_interaction_type = np.zeros(0) - runner.virt_packet_last_line_interaction_in_id = np.zeros(0) - runner.virt_packet_last_line_interaction_out_id = np.zeros(0) - - -# This will be a method of the Simulation object -def formal_integral(self, nu, N): - cdef storage_model_t storage - - initialize_storage_model(self.model, self.plasma, self.runner, &storage) - - res = self.make_source_function() - - storage.no_of_shells_i = len(self.runner.r_inner_i) - storage.r_inner_i = PyArray_DATA(self.runner.r_inner_i) - storage.r_outer_i = PyArray_DATA(self.runner.r_outer_i) - - storage.electron_densities_i = PyArray_DATA( - self.runner.electron_densities_integ) - self.runner.line_lists_tau_sobolevs_i = ( - self.runner.tau_sobolevs_integ.flatten(order='F') - ) - storage.line_lists_tau_sobolevs_i = PyArray_DATA( - self.runner.line_lists_tau_sobolevs_i - ) - - att_S_ul = res[0].flatten(order='F') - Jred_lu = res[1].flatten(order='F') - Jblue_lu = res[2].flatten(order='F') - - cdef double *L = _formal_integral( - &storage, - self.model.t_inner.value, - PyArray_DATA(nu), - nu.shape[0], - PyArray_DATA(att_S_ul), - PyArray_DATA(Jred_lu), - PyArray_DATA(Jblue_lu), - N - ) - return c_array_to_numpy(L, np.NPY_DOUBLE, nu.shape[0]) diff --git a/tardis/montecarlo/montecarlo_configuration.py b/tardis/montecarlo/montecarlo_configuration.py new file mode 100644 index 00000000000..f81b6e9b186 --- /dev/null +++ b/tardis/montecarlo/montecarlo_configuration.py @@ -0,0 +1,5 @@ +full_relativity = True +single_packet_seed = -1 +temporary_v_packet_bins = None +number_of_vpackets = 0 +line_interaction_type = None \ No newline at end of file diff --git a/tardis/montecarlo/montecarlo_numba/base.py b/tardis/montecarlo/montecarlo_numba/base.py index a08d5eb735f..fe5ad56cb22 100644 --- a/tardis/montecarlo/montecarlo_numba/base.py +++ b/tardis/montecarlo/montecarlo_numba/base.py @@ -1,16 +1,21 @@ -from numba import prange, njit +from numba import prange, njit, jit +import logging import numpy as np -from tardis.montecarlo.montecarlo_numba.r_packet import RPacket, PacketStatus + +from tardis.montecarlo.montecarlo_numba.r_packet import ( + RPacket, PacketStatus, MonteCarloException) from tardis.montecarlo.montecarlo_numba.numba_interface import ( PacketCollection, VPacketCollection, NumbaModel, numba_plasma_initialize, - Estimators, MonteCarloConfiguration, configuration_initialize) + Estimators, configuration_initialize) + +from tardis.montecarlo import montecarlo_configuration as montecarlo_configuration from tardis.montecarlo.montecarlo_numba.single_packet_loop import ( single_packet_loop) from tardis.montecarlo.montecarlo_numba import njit_dict -def montecarlo_radial1d(model, plasma, runner, montecarlo_configuration): +def montecarlo_radial1d(model, plasma, runner): packet_collection = PacketCollection( runner.input_nu, runner.input_mu, runner.input_energy, runner._output_nu, runner._output_energy @@ -20,19 +25,18 @@ def montecarlo_radial1d(model, plasma, runner, montecarlo_configuration): model.time_explosion.to('s').value) numba_plasma = numba_plasma_initialize(plasma) estimators = Estimators(runner.j_estimator, runner.nu_bar_estimator, - runner.j_b_lu_estimator, runner.edot_lu_estimator) + runner.j_blue_estimator, runner.Edotlu_estimator) v_packets_energy_hist = montecarlo_main_loop( packet_collection, numba_model, numba_plasma, estimators, - runner.spectrum_frequency.value, montecarlo_configuration) + runner.spectrum_frequency.value) runner._montecarlo_virtual_luminosity.value[:] = v_packets_energy_hist @njit(**njit_dict, nogil=True) def montecarlo_main_loop(packet_collection, numba_model, numba_plasma, - estimators, spectrum_frequency, - montecarlo_configuration): + estimators, spectrum_frequency): """ This is the main loop of the MonteCarlo routine that generates packets and sends them through the ejecta. @@ -50,18 +54,23 @@ def montecarlo_main_loop(packet_collection, numba_model, numba_plasma, delta_nu = spectrum_frequency[1] - spectrum_frequency[0] for i in prange(len(output_nus)): - np.random.seed(r_packet.seed) + if montecarlo_configuration.single_packet_seed != -1: + i = montecarlo_configuration.single_packet_seed r_packet = RPacket(numba_model.r_inner[0], packet_collection.packets_input_mu[i], packet_collection.packets_input_nu[i], packet_collection.packets_input_energy[i], i) - + + # We want to set the seed correctly per user; otherwise, random. + np.random.seed(i) vpacket_collection = VPacketCollection( spectrum_frequency, montecarlo_configuration.number_of_vpackets, montecarlo_configuration.temporary_v_packet_bins) - single_packet_loop(r_packet, numba_model, numba_plasma, estimators, - vpacket_collection, montecarlo_configuration) + loop = single_packet_loop(r_packet, numba_model, numba_plasma, estimators, + vpacket_collection) + # if loop and 'stop' in loop: + # raise MonteCarloException output_nus[i] = r_packet.nu @@ -75,6 +84,9 @@ def montecarlo_main_loop(packet_collection, numba_model, numba_plasma, v_packets_idx = np.floor((vpackets_nu - spectrum_frequency[0]) / delta_nu).astype(np.int64) + # if we're only in a single-packet mode + if montecarlo_configuration.single_packet_seed != -1: + break for j, idx in enumerate(v_packets_idx): if ((vpackets_nu[j] < spectrum_frequency[0]) or (vpackets_nu[j] > spectrum_frequency[-1])): @@ -84,5 +96,4 @@ def montecarlo_main_loop(packet_collection, numba_model, numba_plasma, packet_collection.packets_output_energy[:] = output_energies[:] packet_collection.packets_output_nu[:] = output_nus[:] - return v_packets_energy_hist - + return v_packets_energy_hist \ No newline at end of file diff --git a/tardis/montecarlo/montecarlo_numba/interaction.py b/tardis/montecarlo/montecarlo_numba/interaction.py index fda2c077135..482ec7564a6 100644 --- a/tardis/montecarlo/montecarlo_numba/interaction.py +++ b/tardis/montecarlo/montecarlo_numba/interaction.py @@ -3,9 +3,10 @@ from tardis.montecarlo.montecarlo_numba.numba_interface import ( LineInteractionType) - +from tardis.montecarlo import montecarlo_configuration as montecarlo_configuration from tardis.montecarlo.montecarlo_numba.r_packet import ( - get_doppler_factor, get_random_mu) + get_doppler_factor, get_inverse_doppler_factor, get_random_mu, + angle_aberration_CMF_to_LF) from tardis.montecarlo.montecarlo_numba.macro_atom import macro_atom @@ -24,14 +25,24 @@ def general_scatter(r_packet, time_explosion): distance : [type] [description] """ - old_doppler_factor = get_doppler_factor(r_packet.r, r_packet.mu, time_explosion) + old_doppler_factor = get_doppler_factor( + r_packet.r, + r_packet.mu, + time_explosion) comov_energy = r_packet.energy * old_doppler_factor comov_nu = r_packet.nu * old_doppler_factor r_packet.mu = get_random_mu() - inverse_new_doppler_factor = 1. / get_doppler_factor( + inverse_new_doppler_factor = get_inverse_doppler_factor( r_packet.r, r_packet.mu, time_explosion) r_packet.energy = comov_energy * inverse_new_doppler_factor r_packet.nu = comov_nu * inverse_new_doppler_factor + if montecarlo_configuration.full_relativity: + r_packet.mu = angle_aberration_CMF_to_LF( + r_packet, + time_explosion, + r_packet.mu + ) + """ void @@ -68,7 +79,8 @@ def line_scatter(r_packet, time_explosion, line_interaction_type, numba_plasma): # update last_interaction if line_interaction_type == LineInteractionType.SCATTER: - line_emission(r_packet, r_packet.next_line_id, time_explosion, numba_plasma) + line_emission(r_packet, r_packet.next_line_id, + time_explosion, numba_plasma) else: # includes both macro atom and downbranch - encoded in the transition probabilities emission_line_id = macro_atom(r_packet, numba_plasma) line_emission(r_packet, emission_line_id, time_explosion, @@ -84,6 +96,13 @@ def line_emission(r_packet, emission_line_id, time_explosion, r_packet.nu = numba_plasma.line_list_nu[ emission_line_id] / doppler_factor r_packet.next_line_id = emission_line_id + 1 + if montecarlo_configuration.full_relativity: + r_packet.mu = angle_aberration_CMF_to_LF( + r_packet, + time_explosion, + r_packet.mu + ) + diff --git a/tardis/montecarlo/montecarlo_numba/montecarlo_logger.py b/tardis/montecarlo/montecarlo_numba/montecarlo_logger.py new file mode 100644 index 00000000000..d1b5722dbbd --- /dev/null +++ b/tardis/montecarlo/montecarlo_numba/montecarlo_logger.py @@ -0,0 +1,67 @@ +import logging +from functools import wraps + +DEBUG_MODE = False +LOG_FILE = 'montecarlo_log.log' +BUFFER = 1 +ticker = 1 + +logger = logging.getLogger(__name__) +logger.setLevel(logging.DEBUG) +logger.handlers = [] + +if LOG_FILE: + logger.propagate = False + console_handler = logging.FileHandler(LOG_FILE) +else: + console_handler = logging.StreamHandler() + +console_handler.setLevel(logging.DEBUG) +console_formatter = logging.Formatter( + '%(name)s - %(levelname)s - %(message)s') +console_handler.setFormatter(console_formatter) +logger.addHandler(console_handler) + +def log_decorator(func): + """ + Decorator to log functions while in debug mode, i.e., when + `debug_montecarlo` is True in the config. Works for + `@jit'd and `@njit`'d functions, but with a significant speed + penalty. + + TODO: in nopython mode: do I need a context manager? + + Input: + func : (function) function to be logged. + + Output: + wrapper : (function) wrapper to the function being logged. + """ + + # to ensure that calling `help` on the decorated function still works + # @wraps(func) + def wrapper(*args, **kwargs): + """ + Wrapper to the log_decorator. + + When called, it has the side effect of + logging every `BUFFER` input and output to `func`, if `DEBUG_MODE` is + `True`. + + Input: + *args : arguments to be passed to `func`. + **kwargs : keyword arguments to be passed to `func`. + + Output: + result : result of calling `func` on `*args` and `**kwargs`. + """ + result = func(*args, **kwargs) + if DEBUG_MODE: + global ticker + ticker += 1 + if ticker % BUFFER == 0: # without a buffer, performance suffers + logger.debug(f'Func: {func.__name__}. Input: {(args, kwargs)}') + logger.debug(f'Output: {result}.') + return result + + return wrapper diff --git a/tardis/montecarlo/montecarlo_numba/numba_interface.py b/tardis/montecarlo/montecarlo_numba/numba_interface.py index 82c7631c58e..d52fb377deb 100644 --- a/tardis/montecarlo/montecarlo_numba/numba_interface.py +++ b/tardis/montecarlo/montecarlo_numba/numba_interface.py @@ -5,6 +5,9 @@ from tardis import constants as const +from tardis.montecarlo import montecarlo_configuration as montecarlo_configuration + + C_SPEED_OF_LIGHT = const.c.to('cm/s').value @@ -136,52 +139,36 @@ def __init__(self, spectrum_frequency, number_of_vpackets, estimators_spec = [ ('j_estimator', float64[:]), ('nu_bar_estimator', float64[:]), - ('j_b_lu_estimator', float64[:, :]), - ('edot_lu_estimator', float64[:, :]) + ('j_blue_estimator', float64[:, :]), + ('Edotlu_estimator', float64[:, :]) ] @jitclass(estimators_spec) class Estimators(object): - def __init__(self, j_estimator, nu_bar_estimator, j_b_lu_estimator, - edot_lu_estimator): + def __init__(self, j_estimator, nu_bar_estimator, j_blue_estimator, + Edotlu_estimator): self.j_estimator = j_estimator self.nu_bar_estimator = nu_bar_estimator - self.j_b_lu_estimator = j_b_lu_estimator - self.edot_lu_estimator = edot_lu_estimator - -monte_carlo_configuration_spec = [ - ('line_interaction_type', int64), - ('number_of_vpackets', int64), - ('temporary_v_packet_bins', int64), - ('full_relativity', boolean) -] - - -@jitclass(monte_carlo_configuration_spec) -class MonteCarloConfiguration(object): - def __init__(self, number_of_vpackets, line_interaction_type, - temporary_v_packet_bins, full_relativity=False): - self.line_interaction_type = line_interaction_type - self.number_of_vpackets = number_of_vpackets - self.temporary_v_packet_bins = temporary_v_packet_bins - self.full_relativity = full_relativity + self.j_blue_estimator = j_blue_estimator + self.Edotlu_estimator = Edotlu_estimator def configuration_initialize(runner, number_of_vpackets, temporary_v_packet_bins=20000): if runner.line_interaction_type == 'macroatom': - line_interaction_type = LineInteractionType.MACROATOM + montecarlo_configuration.line_interaction_type = LineInteractionType.MACROATOM elif runner.line_interaction_type == 'downbranch': - line_interaction_type = LineInteractionType.DOWNBRANCH + montecarlo_configuration.line_interaction_type = LineInteractionType.DOWNBRANCH elif runner.line_interaction_type == 'scatter': - line_interaction_type = LineInteractionType.SCATTER + montecarlo_configuration.line_interaction_type = LineInteractionType.SCATTER else: raise ValueError(f'Line interaction type must be one of "macroatom",' f'"downbranch", or "scatter" but is ' f'{runner.line_interaction_type}') - - return MonteCarloConfiguration(number_of_vpackets, line_interaction_type, - temporary_v_packet_bins) + montecarlo_configuration.number_of_vpackets = number_of_vpackets + montecarlo_configuration.temporary_v_packet_bins = temporary_v_packet_bins + montecarlo_configuration.full_relativity = runner.enable_full_relativity + montecarlo_configuration.single_packet_seed = runner.single_packet_seed #class TrackRPacket(object): diff --git a/tardis/montecarlo/montecarlo_numba/r_packet.py b/tardis/montecarlo/montecarlo_numba/r_packet.py index d9fa949096a..70fcb0bed50 100644 --- a/tardis/montecarlo/montecarlo_numba/r_packet.py +++ b/tardis/montecarlo/montecarlo_numba/r_packet.py @@ -5,6 +5,8 @@ from tardis.montecarlo.montecarlo_numba import njit_dict +from tardis.montecarlo import montecarlo_configuration as montecarlo_configuration +from tardis.montecarlo.montecarlo_numba.montecarlo_logger import log_decorator from tardis import constants as const class MonteCarloException(ValueError): @@ -62,9 +64,10 @@ def calculate_distance_boundary(r, mu, r_inner, r_outer): return distance, delta_shell + +# @log_decorator @njit(**njit_dict) -def calculate_distance_line(r_packet, comov_nu, nu_line, time_explosion, - montecarlo_configuration): +def calculate_distance_line(r_packet, comov_nu, nu_line, time_explosion): """ Parameters @@ -84,27 +87,37 @@ def calculate_distance_line(r_packet, comov_nu, nu_line, time_explosion, if nu_line == 0.0: return MISS_DISTANCE - nu_diff = comov_nu - nu_line + + # for numerical reasons, if line is too close, we set the distance to 0. if np.abs(nu_diff / comov_nu) < CLOSE_LINE_THRESHOLD: nu_diff = 0.0 - if nu_diff <= 0: - print('nu difference is less than 0.0', nu_diff, comov_nu, nu, nu_line, time_explosion) - raise MonteCarloException('nu difference is less than 0.0') - if montecarlo_configuration.full_relativity: - nu_r = nu_line / nu - ct = C_SPEED_OF_LIGHT * time_explosion - distance = -r_packet.mu * r_packet.r + ( - ct - nu_r**2 * np.sqrt( - ct**2 - (1 + r_packet**2 * (1 - r_packet.mu**2) * - (1 + 1 / nu_r**2)))) / (1 + nu_r**3) - else: + if nu_diff >= 0: distance = (nu_diff / nu) * C_SPEED_OF_LIGHT * time_explosion + else: + raise MonteCarloException('nu difference is less than 0.0; for more' + ' information, see print statement beforehand') + if montecarlo_configuration.full_relativity: + return calculate_distance_line_full_relativity(nu_line, nu, + time_explosion, + r_packet) return distance +@njit(**njit_dict) +def calculate_distance_line_full_relativity(nu_line, nu, time_explosion, + r_packet): + # distance = - mu * r + (ct - nu_r * nu_r * sqrt(ct * ct - (1 + r * r * (1 - mu * mu) * (1 + pow(nu_r, -2))))) / (1 + nu_r * nu_r); + nu_r = nu_line / nu + ct = C_SPEED_OF_LIGHT * time_explosion + distance = -r_packet.mu * r_packet.r + ( + ct - nu_r ** 2 * np.sqrt( + ct ** 2 - (1 + r_packet.r ** 2 * (1 - r_packet.mu ** 2) * + (1 + pow(nu_r, -2))))) / (1 + nu_r ** 2) + return distance + @njit(**njit_dict) def calculate_distance_electron(electron_density, tau_event): return tau_event / (electron_density * SIGMA_THOMSON) @@ -115,9 +128,37 @@ def calculate_tau_electron(electron_density, distance): @njit(**njit_dict) def get_doppler_factor(r, mu, time_explosion): - beta = (r / time_explosion) / C_SPEED_OF_LIGHT + beta = r / (time_explosion * C_SPEED_OF_LIGHT) + if not montecarlo_configuration.full_relativity: + return get_doppler_factor_partial_relativity(mu, beta) + else: + return get_doppler_factor_full_relativity(mu, beta) + +@njit(**njit_dict) +def get_doppler_factor_partial_relativity(mu, beta): return 1.0 - mu * beta +@njit(**njit_dict) +def get_doppler_factor_full_relativity(mu, beta): + return (1.0 - mu * beta) / np.sqrt(1 - beta * beta) + + +@njit(**njit_dict) +def get_inverse_doppler_factor(r, mu, time_explosion): + beta = (r / time_explosion) / C_SPEED_OF_LIGHT + if not montecarlo_configuration.full_relativity: + return get_inverse_doppler_factor_partial_relativity(mu, beta) + else: + return get_inverse_doppler_factor_full_relativity(mu, beta) + +@njit(**njit_dict) +def get_inverse_doppler_factor_partial_relativity(mu, beta): + return 1.0 / (1.0 - mu * beta) + +@njit(**njit_dict) +def get_inverse_doppler_factor_full_relativity(mu, beta): + return (1.0 + mu * beta) / np.sqrt(1 - beta * beta) + @njit(**njit_dict) def get_random_mu(): return 2.0 * np.random.random() - 1.0 @@ -144,7 +185,7 @@ def initialize_line_id(self, numba_plasma, numba_model): @njit(**njit_dict) def update_line_estimators(estimators, r_packet, cur_line_id, distance_trace, - time_explosion, montecarlo_configuration): + time_explosion): """ Function to update the line estimators @@ -167,21 +208,30 @@ def update_line_estimators(estimators, r_packet, cur_line_id, distance_trace, """ if not montecarlo_configuration.full_relativity: - doppler_factor = 1.0 - ((distance_trace + r_packet.mu * r_packet.r) / - (time_explosion * C_SPEED_OF_LIGHT)) - energy = r_packet.energy * doppler_factor + energy = calc_packet_energy(r_packet, distance_trace, + time_explosion) else: - # accurate to 1 / gamma - according to C. Vogl - energy = r_packet.energy + energy = calc_packet_energy_full_relativity(r_packet) - estimators.j_b_lu_estimator[cur_line_id, r_packet.current_shell_id] += ( + estimators.j_blue_estimator[cur_line_id, r_packet.current_shell_id] += ( energy / r_packet.nu) - estimators.edot_lu_estimator[cur_line_id, r_packet.current_shell_id] += ( + estimators.Edotlu_estimator[cur_line_id, r_packet.current_shell_id] += ( energy) @njit(**njit_dict) -def trace_packet(r_packet, numba_model, numba_plasma, estimators, - montecarlo_configuration): +def calc_packet_energy_full_relativity(r_packet): + # accurate to 1 / gamma - according to C. Vogl + return r_packet.energy + +@njit(**njit_dict) +def calc_packet_energy(r_packet, distance_trace, time_explosion): + doppler_factor = 1.0 - ((distance_trace + r_packet.mu * r_packet.r) / + (time_explosion * C_SPEED_OF_LIGHT)) + energy = r_packet.energy * doppler_factor + return energy + +@njit(**njit_dict) +def trace_packet(r_packet, numba_model, numba_plasma, estimators): """ Parameters @@ -189,7 +239,6 @@ def trace_packet(r_packet, numba_model, numba_plasma, estimators, numba_model: tardis.montecarlo.montecarlo_numba.numba_interface.NumbaModel numba_plasma: tardis.montecarlo.montecarlo_numba.numba_interface.NumbaPlasma estimators: tardis.montecarlo.montecarlo_numba.numba_interface.Estimators - montecarlo_configuration: tardis.montecarlo.montecarlo_numba.numba_interface.MonteCarloConfiguration Returns ------- @@ -239,8 +288,7 @@ def trace_packet(r_packet, numba_model, numba_plasma, estimators, # Calculating the distance until the current photons co-moving nu # redshifts to the line frequency distance_trace = calculate_distance_line( - r_packet, comov_nu, nu_line, numba_model.time_explosion, - montecarlo_configuration) + r_packet, comov_nu, nu_line, numba_model.time_explosion) # calculating the tau electron of how far the trace has progressed tau_trace_electron = calculate_tau_electron(cur_electron_density, @@ -269,7 +317,7 @@ def trace_packet(r_packet, numba_model, numba_plasma, estimators, update_line_estimators( estimators, r_packet, cur_line_id, distance_trace, - numba_model.time_explosion, montecarlo_configuration) + numba_model.time_explosion) if tau_trace_combined > tau_event: interaction_type = InteractionType.LINE # Line @@ -301,8 +349,7 @@ def trace_packet(r_packet, numba_model, numba_plasma, estimators, @njit(**njit_dict) -def move_r_packet(r_packet, distance, time_explosion, numba_estimator, - montecarlo_configuration): +def move_r_packet(r_packet, distance, time_explosion, numba_estimator): """Move packet a distance and recalculate the new angle mu Parameters @@ -318,23 +365,26 @@ def move_r_packet(r_packet, distance, time_explosion, numba_estimator, distance in cm """ - - doppler_factor = get_doppler_factor(r_packet.r, r_packet.mu, time_explosion) + doppler_factor = get_doppler_factor(r_packet.r, + r_packet.mu, + time_explosion) comov_nu = r_packet.nu * doppler_factor comov_energy = r_packet.energy * doppler_factor if montecarlo_configuration.full_relativity: - numba_estimator.j_estimator[r_packet.current_shell_id] += ( - comov_energy * distance * doppler_factor) - numba_estimator.nu_bar_estimator[r_packet.current_shell_id] += ( - comov_energy * distance * comov_nu * doppler_factor) - + distance = distance * doppler_factor + set_estimators_full_relativity(r_packet, + distance, + numba_estimator, + comov_nu, + comov_energy, + doppler_factor) else: - numba_estimator.j_estimator[r_packet.current_shell_id] += ( - comov_energy * distance) - numba_estimator.nu_bar_estimator[r_packet.current_shell_id] += ( - comov_energy * distance * comov_nu) - + set_estimators(r_packet, + distance, + numba_estimator, + comov_nu, + comov_energy) r = r_packet.r if (distance > 0.0): new_r = np.sqrt(r**2 + distance**2 + @@ -342,6 +392,52 @@ def move_r_packet(r_packet, distance, time_explosion, numba_estimator, r_packet.mu = (r_packet.mu * r + distance) / new_r r_packet.r = new_r +@njit(**njit_dict) +def set_estimators(r_packet, distance, numba_estimator, comov_nu, comov_energy): + numba_estimator.j_estimator[r_packet.current_shell_id] += ( + comov_energy * distance) + numba_estimator.nu_bar_estimator[r_packet.current_shell_id] += ( + comov_energy * distance * comov_nu) + +@njit(**njit_dict) +def set_estimators_full_relativity(r_packet, + distance, + numba_estimator, + comov_nu, + comov_energy, + doppler_factor): + numba_estimator.j_estimator[r_packet.current_shell_id] += ( + comov_energy * distance * doppler_factor) + numba_estimator.nu_bar_estimator[r_packet.current_shell_id] += ( + comov_energy * distance * comov_nu * doppler_factor) + +@njit(**njit_dict) +def line_emission(r_packet, emission_line_id, numba_plasma, time_explosion): + """ + + Parameters + ---------- + r_packet: tardis.montecarlo.montecarlo_numba.r_packet.RPacket + emission_line_id: int + numba_plasma + time_explosion + + Returns + ------- + + """ + doppler_factor = get_doppler_factor(r_packet.r, + r_packet.mu, + time_explosion,) + r_packet.nu = numba_plasma.line_list_nu[emission_line_id] / doppler_factor + r_packet.next_line_id = emission_line_id + 1 + if montecarlo_configuration.full_relativity: + r_packet.mu = angle_aberration_CMF_to_LF( + r_packet, + time_explosion + ) + + @njit(**njit_dict) def move_packet_across_shell_boundary(packet, delta_shell, no_of_shells): @@ -349,12 +445,12 @@ def move_packet_across_shell_boundary(packet, delta_shell, Move packet across shell boundary - realizing if we are still in the simulation or have moved out through the inner boundary or outer boundary and updating packet status. - + Parameters ---------- distance : float distance to move to shell boundary - + delta_shell: int is +1 if moving outward or -1 if moving inward @@ -370,25 +466,24 @@ def move_packet_across_shell_boundary(packet, delta_shell, else: packet.current_shell_id = next_shell_id - -def line_emission(r_packet, emission_line_id, numba_plasma, time_explosion): +@njit(**njit_dict) +def angle_aberration_CMF_to_LF(r_packet, time_explosion, mu): """ - - Parameters - ---------- - r_packet: tardis.montecarlo.montecarlo_numba.r_packet.RPacket - emission_line_id: int - numba_plasma - time_explosion - - Returns - ------- - + Converts angle aberration from comoving frame to + laboratory frame. """ - doppler_factor = get_doppler_factor(r_packet.r, r_packet.mu, time_explosion) - r_packet.nu = numba_plasma.line_list_nu[emission_line_id] / doppler_factor - r_packet.next_line_id = emission_line_id + 1 + ct = C_SPEED_OF_LIGHT * time_explosion + beta = r_packet.r / (ct) + return (r_packet.mu + beta) / (1.0 + beta * mu) -def angle_aberration_CMF_to_LF(r_packet, time_explosion): - beta = r_packet.r / (time_explosion * C_SPEED_OF_LIGHT) - return (r_packet.mu + beta) / (1.0 + beta * r_packet.mu) +@njit(**njit_dict) +def angle_aberration_LF_to_CMF(r_packet, time_explosion, mu): + """ + + c code: + double beta = rpacket_get_r (packet) * storage->inverse_time_explosion * INVERSE_C; + return (mu - beta) / (1.0 - beta * mu); + """ + ct = C_SPEED_OF_LIGHT * time_explosion + beta = r_packet.r /(ct) + return (mu - beta) / (1.0 - beta * mu) diff --git a/tardis/montecarlo/montecarlo_numba/single_packet_loop.py b/tardis/montecarlo/montecarlo_numba/single_packet_loop.py index 4955f306c39..e0c0a4b1996 100644 --- a/tardis/montecarlo/montecarlo_numba/single_packet_loop.py +++ b/tardis/montecarlo/montecarlo_numba/single_packet_loop.py @@ -2,20 +2,29 @@ import numpy as np from tardis.montecarlo.montecarlo_numba.r_packet import ( - InteractionType, PacketStatus, get_doppler_factor, trace_packet, - move_packet_across_shell_boundary, move_r_packet) + InteractionType, PacketStatus, get_inverse_doppler_factor, trace_packet, + move_packet_across_shell_boundary, move_r_packet, + MonteCarloException) from tardis.montecarlo.montecarlo_numba.interaction import ( general_scatter, line_scatter) from tardis.montecarlo.montecarlo_numba.numba_interface import \ LineInteractionType +from tardis.montecarlo import montecarlo_configuration as montecarlo_configuration from tardis.montecarlo.montecarlo_numba.vpacket import trace_vpacket_volley +from tardis import constants as const +C_SPEED_OF_LIGHT = const.c.to('cm/s').value + +from tardis.montecarlo.montecarlo_numba.montecarlo_logger import log_decorator +from tardis.montecarlo.montecarlo_numba import ( + montecarlo_logger as mc_logger) + +# @log_decorator @njit def single_packet_loop(r_packet, numba_model, numba_plasma, estimators, - vpacket_collection, montecarlo_configuration, - track_rpackets=False): + vpacket_collection): """ Parameters @@ -33,21 +42,20 @@ def single_packet_loop(r_packet, numba_model, numba_plasma, estimators, This function does not return anything but changes the r_packet object and if virtual packets are requested - also updates the vpacket_collection """ - - np.random.seed(r_packet.index) + flag = None line_interaction_type = montecarlo_configuration.line_interaction_type - doppler_factor = get_doppler_factor(r_packet.r, r_packet.mu, - numba_model.time_explosion) - r_packet.nu /= doppler_factor - r_packet.energy /= doppler_factor + if montecarlo_configuration.full_relativity: + set_packet_props_full_relativity(r_packet, numba_model) + else: + set_packet_props_partial_relativity(r_packet, numba_model) r_packet.initialize_line_id(numba_plasma, numba_model) trace_vpacket_volley(r_packet, vpacket_collection, numba_model, numba_plasma) - if track_rpackets: + if mc_logger.DEBUG_MODE: r_packet_track_nu = [r_packet.nu] r_packet_track_mu = [r_packet.mu] r_packet_track_r = [r_packet.r] @@ -55,8 +63,12 @@ def single_packet_loop(r_packet, numba_model, numba_plasma, estimators, r_packet_track_distance = [0.] while r_packet.status == PacketStatus.IN_PROCESS: + # try: distance, interaction_type, delta_shell = trace_packet( - r_packet, numba_model, numba_plasma, estimators=estimators) + r_packet, numba_model, numba_plasma, estimators) + # except MonteCarloException: + # flag = 'stop' + # break if interaction_type == InteractionType.BOUNDARY: move_r_packet(r_packet, distance, numba_model.time_explosion, @@ -69,9 +81,13 @@ def single_packet_loop(r_packet, numba_model, numba_plasma, estimators, estimators) line_scatter(r_packet, numba_model.time_explosion, line_interaction_type, numba_plasma) - + # try: trace_vpacket_volley( - r_packet, vpacket_collection, numba_model, numba_plasma) + r_packet, vpacket_collection, numba_model, numba_plasma, + ) + # except MonteCarloException: + # flag = 'stop' + # break elif interaction_type == InteractionType.ESCATTERING: move_r_packet(r_packet, distance, numba_model.time_explosion, @@ -80,8 +96,7 @@ def single_packet_loop(r_packet, numba_model, numba_plasma, estimators, trace_vpacket_volley(r_packet, vpacket_collection, numba_model, numba_plasma) - - if track_rpackets: + if mc_logger.DEBUG_MODE: r_packet_track_nu.append(r_packet.nu) r_packet_track_mu.append(r_packet.mu) r_packet_track_r.append(r_packet.r) @@ -89,7 +104,30 @@ def single_packet_loop(r_packet, numba_model, numba_plasma, estimators, r_packet_track_distance.append(distance) - if track_rpackets is True: + if mc_logger.DEBUG_MODE: return (r_packet_track_nu, r_packet_track_mu, r_packet_track_r, - r_packet_track_interaction, r_packet_track_distance) + r_packet_track_interaction, r_packet_track_distance, flag) + + # check where else initialize line ID happens! + +@njit +def set_packet_props_partial_relativity(r_packet, numba_model): + inverse_doppler_factor = get_inverse_doppler_factor(r_packet.r, r_packet.mu, + numba_model.time_explosion, + ) + r_packet.nu *= inverse_doppler_factor + r_packet.energy *= inverse_doppler_factor + +@njit +def set_packet_props_full_relativity(r_packet, numba_model): + beta = (r_packet.r / numba_model.time_explosion) / C_SPEED_OF_LIGHT + + inverse_doppler_factor = get_inverse_doppler_factor(r_packet.r, r_packet.mu, + numba_model.time_explosion, + ) + + r_packet.nu *= inverse_doppler_factor + r_packet.energy *= inverse_doppler_factor + r_packet.mu = (r_packet.mu + beta) / (1 + beta * r_packet.mu) + diff --git a/tardis/montecarlo/montecarlo_numba/vpacket.py b/tardis/montecarlo/montecarlo_numba/vpacket.py index 08f78a0ba68..9eb3efaa94b 100644 --- a/tardis/montecarlo/montecarlo_numba/vpacket.py +++ b/tardis/montecarlo/montecarlo_numba/vpacket.py @@ -2,12 +2,14 @@ from numba import jitclass, njit, gdb from tardis.montecarlo.montecarlo_numba import njit_dict +from tardis.montecarlo import montecarlo_configuration as montecarlo_configuration import numpy as np from tardis.montecarlo.montecarlo_numba.r_packet import ( calculate_distance_boundary, get_doppler_factor, calculate_distance_line, - calculate_tau_electron, PacketStatus, move_packet_across_shell_boundary) + calculate_tau_electron, PacketStatus, move_packet_across_shell_boundary, + angle_aberration_LF_to_CMF, angle_aberration_CMF_to_LF) vpacket_spec = [ ('r', float64), @@ -16,12 +18,14 @@ ('energy', float64), ('next_line_id', int64), ('current_shell_id', int64), - ('status', int64) + ('status', int64), + ('index', int64) ] @jitclass(vpacket_spec) class VPacket(object): - def __init__(self, r, mu, nu, energy, current_shell_id, next_line_id): + def __init__(self, r, mu, nu, energy, current_shell_id, next_line_id, + index=0): self.r = r self.mu = mu self.nu = nu @@ -29,6 +33,7 @@ def __init__(self, r, mu, nu, energy, current_shell_id, next_line_id): self.current_shell_id = current_shell_id self.next_line_id = next_line_id self.status = PacketStatus.IN_PROCESS + self.index = index @@ -66,7 +71,7 @@ def trace_vpacket_within_shell(v_packet, numba_model, numba_plasma): cur_line_id, v_packet.current_shell_id] distance_trace_line = calculate_distance_line( - v_packet.nu, comov_nu, nu_line, numba_model.time_explosion) + v_packet, comov_nu, nu_line, numba_model.time_explosion) if distance_boundary <= distance_trace_line: break @@ -118,7 +123,8 @@ def trace_vpacket(v_packet, numba_model, numba_plasma): return tau_trace_combined @njit(**njit_dict) -def trace_vpacket_volley(r_packet, vpacket_collection, numba_model, numba_plasma): +def trace_vpacket_volley(r_packet, vpacket_collection, numba_model, + numba_plasma): """ Shoot a volley of vpackets (the vpacket collection specifies how many) from the current position of the rpacket. @@ -140,7 +146,6 @@ def trace_vpacket_volley(r_packet, vpacket_collection, numba_model, numba_plasma return - no_of_vpackets = vpacket_collection.number_of_vpackets if no_of_vpackets == 0: @@ -150,6 +155,10 @@ def trace_vpacket_volley(r_packet, vpacket_collection, numba_model, numba_plasma if r_packet.r > numba_model.r_inner[0]: # not on inner_boundary mu_min = -np.sqrt(1 - (numba_model.r_inner[0] / r_packet.r) ** 2) v_packet_on_inner_boundary = False + if montecarlo_configuration.full_relativity: + mu_min = angle_aberration_LF_to_CMF(r_packet, + numba_model.time_explosion, + mu_min) else: v_packet_on_inner_boundary = True mu_min = 0.0 @@ -165,10 +174,18 @@ def trace_vpacket_volley(r_packet, vpacket_collection, numba_model, numba_plasma else: weight = (1 - mu_min) / (2 * no_of_vpackets) + # C code: next line, angle_aberration_CMF_to_LF( & virt_packet, storage); + if montecarlo_configuration.full_relativity: + v_packet_mu = angle_aberration_CMF_to_LF( + r_packet, + numba_model.time_explosion, + v_packet_mu + ) v_packet_doppler_factor = get_doppler_factor( r_packet.r, v_packet_mu, numba_model.time_explosion) # transform between r_packet mu and v_packet_mu + doppler_factor_ratio = ( r_packet_doppler_factor / v_packet_doppler_factor) @@ -177,8 +194,8 @@ def trace_vpacket_volley(r_packet, vpacket_collection, numba_model, numba_plasma v_packet = VPacket(r_packet.r, v_packet_mu, v_packet_nu, v_packet_energy, r_packet.current_shell_id, - r_packet.next_line_id) - + r_packet.next_line_id, i) + tau_vpacket = trace_vpacket(v_packet, numba_model, numba_plasma) v_packet.energy *= np.exp(-tau_vpacket) diff --git a/tardis/montecarlo/setup_package.py b/tardis/montecarlo/setup_package.py index 44dcc1913f9..2e0a5a107f9 100644 --- a/tardis/montecarlo/setup_package.py +++ b/tardis/montecarlo/setup_package.py @@ -3,7 +3,7 @@ from setuptools import Extension import os from astropy_helpers.distutils_helpers import get_distutils_option -from Cython.Build import cythonize +# from Cython.Build import cythonize import yaml @@ -30,27 +30,27 @@ yaml.dump(vpacket_config, open(vpacket_config_file_path, "w"), default_flow_style=False) -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'))] - 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'))] - - return cythonize( - Extension('tardis.montecarlo.montecarlo', sources, - include_dirs=['tardis/montecarlo/src', - 'tardis/montecarlo/src/randomkit', - 'numpy'], - depends=deps, - extra_compile_args=compile_args, - extra_link_args=link_args, - define_macros=define_macros) - ) +# 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'))] +# 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'))] + + # return cythonize( + # Extension('tardis.montecarlo.montecarlo', sources, + # include_dirs=['tardis/montecarlo/src', + # 'tardis/montecarlo/src/randomkit', + # 'numpy'], + # depends=deps, + # extra_compile_args=compile_args, + # extra_link_args=link_args, + # define_macros=define_macros) + # ) def get_package_data(): diff --git a/tardis/montecarlo/src/abbrev.h b/tardis/montecarlo/src/abbrev.h deleted file mode 100644 index fff7d07036c..00000000000 --- a/tardis/montecarlo/src/abbrev.h +++ /dev/null @@ -1,24 +0,0 @@ -#ifndef ABBREV_H -#define ABBREV_H - -#include - -/** - * @brief safe malloc; checks for NULL and aborts if encountered - */ -static inline void* safe_malloc(size_t size) { - void *mem = malloc(size); - if (mem == NULL && size != 0) abort(); - return mem; -} - -/** - * @brief safe realloc; checks for NULL and aborts if encountered - */ -static inline void* safe_realloc(void *ptr, size_t size) { - void *mem = realloc(ptr, size); - if (mem == NULL && size != 0) abort(); - return mem; -} - -#endif /* ABBREV_H */ diff --git a/tardis/montecarlo/src/cmontecarlo.c b/tardis/montecarlo/src/cmontecarlo.c deleted file mode 100644 index 41441d8d363..00000000000 --- a/tardis/montecarlo/src/cmontecarlo.c +++ /dev/null @@ -1,1291 +0,0 @@ - -#include -#include -#include -#include -#ifdef WITHOPENMP -#include -#endif -#include "io.h" -#include "abbrev.h" -#include "status.h" -#include "rpacket.h" -#include "cmontecarlo.h" - - -/** Look for a place to insert a value in an inversely sorted float array. - * - * @param x an inversely (largest to lowest) sorted float array - * @param x_insert a value to insert - * @param imin lower bound - * @param imax upper bound - * - * @return index of the next boundary to the left - */ -tardis_error_t -reverse_binary_search (const double *x, double x_insert, - int64_t imin, int64_t imax, int64_t * result) -{ - /* - Have in mind that *x points to a reverse sorted array. - That is large values will have small indices and small ones - will have large indices. - */ - tardis_error_t ret_val = TARDIS_ERROR_OK; - if (x_insert > x[imin] || x_insert < x[imax]) - { - ret_val = TARDIS_ERROR_BOUNDS_ERROR; - } - else - { - int 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 && x_insert < x[imin + 1]) - { - *result = imin + 1; - } - else - { - *result = imin; - } - } - return ret_val; -} - -/** Insert a value in to an array of line frequencies - * - * @param nu array of line frequencies - * @param nu_insert value of nu key - * @param number_of_lines number of lines in the line list - * - * @return index of the next line ot the red. If the key value is redder than the reddest line returns number_of_lines. - */ -tardis_error_t -line_search (const double *nu, double nu_insert, int64_t number_of_lines, - int64_t * result) -{ - tardis_error_t ret_val = TARDIS_ERROR_OK; - int64_t imin = 0; - int64_t imax = number_of_lines - 1; - if (nu_insert > nu[imin]) - { - *result = imin; - } - else if (nu_insert < nu[imax]) - { - *result = imax + 1; - } - else - { - ret_val = reverse_binary_search (nu, nu_insert, imin, imax, result); - *result = *result + 1; - } - return ret_val; -} - -tardis_error_t -binary_search (const double *x, double x_insert, int64_t imin, - int64_t imax, int64_t * result) -{ - /* - Have in mind that *x points to a sorted array. - Like [1,2,3,4,5,...] - */ - int imid; - tardis_error_t ret_val = TARDIS_ERROR_OK; - if (x_insert < x[imin] || x_insert > x[imax]) - { - ret_val = TARDIS_ERROR_BOUNDS_ERROR; - } - else - { - while (imax >= imin) - { - imid = (imin + imax) / 2; - if (x[imid] == x_insert) - { - *result = imid; - break; - } - else if (x[imid] < x_insert) - { - imin = imid + 1; - } - else - { - imax = imid - 1; - } - } - if (imax - imid == 2 && x_insert < x[imin + 1]) - { - *result = imin; - } - else - { - *result = imin; - } - } - return ret_val; -} - -void -angle_aberration_CMF_to_LF (rpacket_t *packet, const storage_model_t *storage) -{ - if (storage->full_relativity) - { - double beta = rpacket_get_r (packet) * storage->inverse_time_explosion * INVERSE_C; - double mu_0 = rpacket_get_mu (packet); - rpacket_set_mu (packet, (mu_0 + beta) / (1.0 + beta * mu_0)); - } -} - -/** Transform the lab frame direction cosine to the CMF - * - * @param packet - * @param storage - * @param mu lab frame direction cosine - * - * @return CMF direction cosine - */ -double -angle_aberration_LF_to_CMF (rpacket_t *packet, const storage_model_t *storage, double mu) -{ - double beta = rpacket_get_r (packet) * storage->inverse_time_explosion * INVERSE_C; - return (mu - beta) / (1.0 - beta * mu); -} - -double -rpacket_doppler_factor (const rpacket_t *packet, const storage_model_t *storage) -{ - double beta = rpacket_get_r (packet) * storage->inverse_time_explosion * INVERSE_C; - if (!storage->full_relativity) - { - return 1.0 - rpacket_get_mu (packet) * beta; - } - else - { - return (1.0 - rpacket_get_mu (packet) * beta) / sqrt (1 - beta * beta); - } -} - -double -rpacket_inverse_doppler_factor (const rpacket_t *packet, const storage_model_t *storage) -{ - double beta = rpacket_get_r (packet) * storage->inverse_time_explosion * INVERSE_C; - if (!storage->full_relativity) - { - return 1.0 / (1.0 - rpacket_get_mu (packet) * beta); - } - else - { - return (1.0 + rpacket_get_mu (packet) * beta) / sqrt (1 - beta * beta); - } -} - -double -bf_cross_section (const storage_model_t * storage, int64_t continuum_id, double comov_nu) -{ - double bf_xsect; - double *x_sect = storage->photo_xsect[continuum_id]->x_sect; - double *nu = storage->photo_xsect[continuum_id]->nu; - - switch (storage->bf_treatment) - { - case LIN_INTERPOLATION: - { - int64_t result; - tardis_error_t error = binary_search (nu, comov_nu, 0, - storage->photo_xsect[continuum_id]->no_of_points - 1, &result); - if (error == TARDIS_ERROR_BOUNDS_ERROR) - { - bf_xsect = 0.0; - } - else - { - bf_xsect = x_sect[result-1] + (comov_nu - nu[result-1]) / (nu[result] - nu[result-1]) - * (x_sect[result] - x_sect[result-1]); - } - break; - } - - case HYDROGENIC: - { - double nu_ratio = nu[0] / comov_nu; - bf_xsect = x_sect[0] * nu_ratio * nu_ratio * nu_ratio; - break; - } - - default: - fprintf (stderr, "(%d) is not a valid bound-free cross section treatment.\n", storage->bf_treatment); - exit(1); - } - return bf_xsect; -} - -void calculate_chi_bf (rpacket_t * packet, storage_model_t * storage) -{ - double doppler_factor = rpacket_doppler_factor (packet, storage); - double comov_nu = rpacket_get_nu (packet) * doppler_factor; - - int64_t no_of_continuum_edges = storage->no_of_edges; - int64_t current_continuum_id; - line_search(storage->continuum_list_nu, comov_nu, no_of_continuum_edges, ¤t_continuum_id); - rpacket_set_current_continuum_id (packet, current_continuum_id); - - int64_t shell_id = rpacket_get_current_shell_id (packet); - double T = storage->t_electrons[shell_id]; - double boltzmann_factor = exp (-(H * comov_nu) / (KB * T)); - - double bf_helper = 0; - for(int64_t i = current_continuum_id; i < no_of_continuum_edges; i++) - { - // get the level population for the level ijk in the current shell: - double l_pop = storage->l_pop[shell_id * no_of_continuum_edges + i]; - // get the level population ratio \frac{n_{0,j+1,k}}{n_{i,j,k}} \frac{n_{i,j,k}}{n_{0,j+1,k}}^{*}: - double l_pop_r = storage->l_pop_r[shell_id * no_of_continuum_edges + i]; - double bf_x_sect = bf_cross_section (storage, i, comov_nu); - if (bf_x_sect == 0.0) - { - break; - } - bf_helper += l_pop * bf_x_sect * (1.0 - l_pop_r * boltzmann_factor) * doppler_factor; - - packet->chi_bf_tmp_partial[i] = bf_helper; - } - - rpacket_set_chi_boundfree (packet, bf_helper); -} - -void calculate_chi_ff (rpacket_t * packet, const storage_model_t * storage) -{ - double doppler_factor = rpacket_doppler_factor (packet, storage); - double comov_nu = rpacket_get_nu (packet) * doppler_factor; - int64_t shell_id = rpacket_get_current_shell_id (packet); - double T = storage->t_electrons[shell_id]; - double boltzmann_factor = exp (-(H * comov_nu) / KB / T); - double chi_ff_factor = storage->chi_ff_factor[shell_id]; - - double chi_ff = chi_ff_factor * (1 - boltzmann_factor) * pow (comov_nu, -3); - - rpacket_set_chi_freefree (packet, chi_ff * doppler_factor); -} - -void -compute_distance2boundary (rpacket_t * packet, const storage_model_t * storage) -{ - double r = rpacket_get_r (packet); - double mu = rpacket_get_mu (packet); - double r_outer = storage->r_outer[rpacket_get_current_shell_id (packet)]; - double r_inner = storage->r_inner[rpacket_get_current_shell_id (packet)]; - double check, distance; - if (mu > 0.0) - { // direction outward - rpacket_set_next_shell_id (packet, 1); - distance = sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu); - } - else - { // going inward - if ( (check = r_inner * r_inner + (r * r * (mu * mu - 1.0)) )>= 0.0) - { // hit inner boundary - rpacket_set_next_shell_id (packet, -1); - distance = - r * mu - sqrt (check); - } - else - { // miss inner boundary - rpacket_set_next_shell_id (packet, 1); - distance = sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu); - } - } - rpacket_set_d_boundary (packet, distance); -} - -tardis_error_t -compute_distance2line (rpacket_t * packet, const storage_model_t * storage) -{ - if (!rpacket_get_last_line (packet)) - { - double r = rpacket_get_r (packet); - double mu = rpacket_get_mu (packet); - double nu = rpacket_get_nu (packet); - double nu_line = rpacket_get_nu_line (packet); - double distance, nu_diff; - double ct = storage->time_explosion * C; - double doppler_factor = rpacket_doppler_factor (packet, storage); - double comov_nu = nu * doppler_factor; - if ( (nu_diff = comov_nu - nu_line) >= 0) - { - if (!storage->full_relativity) - { - distance = (nu_diff / nu) * ct; - } - else - { - double nu_r = nu_line / nu; - distance = - mu * r + (ct - nu_r * nu_r * sqrt(ct * ct - - (1 + r * r * (1 - mu * mu) * (1 + pow (nu_r, -2))))) / (1 + nu_r * nu_r); - } - rpacket_set_d_line (packet, distance); - return TARDIS_ERROR_OK; - } - else - { - if (rpacket_get_next_line_id (packet) == storage->no_of_lines - 1) - { - fprintf (stderr, "last_line = %f\n", - storage-> - line_list_nu[rpacket_get_next_line_id (packet) - 1]); - fprintf (stderr, "Last line in line list reached!"); - } - else if (rpacket_get_next_line_id (packet) == 0) - { - fprintf (stderr, "First line in line list!"); - fprintf (stderr, "next_line = %f\n", - storage-> - line_list_nu[rpacket_get_next_line_id (packet) + 1]); - } - else - { - fprintf (stderr, "last_line = %f\n", - storage-> - line_list_nu[rpacket_get_next_line_id (packet) - 1]); - fprintf (stderr, "next_line = %f\n", - storage-> - line_list_nu[rpacket_get_next_line_id (packet) + 1]); - } - fprintf (stderr, "ERROR: Comoving nu less than nu_line!\n"); - fprintf (stderr, "comov_nu = %f\n", comov_nu); - fprintf (stderr, "nu_line = %f\n", nu_line); - fprintf (stderr, "(comov_nu - nu_line) / nu_line = %f\n", - (comov_nu - nu_line) / nu_line); - fprintf (stderr, "r = %f\n", r); - fprintf (stderr, "mu = %f\n", mu); - fprintf (stderr, "nu = %f\n", nu); - fprintf (stderr, "doppler_factor = %f\n", doppler_factor); - fprintf (stderr, "cur_zone_id = %" PRIi64 "\n", rpacket_get_current_shell_id (packet)); - return TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE; - } - } - else - { - rpacket_set_d_line (packet, MISS_DISTANCE); - return TARDIS_ERROR_OK; - } -} - -void -compute_distance2continuum (rpacket_t * packet, storage_model_t * storage) -{ - double chi_continuum, d_continuum; - double chi_electron = storage->electron_densities[rpacket_get_current_shell_id(packet)] * - storage->sigma_thomson; - if (storage->full_relativity) - { - chi_electron *= rpacket_doppler_factor (packet, storage); - } - - if (storage->cont_status == CONTINUUM_ON) - { - if (packet->compute_chi_bf) - { - calculate_chi_bf (packet, storage); - calculate_chi_ff (packet, storage); - } - else - { - packet->compute_chi_bf=true; - } - chi_continuum = rpacket_get_chi_boundfree (packet) + rpacket_get_chi_freefree (packet) + chi_electron; - d_continuum = rpacket_get_tau_event (packet) / chi_continuum; - } - else - { - chi_continuum = chi_electron; - d_continuum = storage->inverse_electron_densities[rpacket_get_current_shell_id (packet)] * - storage->inverse_sigma_thomson * rpacket_get_tau_event (packet); - } - - if (rpacket_get_virtual_packet(packet) > 0) - { - //Set all continuum distances to MISS_DISTANCE in case of an virtual_packet - d_continuum = MISS_DISTANCE; - packet->compute_chi_bf = false; - } - else - { - - // fprintf(stderr, "--------\n"); - // fprintf(stderr, "nu = %e \n", rpacket_get_nu(packet)); - // fprintf(stderr, "chi_electron = %e\n", chi_electron); - // fprintf(stderr, "chi_boundfree = %e\n", calculate_chi_bf(packet, storage)); - // fprintf(stderr, "chi_line = %e \n", rpacket_get_tau_event(packet) / rpacket_get_d_line(packet)); - // fprintf(stderr, "--------\n"); - - //rpacket_set_chi_freefree(packet, chi_freefree); - rpacket_set_chi_electron (packet, chi_electron); - } - rpacket_set_chi_continuum (packet, chi_continuum); - rpacket_set_d_continuum (packet, d_continuum); -} - -void -macro_atom (rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state) -{ - int emit = 0, i = 0, offset = -1; - uint64_t activate_level = rpacket_get_macro_atom_activation_level (packet); - while (emit >= 0) - { - double event_random = rk_double (mt_state); - i = storage->macro_block_references[activate_level] - 1; - double p = 0.0; - offset = storage->transition_probabilities_nd * - rpacket_get_current_shell_id (packet); - do - { - ++i; - p += storage->transition_probabilities[offset + i]; - } - while (p <= event_random); - emit = storage->transition_type[i]; - activate_level = storage->destination_level_id[i]; - } - switch (emit) - { - case BB_EMISSION: - line_emission (packet, storage, storage->transition_line_id[i], mt_state); - break; - - case BF_EMISSION: - rpacket_set_current_continuum_id (packet, storage->transition_line_id[i]); - storage->last_line_interaction_out_id[rpacket_get_id (packet)] = - rpacket_get_current_continuum_id (packet); - - continuum_emission (packet, storage, mt_state, sample_nu_free_bound, 3); - break; - - case FF_EMISSION: - continuum_emission (packet, storage, mt_state, sample_nu_free_free, 4); - break; - - case ADIABATIC_COOLING: - storage->last_interaction_type[rpacket_get_id (packet)] = 5; - rpacket_set_status (packet, TARDIS_PACKET_STATUS_REABSORBED); - break; - - default: - fprintf (stderr, "This process for macro-atom deactivation should not exist! (emit = %d)\n", emit); - exit(1); - } -} - -void -move_packet (rpacket_t * packet, storage_model_t * storage, double distance) -{ - double doppler_factor = rpacket_doppler_factor (packet, storage); - if (distance > 0.0) - { - double r = rpacket_get_r (packet); - double new_r = - sqrt (r * r + distance * distance + - 2.0 * r * distance * rpacket_get_mu (packet)); - rpacket_set_mu (packet, - (rpacket_get_mu (packet) * r + distance) / new_r); - rpacket_set_r (packet, new_r); - if (rpacket_get_virtual_packet (packet) <= 0) - { - double comov_energy = rpacket_get_energy (packet) * doppler_factor; - double comov_nu = rpacket_get_nu (packet) * doppler_factor; - if (storage->full_relativity) - { - distance *= doppler_factor; - } -#ifdef WITHOPENMP -#pragma omp atomic -#endif - storage->js[rpacket_get_current_shell_id (packet)] += - comov_energy * distance; -#ifdef WITHOPENMP -#pragma omp atomic -#endif - storage->nubars[rpacket_get_current_shell_id (packet)] += - comov_energy * distance * comov_nu; - - if (storage->cont_status) - { - increment_continuum_estimators(packet, storage, distance, comov_nu, comov_energy); - } - } - } -} - -void -increment_continuum_estimators (const rpacket_t * packet, storage_model_t * storage, double distance, - double comov_nu, double comov_energy) -{ - int64_t current_continuum_id; - int64_t no_of_continuum_edges = storage->no_of_edges; - int64_t shell_id = rpacket_get_current_shell_id (packet); - line_search(storage->continuum_list_nu, comov_nu, no_of_continuum_edges, ¤t_continuum_id); - double T = storage->t_electrons[shell_id]; - double boltzmann_factor = exp (-(H * comov_nu) / (KB * T)); - - #ifdef WITHOPENMP - #pragma omp atomic - #endif - storage->ff_heating_estimator[shell_id] += comov_energy * distance * rpacket_get_chi_freefree (packet); - - for(int64_t i = current_continuum_id; i < no_of_continuum_edges; i++) - { - double bf_xsect = bf_cross_section (storage, i, comov_nu); - int64_t photo_ion_idx = i * storage->no_of_shells + shell_id; - double photo_ion_estimator_helper = comov_energy * distance * bf_xsect / comov_nu; - double bf_heating_estimator_helper = - comov_energy * distance * bf_xsect * (1. - storage->continuum_list_nu[i] / comov_nu); - - #ifdef WITHOPENMP - #pragma omp atomic - #endif - storage->photo_ion_estimator[photo_ion_idx] += photo_ion_estimator_helper; - - #ifdef WITHOPENMP - #pragma omp atomic - #endif - storage->stim_recomb_estimator[photo_ion_idx] += photo_ion_estimator_helper * boltzmann_factor; - - #ifdef WITHOPENMP - #pragma omp atomic - #endif - storage->bf_heating_estimator[photo_ion_idx] += bf_heating_estimator_helper; - - #ifdef WITHOPENMP - #pragma omp atomic - #endif - storage->stim_recomb_cooling_estimator[photo_ion_idx] += bf_heating_estimator_helper * boltzmann_factor; - - if (photo_ion_estimator_helper != 0.0) - { - #ifdef WITHOPENMP - #pragma omp atomic - #endif - storage->photo_ion_estimator_statistics[photo_ion_idx] += 1; - } - else - { - break; - } - } -} - -double -get_increment_j_blue_estimator_energy (const rpacket_t * packet, - const storage_model_t * storage, - double d_line) -{ - double energy; - if (storage->full_relativity) - { - // Accurate up to a factor 1 / gamma - energy = rpacket_get_energy (packet); - } - else - { - double r = rpacket_get_r (packet); - double r_interaction = sqrt (r * r + d_line * d_line + - 2.0 * r * d_line * rpacket_get_mu (packet)); - double mu_interaction = (rpacket_get_mu (packet) * r + d_line) / r_interaction; - double doppler_factor = 1.0 - mu_interaction * r_interaction * - storage->inverse_time_explosion * INVERSE_C; - energy = rpacket_get_energy (packet) * doppler_factor; - } - return energy; -} - -void -increment_j_blue_estimator (const rpacket_t * packet, storage_model_t * storage, - double d_line, int64_t j_blue_idx) -{ - if (storage->line_lists_j_blues != NULL) - { - double energy = get_increment_j_blue_estimator_energy (packet, storage, - d_line); - #ifdef WITHOPENMP - #pragma omp atomic - #endif - storage->line_lists_j_blues[j_blue_idx] += - energy / rpacket_get_nu (packet); - } -} - -void -increment_Edotlu_estimator (const rpacket_t * packet, storage_model_t * storage, - double d_line, int64_t line_idx) -{ - if (storage->line_lists_Edotlu != NULL) - { - double energy = get_increment_j_blue_estimator_energy (packet, storage, - d_line); - #ifdef WITHOPENMP - #pragma omp atomic - #endif - storage->line_lists_Edotlu[line_idx] += energy; - } -} - - -int64_t -montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet, - int64_t virtual_mode, rk_state *mt_state) -{ - int64_t reabsorbed=-1; - if (virtual_mode == 0) - { - reabsorbed = montecarlo_one_packet_loop (storage, packet, 0, mt_state); - } - else - { - if ((rpacket_get_nu (packet) > storage->spectrum_virt_start_nu) && (rpacket_get_nu(packet) < storage->spectrum_virt_end_nu)) - { - for (int64_t i = 0; i < rpacket_get_virtual_packet_flag (packet); i++) - { - double weight; - rpacket_t virt_packet = *packet; - double mu_min; - if (rpacket_get_r(&virt_packet) > storage->r_inner[0]) - { - mu_min = - -1.0 * sqrt (1.0 - - (storage->r_inner[0] / rpacket_get_r(&virt_packet)) * - (storage->r_inner[0] / rpacket_get_r(&virt_packet))); - - if (storage->full_relativity) - { - // Need to transform the angular size of the photosphere into the CMF - mu_min = angle_aberration_LF_to_CMF (&virt_packet, storage, mu_min); - } - } - else - { - mu_min = 0.0; - } - double mu_bin = (1.0 - mu_min) / rpacket_get_virtual_packet_flag (packet); - rpacket_set_mu(&virt_packet,mu_min + (i + rk_double (mt_state)) * mu_bin); - switch (virtual_mode) - { - case -2: - weight = 1.0 / rpacket_get_virtual_packet_flag (packet); - break; - case -1: - weight = - 2.0 * rpacket_get_mu(&virt_packet) / - rpacket_get_virtual_packet_flag (packet); - break; - case 1: - weight = - (1.0 - - mu_min) / 2.0 / rpacket_get_virtual_packet_flag (packet); - break; - default: - fprintf (stderr, "Something has gone horribly wrong!\n"); - // FIXME MR: we need to somehow signal an error here - // I'm adding an exit() here to inform the compiler about the impossible path - exit(1); - } - angle_aberration_CMF_to_LF (&virt_packet, storage); - double doppler_factor_ratio = - rpacket_doppler_factor (packet, storage) / - rpacket_doppler_factor (&virt_packet, storage); - rpacket_set_energy(&virt_packet, - rpacket_get_energy (packet) * doppler_factor_ratio); - rpacket_set_nu(&virt_packet,rpacket_get_nu (packet) * doppler_factor_ratio); - reabsorbed = montecarlo_one_packet_loop (storage, &virt_packet, 1, mt_state); -#ifdef WITH_VPACKET_LOGGING -#ifdef WITHOPENMP -#pragma omp critical - { -#endif // WITHOPENMP - if (storage->virt_packet_count >= storage->virt_array_size) - { - storage->virt_array_size *= 2; - storage->virt_packet_nus = safe_realloc(storage->virt_packet_nus, sizeof(double) * storage->virt_array_size); - storage->virt_packet_energies = safe_realloc(storage->virt_packet_energies, sizeof(double) * storage->virt_array_size); - storage->virt_packet_last_interaction_in_nu = safe_realloc(storage->virt_packet_last_interaction_in_nu, sizeof(double) * storage->virt_array_size); - storage->virt_packet_last_interaction_type = safe_realloc(storage->virt_packet_last_interaction_type, sizeof(int64_t) * storage->virt_array_size); - storage->virt_packet_last_line_interaction_in_id = safe_realloc(storage->virt_packet_last_line_interaction_in_id, sizeof(int64_t) * storage->virt_array_size); - storage->virt_packet_last_line_interaction_out_id = safe_realloc(storage->virt_packet_last_line_interaction_out_id, sizeof(int64_t) * storage->virt_array_size); - } - storage->virt_packet_nus[storage->virt_packet_count] = rpacket_get_nu(&virt_packet); - storage->virt_packet_energies[storage->virt_packet_count] = rpacket_get_energy(&virt_packet) * weight; - storage->virt_packet_last_interaction_in_nu[storage->virt_packet_count] = storage->last_interaction_in_nu[rpacket_get_id (packet)]; - storage->virt_packet_last_interaction_type[storage->virt_packet_count] = storage->last_interaction_type[rpacket_get_id (packet)]; - storage->virt_packet_last_line_interaction_in_id[storage->virt_packet_count] = storage->last_line_interaction_in_id[rpacket_get_id (packet)]; - storage->virt_packet_last_line_interaction_out_id[storage->virt_packet_count] = storage->last_line_interaction_out_id[rpacket_get_id (packet)]; - storage->virt_packet_count += 1; -#ifdef WITHOPENMP - } -#endif // WITHOPENMP -#endif // WITH_VPACKET_LOGGING - if ((rpacket_get_nu(&virt_packet) < storage->spectrum_end_nu) && - (rpacket_get_nu(&virt_packet) > storage->spectrum_start_nu)) - { -#ifdef WITHOPENMP -#pragma omp critical - { -#endif // WITHOPENMP - int64_t virt_id_nu = - floor ((rpacket_get_nu(&virt_packet) - - storage->spectrum_start_nu) / - storage->spectrum_delta_nu); - storage->spectrum_virt_nu[virt_id_nu] += - rpacket_get_energy(&virt_packet) * weight; -#ifdef WITHOPENMP - } -#endif // WITHOPENMP - } - } - } - else - { - return 1; - } - } - return reabsorbed; -} - -void -move_packet_across_shell_boundary (rpacket_t * packet, - storage_model_t * storage, double distance, rk_state *mt_state) -{ - move_packet (packet, storage, distance); - if (rpacket_get_virtual_packet (packet) > 0) - { - double delta_tau_event = rpacket_get_chi_continuum(packet) * distance; - rpacket_set_tau_event (packet, - rpacket_get_tau_event (packet) + - delta_tau_event); - packet->compute_chi_bf = true; - } - else - { - rpacket_reset_tau_event (packet, mt_state); - } - if ((rpacket_get_current_shell_id (packet) < storage->no_of_shells - 1 - && rpacket_get_next_shell_id (packet) == 1) - || (rpacket_get_current_shell_id (packet) > 0 - && rpacket_get_next_shell_id (packet) == -1)) - { - rpacket_set_current_shell_id (packet, - rpacket_get_current_shell_id (packet) + - rpacket_get_next_shell_id (packet)); - } - else if (rpacket_get_next_shell_id (packet) == 1) - { - rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED); - } - else if ((storage->reflective_inner_boundary == 0) || - (rk_double (mt_state) > storage->inner_boundary_albedo)) - { - rpacket_set_status (packet, TARDIS_PACKET_STATUS_REABSORBED); - } - else - { - double doppler_factor = rpacket_doppler_factor (packet, storage); - double comov_nu = rpacket_get_nu (packet) * doppler_factor; - double comov_energy = rpacket_get_energy (packet) * doppler_factor; - // TODO: correct - rpacket_set_mu (packet, rk_double (mt_state)); - double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); - rpacket_set_nu (packet, comov_nu * inverse_doppler_factor); - rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); - if (rpacket_get_virtual_packet_flag (packet) > 0) - { - montecarlo_one_packet (storage, packet, -2, mt_state); - } - } -} - -void -montecarlo_thomson_scatter (rpacket_t * packet, storage_model_t * storage, - double distance, rk_state *mt_state) -{ - move_packet (packet, storage, distance); - double doppler_factor = rpacket_doppler_factor (packet, storage); - double comov_nu = rpacket_get_nu (packet) * doppler_factor; - double comov_energy = rpacket_get_energy (packet) * doppler_factor; - rpacket_set_mu (packet, 2.0 * rk_double (mt_state) - 1.0); - double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); - rpacket_set_nu (packet, comov_nu * inverse_doppler_factor); - rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); - rpacket_reset_tau_event (packet, mt_state); - storage->last_interaction_type[rpacket_get_id (packet)] = 1; - - angle_aberration_CMF_to_LF (packet, storage); - - if (rpacket_get_virtual_packet_flag (packet) > 0) - { - create_vpacket (storage, packet, mt_state); - } -} - -void -montecarlo_bound_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance, rk_state *mt_state) -{ - // current position in list of continuum edges -> indicates which bound-free processes are possible - int64_t ccontinuum = rpacket_get_current_continuum_id (packet); - - // Determine in which continuum the bf-absorption occurs - double chi_bf = rpacket_get_chi_boundfree (packet); - double zrand = rk_double (mt_state); - double zrand_x_chibf = zrand * chi_bf; - - while ((ccontinuum < storage->no_of_edges - 1) && (packet->chi_bf_tmp_partial[ccontinuum] <= zrand_x_chibf)) - { - ccontinuum++; - } - rpacket_set_current_continuum_id (packet, ccontinuum); - - /* For consistency reasons the branching between ionization and thermal energy is determined using the - comoving frequency at the initial position instead of the frequency at the point of interaction */ - double comov_nu = rpacket_get_nu (packet) * rpacket_doppler_factor (packet, storage); - - /* Move the packet to the place of absorption, select a direction for re-emission and impose energy conservation - in the co-moving frame. */ - move_packet (packet, storage, distance); - double old_doppler_factor = rpacket_doppler_factor (packet, storage); - rpacket_set_mu (packet, 2.0 * rk_double (mt_state) - 1.0); - double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); - double comov_energy = rpacket_get_energy (packet) * old_doppler_factor; - rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); - storage->last_interaction_type[rpacket_get_id (packet)] = 3; // last interaction was a bf-absorption - storage->last_line_interaction_in_id[rpacket_get_id (packet)] = ccontinuum; - - // Convert the rpacket to thermal or ionization energy - zrand = rk_double (mt_state); - int64_t activate_level = (zrand < storage->continuum_list_nu[ccontinuum] / comov_nu) ? - storage->cont_edge2macro_level[ccontinuum] : storage->kpacket2macro_level; - - rpacket_set_macro_atom_activation_level (packet, activate_level); - macro_atom (packet, storage, mt_state); -} - -void -montecarlo_free_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance, rk_state *mt_state) -{ - /* Move the packet to the place of absorption, select a direction for re-emission and impose energy conservation - in the co-moving frame. */ - move_packet (packet, storage, distance); - double old_doppler_factor = rpacket_doppler_factor (packet, storage); - rpacket_set_mu (packet, 2.0 * rk_double (mt_state) - 1.0); - double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); - double comov_energy = rpacket_get_energy (packet) * old_doppler_factor; - rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); - storage->last_interaction_type[rpacket_get_id (packet)] = 4; // last interaction was a ff-absorption - - // Create a k-packet - rpacket_set_macro_atom_activation_level (packet, storage->kpacket2macro_level); - macro_atom (packet, storage, mt_state); -} - -double -sample_nu_free_free (const rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state) -{ - int64_t shell_id = rpacket_get_current_shell_id (packet); - double T = storage->t_electrons[shell_id]; - double zrand = rk_double (mt_state); - return -KB * T / H * log(zrand); // Lucy 2003 MC II Eq.41 -} - -double -sample_nu_free_bound (const rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state) -{ - int64_t continuum_id = rpacket_get_current_continuum_id (packet); - double th_frequency = storage->continuum_list_nu[continuum_id]; - int64_t shell_id = rpacket_get_current_shell_id (packet); - double T = storage->t_electrons[shell_id]; - double zrand = rk_double (mt_state); - return th_frequency * (1 - (KB * T / H / th_frequency * log(zrand))); // Lucy 2003 MC II Eq.26 -} - -void -montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage, - double distance, rk_state *mt_state) -{ - uint64_t next_line_id = rpacket_get_next_line_id (packet); - uint64_t line2d_idx = next_line_id + - storage->no_of_lines * rpacket_get_current_shell_id (packet); - if (rpacket_get_virtual_packet (packet) == 0) - { - increment_j_blue_estimator (packet, storage, distance, line2d_idx); - increment_Edotlu_estimator (packet, storage, distance, line2d_idx); - } - double tau_line = - storage->line_lists_tau_sobolevs[line2d_idx]; - double tau_continuum = rpacket_get_chi_continuum(packet) * distance; - double tau_combined = tau_line + tau_continuum; - //rpacket_set_next_line_id (packet, rpacket_get_next_line_id (packet) + 1); - - if (next_line_id + 1 == storage->no_of_lines) - { - rpacket_set_last_line (packet, true); - } - if (rpacket_get_virtual_packet (packet) > 0) - { - rpacket_set_tau_event (packet, - rpacket_get_tau_event (packet) + tau_line); - rpacket_set_next_line_id (packet, next_line_id + 1); - test_for_close_line (packet, storage); - } - else if (rpacket_get_tau_event (packet) < tau_combined) - { // Line absorption occurs - move_packet (packet, storage, distance); - double old_doppler_factor = rpacket_doppler_factor (packet, storage); - rpacket_set_mu (packet, 2.0 * rk_double (mt_state) - 1.0); - double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); - double comov_energy = rpacket_get_energy (packet) * old_doppler_factor; - rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); - storage->last_interaction_in_nu[rpacket_get_id (packet)] = - rpacket_get_nu (packet); - storage->last_line_interaction_in_id[rpacket_get_id (packet)] = - next_line_id; - storage->last_line_interaction_shell_id[rpacket_get_id (packet)] = - rpacket_get_current_shell_id (packet); - storage->last_interaction_type[rpacket_get_id (packet)] = 2; - if (storage->line_interaction_id == 0) - { - line_emission (packet, storage, next_line_id, mt_state); - } - else if (storage->line_interaction_id >= 1) - { - rpacket_set_macro_atom_activation_level (packet, - storage->line2macro_level_upper[next_line_id]); - macro_atom (packet, storage, mt_state); - } - } - else - { // Packet passes line without interacting - rpacket_set_tau_event (packet, - rpacket_get_tau_event (packet) - tau_line); - rpacket_set_next_line_id (packet, next_line_id + 1); - packet->compute_chi_bf = false; - test_for_close_line (packet, storage); - } -} - -void -line_emission (rpacket_t * packet, storage_model_t * storage, int64_t emission_line_id, rk_state *mt_state) -{ - double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); - storage->last_line_interaction_out_id[rpacket_get_id (packet)] = emission_line_id; - if (storage->cont_status == CONTINUUM_ON) - { - storage->last_interaction_out_type[rpacket_get_id (packet)] = 2; - } - - rpacket_set_nu (packet, - storage->line_list_nu[emission_line_id] * inverse_doppler_factor); - rpacket_set_nu_line (packet, storage->line_list_nu[emission_line_id]); - rpacket_set_next_line_id (packet, emission_line_id + 1); - rpacket_reset_tau_event (packet, mt_state); - - angle_aberration_CMF_to_LF (packet, storage); - - if (rpacket_get_virtual_packet_flag (packet) > 0) - { - bool virtual_close_line = false; - if (!rpacket_get_last_line (packet) && - fabs (storage->line_list_nu[rpacket_get_next_line_id (packet)] - - rpacket_get_nu_line (packet)) < - (rpacket_get_nu_line (packet)* 1e-7)) - { - virtual_close_line = true; - } - // QUESTIONABLE!!! - bool old_close_line = rpacket_get_close_line (packet); - rpacket_set_close_line (packet, virtual_close_line); - create_vpacket (storage, packet, mt_state); - rpacket_set_close_line (packet, old_close_line); - virtual_close_line = false; - } - test_for_close_line (packet, storage); -} - -void test_for_close_line (rpacket_t * packet, const storage_model_t * storage) -{ - if (!rpacket_get_last_line (packet) && - fabs (storage->line_list_nu[rpacket_get_next_line_id (packet)] - - rpacket_get_nu_line (packet)) < (rpacket_get_nu_line (packet)* - 1e-7)) - { - rpacket_set_close_line (packet, true); - } -} - -void -continuum_emission (rpacket_t * packet, storage_model_t * storage, rk_state *mt_state, - pt2sample_nu sample_nu_continuum, int64_t emission_type_id) -{ - double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); - double nu_comov = sample_nu_continuum (packet, storage, mt_state); - rpacket_set_nu (packet, nu_comov * inverse_doppler_factor); - rpacket_reset_tau_event (packet, mt_state); - - storage->last_interaction_out_type[rpacket_get_id (packet)] = emission_type_id; - - // Have to find current position in line list - int64_t current_line_id; - line_search (storage->line_list_nu, nu_comov, storage->no_of_lines, ¤t_line_id); - - bool last_line = (current_line_id == storage->no_of_lines); - rpacket_set_last_line (packet, last_line); - rpacket_set_next_line_id (packet, current_line_id); - - angle_aberration_CMF_to_LF (packet, storage); - - if (rpacket_get_virtual_packet_flag (packet) > 0) - { - create_vpacket (storage, packet, mt_state); - } -} - -static void -montecarlo_compute_distances (rpacket_t * packet, storage_model_t * storage) -{ - // Check if the last line was the same nu as the current line. - if (rpacket_get_close_line (packet)) - { - // If so set the distance to the line to 0.0 - rpacket_set_d_line (packet, 0.0); - // Reset close_line. - rpacket_set_close_line (packet, false); - } - else - { - compute_distance2boundary (packet, storage); - compute_distance2line (packet, storage); - // FIXME MR: return status of compute_distance2line() is ignored - compute_distance2continuum (packet, storage); - } -} - -montecarlo_event_handler_t -get_event_handler (rpacket_t * packet, storage_model_t * storage, - double *distance, rk_state *mt_state) -{ - montecarlo_compute_distances (packet, storage); - double d_boundary = rpacket_get_d_boundary (packet); - double d_continuum = rpacket_get_d_continuum (packet); - double d_line = rpacket_get_d_line (packet); - montecarlo_event_handler_t handler; - if (d_line <= d_boundary && d_line <= d_continuum) - { - *distance = d_line; - handler = &montecarlo_line_scatter; - } - else if (d_boundary <= d_continuum) - { - *distance = d_boundary; - handler = &move_packet_across_shell_boundary; - } - else - { - *distance = d_continuum; - handler = montecarlo_continuum_event_handler (packet, storage, mt_state); - } - return handler; -} - -montecarlo_event_handler_t -montecarlo_continuum_event_handler (rpacket_t * packet, storage_model_t * storage, rk_state *mt_state) -{ - if (storage->cont_status) - { - double zrand_x_chi_cont = rk_double (mt_state) * rpacket_get_chi_continuum (packet); - double chi_th = rpacket_get_chi_electron (packet); - double chi_bf = rpacket_get_chi_boundfree (packet); - - if (zrand_x_chi_cont < chi_th) - { - return &montecarlo_thomson_scatter; - } - else if (zrand_x_chi_cont < chi_th + chi_bf) - { - return &montecarlo_bound_free_scatter; - } - else - { - return &montecarlo_free_free_scatter; - } - } - else - { - return &montecarlo_thomson_scatter; - } -} - -int64_t -montecarlo_one_packet_loop (storage_model_t * storage, rpacket_t * packet, - int64_t virtual_packet, rk_state *mt_state) -{ - rpacket_set_tau_event (packet, 0.0); - rpacket_set_nu_line (packet, 0.0); - rpacket_set_virtual_packet (packet, virtual_packet); - rpacket_set_status (packet, TARDIS_PACKET_STATUS_IN_PROCESS); - // Initializing tau_event if it's a real packet. - if (virtual_packet == 0) - { - rpacket_reset_tau_event (packet,mt_state); - } - // For a virtual packet tau_event is the sum of all the tau's that the packet passes. - while (rpacket_get_status (packet) == TARDIS_PACKET_STATUS_IN_PROCESS) - { - // Check if we are at the end of line list. - if (!rpacket_get_last_line (packet)) - { - rpacket_set_nu_line (packet, - storage-> - line_list_nu[rpacket_get_next_line_id - (packet)]); - } - double distance; - get_event_handler (packet, storage, &distance, mt_state) (packet, storage, - distance, mt_state); - if (virtual_packet > 0 && rpacket_get_tau_event (packet) > storage->tau_russian) - { - double event_random = rk_double (mt_state); - if (event_random > storage->survival_probability) - { - rpacket_set_energy(packet, 0.0); - rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED); - } - else - { - rpacket_set_energy(packet, - rpacket_get_energy (packet) / storage->survival_probability * - exp (-1.0 * rpacket_get_tau_event (packet))); - rpacket_set_tau_event (packet, 0.0); - } - } - } - if (virtual_packet > 0) - { - rpacket_set_energy (packet, - rpacket_get_energy (packet) * exp (-1.0 * - rpacket_get_tau_event - (packet))); - } - return rpacket_get_status (packet) == - TARDIS_PACKET_STATUS_REABSORBED ? 1 : 0; -} - -void -montecarlo_main_loop(storage_model_t * storage, int64_t virtual_packet_flag, int nthreads, unsigned long seed) -{ - int64_t finished_packets = 0; - storage->virt_packet_count = 0; -#ifdef WITH_VPACKET_LOGGING - storage->virt_packet_nus = (double *)safe_malloc(sizeof(double) * storage->no_of_packets); - storage->virt_packet_energies = (double *)safe_malloc(sizeof(double) * storage->no_of_packets); - storage->virt_packet_last_interaction_in_nu = (double *)safe_malloc(sizeof(double) * storage->no_of_packets); - storage->virt_packet_last_interaction_type = (int64_t *)safe_malloc(sizeof(int64_t) * storage->no_of_packets); - storage->virt_packet_last_line_interaction_in_id = (int64_t *)safe_malloc(sizeof(int64_t) * storage->no_of_packets); - storage->virt_packet_last_line_interaction_out_id = (int64_t *)safe_malloc(sizeof(int64_t) * storage->no_of_packets); - storage->virt_array_size = storage->no_of_packets; -#endif // WITH_VPACKET_LOGGING -#ifdef WITHOPENMP - omp_set_dynamic(0); - if (nthreads > 0) - { - omp_set_num_threads(nthreads); - } - -#pragma omp parallel firstprivate(finished_packets) - { - rk_state mt_state; - rk_seed (seed + omp_get_thread_num(), &mt_state); -#pragma omp master - { - fprintf(stderr, "Running with OpenMP - %d threads\n", omp_get_num_threads()); - print_progress(0, storage->no_of_packets); - } -#else - rk_state mt_state; - rk_seed (seed, &mt_state); - fprintf(stderr, "Running without OpenMP\n"); -#endif - int64_t chi_bf_tmp_size = (storage->cont_status) ? storage->no_of_edges : 0; - double *chi_bf_tmp_partial = safe_malloc(sizeof(double) * chi_bf_tmp_size); - - #pragma omp for - for (int64_t packet_index = 0; packet_index < storage->no_of_packets; ++packet_index) - { - int reabsorbed = 0; - rpacket_t packet; - rpacket_set_id(&packet, packet_index); - rpacket_init(&packet, storage, packet_index, virtual_packet_flag, chi_bf_tmp_partial); - if (virtual_packet_flag > 0) - { - reabsorbed = montecarlo_one_packet(storage, &packet, -1, &mt_state); - } - reabsorbed = montecarlo_one_packet(storage, &packet, 0, &mt_state); - storage->output_nus[packet_index] = rpacket_get_nu(&packet); - if (reabsorbed == 1) - { - storage->output_energies[packet_index] = -rpacket_get_energy(&packet); - } - else - { - storage->output_energies[packet_index] = rpacket_get_energy(&packet); - } - if ( ++finished_packets%100 == 0 ) - { -#ifdef WITHOPENMP - // WARNING: This only works with a static sheduler and gives an approximation of progress. - // The alternative would be to have a shared variable but that could potentially decrease performance when using many threads. - if (omp_get_thread_num() == 0 ) - print_progress(finished_packets * omp_get_num_threads(), storage->no_of_packets); -#else - print_progress(finished_packets, storage->no_of_packets); -#endif - } - } - free(chi_bf_tmp_partial); -#ifdef WITHOPENMP - } -#endif - print_progress(storage->no_of_packets, storage->no_of_packets); - fprintf(stderr,"\n"); -} - -void -create_vpacket (storage_model_t * storage, rpacket_t * packet, - rk_state *mt_state) -{ - if (storage->enable_biasing) - { - int64_t shell_id = rpacket_get_current_shell_id(packet); - double tau_bias = (storage->tau_bias[shell_id + 1] + - (storage->tau_bias[shell_id] - storage->tau_bias[shell_id + 1]) * - (storage->r_outer[shell_id] - rpacket_get_r (packet)) / - (storage->r_outer[shell_id] - storage->r_inner[shell_id])); - double vpacket_prob = exp(-tau_bias); - double event_random = rk_double (mt_state); - if (event_random < vpacket_prob) - { - packet->vpacket_weight = 1. / vpacket_prob; - montecarlo_one_packet (storage, packet, 1, mt_state); - } - } - else - { - montecarlo_one_packet (storage, packet, 1, mt_state); - } -} diff --git a/tardis/montecarlo/src/cmontecarlo.h b/tardis/montecarlo/src/cmontecarlo.h deleted file mode 100644 index 767017b7303..00000000000 --- a/tardis/montecarlo/src/cmontecarlo.h +++ /dev/null @@ -1,159 +0,0 @@ -#ifndef TARDIS_CMONTECARLO_H -#define TARDIS_CMONTECARLO_H - -#include -#include "randomkit/randomkit.h" -#include "rpacket.h" -#include "status.h" -#include "storage.h" - -#ifdef WITH_VPACKET_LOGGING -#define LOG_VPACKETS 1 -#else -#define LOG_VPACKETS 0 -#endif - -typedef void (*montecarlo_event_handler_t) (rpacket_t *packet, - storage_model_t *storage, - double distance, rk_state *mt_state); - -typedef double (*pt2sample_nu) (const rpacket_t *packet, - const storage_model_t *storage, - rk_state *mt_state); - -void initialize_random_kit (unsigned long seed); - -tardis_error_t line_search (const double *nu, double nu_insert, - int64_t number_of_lines, int64_t * result); - -tardis_error_t -reverse_binary_search (const double *x, double x_insert, - int64_t imin, int64_t imax, int64_t * result); - -tardis_error_t -binary_search (const double *x, double x_insert, - int64_t imin, int64_t imax, int64_t * result); - -double rpacket_doppler_factor (const rpacket_t *packet, const storage_model_t *storage); - -double rpacket_inverse_doppler_factor (const rpacket_t *packet, const storage_model_t *storage); - -void -angle_aberration_CMF_to_LF (rpacket_t * packet, const storage_model_t * storage); - -double -angle_aberration_LF_to_CMF (rpacket_t *packet, const storage_model_t *storage, double mu); - -/** Calculate the distance to shell boundary. - * - * @param packet rpacket structure with packet information - * @param storage storage model data - * - * @return distance to shell boundary - */ -void compute_distance2boundary (rpacket_t * packet, - const storage_model_t * storage); - -/** Calculate the distance the packet has to travel until it redshifts to the first spectral line. - * - * @param packet rpacket structure with packet information - * @param storage storage model data - * - * @return distance to the next spectral line - */ -tardis_error_t compute_distance2line (rpacket_t * packet, - const storage_model_t * storage); - -/** Calculate the distance to the next continuum event, which can be a Thomson scattering, bound-free absorption or - free-free transition. - * - * @param packet rpacket structure with packet information - * @param storage storage model data - * - * sets distance to the next continuum event (in centimeters) in packet rpacket structure - */ -void compute_distance2continuum (rpacket_t *packet, storage_model_t *storage); - -void macro_atom (rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state); - -void move_packet (rpacket_t * packet, storage_model_t * storage, - double distance); - -void increment_j_blue_estimator (const rpacket_t * packet, - storage_model_t * storage, - double d_line, - int64_t j_blue_idx); - -void increment_Edotlu_estimator (const rpacket_t * packet, - storage_model_t * storage, - double d_line, - int64_t j_blue_idx); - -double get_increment_j_blue_estimator_energy (const rpacket_t * packet, - const storage_model_t * storage, - double d_line); - -void -increment_continuum_estimators (const rpacket_t * packet, storage_model_t * storage, double distance, - double comov_nu, double comov_energy); - -int64_t montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet, - int64_t virtual_mode, rk_state *mt_state); - -int64_t montecarlo_one_packet_loop (storage_model_t * storage, - rpacket_t * packet, - int64_t virtual_packet, rk_state *mt_state); - -void montecarlo_main_loop(storage_model_t * storage, - int64_t virtual_packet_flag, - int nthreads, - unsigned long seed); - -montecarlo_event_handler_t get_event_handler (rpacket_t * packet, storage_model_t * storage, - double *distance, rk_state *mt_state); - -void test_for_close_line(rpacket_t * packet, const storage_model_t * storage); - -/* New handlers for continuum implementation */ - -montecarlo_event_handler_t montecarlo_continuum_event_handler(rpacket_t * packet, storage_model_t * storage, rk_state *mt_state); - -void montecarlo_free_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance, rk_state *mt_state); - -void montecarlo_bound_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance, rk_state *mt_state); - -double -bf_cross_section(const storage_model_t * storage, int64_t continuum_id, double comov_nu); - -void calculate_chi_bf(rpacket_t * packet, storage_model_t * storage); - -void calculate_chi_ff(rpacket_t * packet, const storage_model_t * storage); - -void -move_packet_across_shell_boundary (rpacket_t * packet, - storage_model_t * storage, double distance, rk_state *mt_state); - -void -montecarlo_thomson_scatter (rpacket_t * packet, storage_model_t * storage, - double distance, rk_state *mt_state); - -void -montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage, - double distance, rk_state *mt_state); - -void -line_emission(rpacket_t * packet, storage_model_t * storage, - int64_t emission_line_id, rk_state *mt_state); - -double sample_nu_free_free(const rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state); - -double sample_nu_free_bound(const rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state); - -void continuum_emission(rpacket_t * packet, storage_model_t * storage, rk_state *mt_state, - pt2sample_nu sample_nu_continuum, int64_t emission_type_id); - -void -create_vpacket (storage_model_t * storage, rpacket_t * packet, - rk_state *mt_state); - -#endif // TARDIS_CMONTECARLO_H diff --git a/tardis/montecarlo/src/integrator.c b/tardis/montecarlo/src/integrator.c deleted file mode 100644 index c885257d458..00000000000 --- a/tardis/montecarlo/src/integrator.c +++ /dev/null @@ -1,343 +0,0 @@ -#define _USE_MATH_DEFINES - -#include -#include -#include -#include -#include -#include - -#include "io.h" -#include "storage.h" -#include "integrator.h" -#include "cmontecarlo.h" - -#include "omp_helper.h" - -#define NULEN 0 -#define LINELEN 1 -#define PLEN 2 -#define SHELLEN 3 - -#define C_INV 3.33564e-11 -#define M_PI acos (-1) -#define KB_CGS 1.3806488e-16 -#define H_CGS 6.62606957e-27 - -/** - * Calculate the intensity of a black-body according to the following formula - * .. math:: - * I(\\nu, T) = \\frac{2h\\nu^3}{c^2}\frac{1}{e^{h\\nu \\beta_\\textrm{rad}} - 1} -*/ -double -intensity_black_body (double nu, double T) -{ - double beta_rad = 1 / (KB_CGS * T); - double coefficient = 2 * H_CGS * C_INV * C_INV; - return coefficient * nu * nu * nu / (exp(H_CGS * nu * beta_rad) - 1 ); -} - - -/*! @brief Algorithm to integrate an array using the trapezoid integration rule - * -*/ -double -trapezoid_integration (const double* array, const double h, int N) -{ - double result = (array[0] + array[N-1])/2; - for (int idx = 1; idx < N-1; ++idx) - { - result += array[idx]; - } - return result * h; -} - -/*! @brief Calculate distance to p line - * - * Calculate half of the length of the p-line inside a shell - * of radius r in terms of unit length (c * t_exp). - * If shell and p-line do not intersect, return 0. - * - * @param r radius of the shell - * @param p distance of the p-line to the center of the supernova - * @param inv_t inverse time_explosio is needed to norm to unit-length - * @return half the lenght inside the shell or zero - */ -static inline double -calculate_z(double r, double p, double inv_t) -{ - return (r > p) ? sqrt(r * r - p * p) * C_INV * inv_t : 0; -} - - -/*! - * @brief Calculate p line intersections - * - * This function calculates the intersection points of the p-line with each shell - * - * @param storage (INPUT) A storage model containing the environment - * @param p (INPUT) distance of the integration line to the center - * @param oz (OUTPUT) will be set with z values. The array is truncated by the - * value `1`. - * @param oshell_id (OUTPUT) will be set with the corresponding shell_ids - * @return number of shells intersected by the p-line - */ -int64_t -populate_z(const storage_model_t *storage, const double p, double *oz, int64_t *oshell_id) -{ - - // Abbreviations - double *r = storage->r_outer_i; - const int64_t N = storage->no_of_shells_i; - double inv_t = storage->inverse_time_explosion; - double z = 0; - - int64_t i = 0, offset = N, i_low, i_up; - - if (p <= storage->r_inner_i[0]) - { - // Intersect the photosphere - for(i = 0; i < N; ++i) - { // Loop from inside to outside - oz[i] = 1 - calculate_z(r[i], p, inv_t); - oshell_id[i] = i; - } - return N; - } - else - { - // No intersection with the photosphere - // that means we intersect each shell twice - for(i = 0; i < N; ++i) - { // Loop from inside to outside - z = calculate_z(r[i], p, inv_t); - if (z == 0) - continue; - if (offset == N) - { - offset = i; - } - // Calculate the index in the resulting array - i_low = N - i - 1; // the far intersection with the shell - i_up = N + i - 2 * offset; // the nearer intersection with the shell - - // Setting the arrays - oz[i_low] = 1 + z; - oshell_id[i_low] = i; - oz[i_up] = 1 - z; - oshell_id[i_up] = i; - } - return 2 * (N - offset); - } -} - - -/*! @brief Calculate integration points - * - */ -void -calculate_p_values(double R_max, int64_t N, double *opp) -{ - for(int i = 0; ino_of_lines, - size_shell = storage->no_of_shells_i, - size_tau = size_line * size_shell, - finished_nus = 0; - - double R_ph = storage->r_inner_i[0]; - double R_max = storage->r_outer_i[size_shell - 1]; - double pp[N]; - double *exp_tau = calloc(size_tau, sizeof(double)); -#pragma omp parallel firstprivate(L, exp_tau) - { - -#pragma omp master - { - if (omp_get_num_threads() > 1) { - fprintf(stderr, "Doing the formal integral\nRunning with OpenMP - %d threads\n", omp_get_num_threads()); - } else { - fprintf(stderr, "Doing the formal integral\nRunning without OpenMP\n"); - } - print_progress_fi(0, inu_size); - } - - // Initializing all the thread-local variables - int64_t offset = 0, i = 0, - size_z = 0, - idx_nu_start = 0, - direction = 0, - first = 0; - - double I_nu[N], - //I_nu_b[N], - //I_nu_r[N], - z[2 * storage->no_of_shells_i], - p = 0, - nu_start, - nu_end, - nu, - zstart, - zend, - escat_contrib, - escat_op, - Jkkp; - int64_t shell_id[2 * storage->no_of_shells_i]; - - double *pexp_tau, *patt_S_ul, *pline, *pJred_lu, *pJblue_lu; - - // Prepare exp_tau -#pragma omp for - for (i = 0; i < size_tau; ++i) { - exp_tau[i] = exp( -storage->line_lists_tau_sobolevs_i[i]); - } - calculate_p_values(R_max, N, pp); - // Done with the initialization - - // Loop over wavelengths in spectrum -#pragma omp for - for (int nu_idx = 0; nu_idx < inu_size ; ++nu_idx) - { - nu = inu[nu_idx]; - - // Loop over discrete values along line - for (int p_idx = 1; p_idx < N; ++p_idx) - { - escat_contrib = 0; - p = pp[p_idx]; - - // initialize z intersections for p values - size_z = populate_z(storage, p, z, shell_id); - - // initialize I_nu - if (p <= R_ph) - I_nu[p_idx] = intensity_black_body(nu * z[0], iT); - else - I_nu[p_idx] = 0; - - // Find first contributing line - nu_start = nu * z[0]; - nu_end = nu * z[1]; - line_search( - storage->line_list_nu, - nu_start, - size_line, - &idx_nu_start - ); - offset = shell_id[0] * size_line; - - // start tracking accumulated e-scattering optical depth - zstart = storage->time_explosion / C_INV * (1. - z[0]); - - // Initialize pointers - pline = storage->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; - - // TODO: Ugly loop - // Loop over all intersections - // TODO: replace by number of intersections and remove break - for (i = 0; i < size_z - 1; ++i) - { - escat_op = storage->electron_densities_i[shell_id[i]] * storage->sigma_thomson; - nu_end = nu * z[i+1]; - - // TODO: e-scattering: in principle we also have to check - // that dtau is <<1 (as assumed in Lucy 1999); if not, there - // is the chance that I_nu_b becomes negative - for (;pline < storage->line_list_nu + size_line; - // We have to increment all pointers simultaneously - ++pline, - ++pexp_tau, - ++patt_S_ul, - ++pJblue_lu) - { - if (*pline < nu_end) - { - // next resonance not in current shell - break; - } - - // Calculate e-scattering optical depth to next resonance point - zend = storage->time_explosion / C_INV * (1. - *pline / nu); - - if (first == 1){ - // First contribution to integration - // NOTE: this treatment of I_nu_b (given 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]) ; - first = 0; - } - else{ - // Account for e-scattering, c.f. Eqs 27, 28 in Lucy 1999 - Jkkp = 0.5 * (*pJred_lu + *pJblue_lu); - escat_contrib += (zend - zstart) * escat_op * (Jkkp - I_nu[p_idx]) ; - // this introduces the necessary offset of one element between pJblue_lu and pJred_lu - pJred_lu += 1; - } - I_nu[p_idx] = I_nu[p_idx] + escat_contrib; - - // Lucy 1999, Eq 26 - I_nu[p_idx] = I_nu[p_idx] * (*pexp_tau) + *patt_S_ul; - - // 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); - zend = storage->time_explosion / C_INV * (1. - nu_end / nu); - escat_contrib += (zend - zstart) * escat_op * (Jkkp - I_nu[p_idx]); - zstart = zend; - - if (i < size_z-1){ - // advance pointers - direction = shell_id[i+1] - shell_id[i]; - pexp_tau += direction * size_line; - patt_S_ul += direction * size_line; - pJred_lu += direction * size_line; - pJblue_lu += direction * size_line; - } - } - I_nu[p_idx] *= p; - } - // TODO: change integration to match the calculation of p values - L[nu_idx] = 8 * M_PI * M_PI * trapezoid_integration(I_nu, R_max/N, N); -#pragma omp atomic update - ++finished_nus; - if (finished_nus%10 == 0){ - print_progress_fi(finished_nus, inu_size); - } - } - // Free everything allocated on heap - free(exp_tau); - printf("\n"); - } - return L; -} \ No newline at end of file diff --git a/tardis/montecarlo/src/integrator.h b/tardis/montecarlo/src/integrator.h deleted file mode 100644 index c06e4a286a0..00000000000 --- a/tardis/montecarlo/src/integrator.h +++ /dev/null @@ -1,19 +0,0 @@ -#ifndef TARDIS_INTEGRATOR_H -#define TARDIS_INTEGRATOR_H - -#include "storage.h" - -double -integrate_intensity(const double* I_nu, const double h, int N); - -int64_t -populate_z(const storage_model_t *storage, const double p, double *oz, int64_t *oshell_id); - -void -calculate_p_values(double R_max, int64_t N, double *opp); - -double -*_formal_integral( - const storage_model_t *storage, double T, double *nu, int64_t nu_size, double *att_S_ul, double *Jred_lu, double *Jblue_lu, int N); - -#endif diff --git a/tardis/montecarlo/src/io.h b/tardis/montecarlo/src/io.h deleted file mode 100644 index 122af1e5967..00000000000 --- a/tardis/montecarlo/src/io.h +++ /dev/null @@ -1,32 +0,0 @@ -#define _POSIX_C_SOURCE 1 - -#include -#include -#include - -#define STATUS_FORMAT "\r\033[2K\t[%" PRId64 "%%] Packets(finished/total): %" PRId64 "/%" PRId64 -#define STATUS_FORMAT_FI "\r\033[2K\t[%" PRId64 "%%] Bins(finished/total): %" PRId64 "/%" PRId64 - -static inline void -print_progress (const int64_t current, const int64_t total) -{ - if (isatty(fileno(stderr))) - { - fprintf(stderr, STATUS_FORMAT, - current * 100 / total, - current, - total); - } -} - -static inline void -print_progress_fi (const int64_t current, const int64_t total) -{ - if (isatty(fileno(stderr))) - { - fprintf(stderr, STATUS_FORMAT_FI, - current * 100 / total, - current, - total); - } -} diff --git a/tardis/montecarlo/src/omp_helper.h b/tardis/montecarlo/src/omp_helper.h deleted file mode 100644 index abb9c2acc7c..00000000000 --- a/tardis/montecarlo/src/omp_helper.h +++ /dev/null @@ -1,7 +0,0 @@ -#ifdef WITHOPENMP - #include -#else -int omp_get_num_threads(){ - return 1; -} -#endif diff --git a/tardis/montecarlo/src/randomkit/LICENSE b/tardis/montecarlo/src/randomkit/LICENSE deleted file mode 100644 index c0d13eb1acc..00000000000 --- a/tardis/montecarlo/src/randomkit/LICENSE +++ /dev/null @@ -1,20 +0,0 @@ -Copyright (c) 2003-2006, Jean-Sebastien Roy (js@jeannot.org) - -Permission is hereby granted, free of charge, to any person obtaining a -copy of this software and associated documentation files (the -"Software"), to deal in the Software without restriction, including -without limitation the rights to use, copy, modify, merge, publish, -distribute, sublicense, and/or sell copies of the Software, and to -permit persons to whom the Software is furnished to do so, subject to -the following conditions: - -The above copyright notice and this permission notice shall be included -in all copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS -OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF -MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. -IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY -CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, -TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE -SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/tardis/montecarlo/src/randomkit/randomkit.h b/tardis/montecarlo/src/randomkit/randomkit.h deleted file mode 100644 index d5220d01857..00000000000 --- a/tardis/montecarlo/src/randomkit/randomkit.h +++ /dev/null @@ -1,42 +0,0 @@ -/* Random kit 1.6 */ - -/* - * Anyone who attempts to generate random numbers by deterministic means is, - * of course, living in a state of sin. - * -- John von Neumann - */ - -/* - * Copyright (c) 2003-2006, Jean-Sebastien Roy (js@jeannot.org) - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -/* @(#) $Jeannot: randomkit.h,v 1.29 2006/02/19 14:44:14 js Exp $ */ - -#ifndef _RANDOMKIT_ -#define _RANDOMKIT_ - -#include "rk_mt.h" -#include "rk_sobol.h" -#include "rk_primitive.h" -#include "rk_isaac.h" - -#endif /* _RANDOMKIT_ */ diff --git a/tardis/montecarlo/src/randomkit/rk_isaac.c b/tardis/montecarlo/src/randomkit/rk_isaac.c deleted file mode 100644 index 9f52fafae4c..00000000000 --- a/tardis/montecarlo/src/randomkit/rk_isaac.c +++ /dev/null @@ -1,258 +0,0 @@ -/* Random kit 1.6 */ - -/* - * Copyright (c) 2004-2006, Jean-Sebastien Roy (js@jeannot.org) - * - * ISAAC RNG by By Bob Jenkins. Based on Bob Jenkins public domain code. - * - * Original algorithm for the implementation of rk_interval function from - * Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by - * Magnus Jonsson. - * - * Constants used in the rk_double implementation by Isaku Wada. - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -static char const rcsid[] = - "@(#) $Jeannot: rk_isaac.c,v 1.4 2006/02/20 19:13:37 js Exp $"; - -#include "rk_isaac.h" -#include -#include -#include - -/* use the contents of randrsl[0..RK_ISAAC_STATE_LEN-1] as the seed. */ -static void isaac_init(rk_isaac_state *rg); - -void rk_isaac_seed(unsigned long seed, rk_isaac_state *state) -{ - rk_knuth_fill(seed, state->randrsl, RK_ISAAC_STATE_LEN); - isaac_init(state); -} - -rk_error rk_isaac_randomseed(rk_isaac_state *state) -{ - if(rk_devfill(state->randrsl, sizeof(state->randrsl), 1) == RK_NOERR) - { - isaac_init(state); - return RK_NOERR; - } - - rk_isaac_seed(rk_seedfromsystem(), state); - - return RK_ENODEV; -} - -#define ind(mm,x) ((mm)[(x>>2)&(RK_ISAAC_STATE_LEN-1)]) -#define rngstep(mix,a,b,mm,m,m2,r,x) \ -{ \ - x = *m; \ - a = ((a^(mix)) + *(m2++)) & 0xffffffff; \ - *(m++) = y = (ind(mm,x) + a + b) & 0xffffffff; \ - *(r++) = b = (ind(mm,y>>RK_ISAAC_STATE_POW) + x) & 0xffffffff; \ -} - -/* Call rand(rk_isaac_state *r) to retrieve a single 32-bit random value */ -unsigned long rk_isaac_random(rk_isaac_state *state) -{ - if (!state->randcnt--) - { - register unsigned long a,b,x,y,*m,*mm,*m2,*r,*mend; - mm=state->randmem; r=state->randrsl; - a = state->randa; b = (state->randb + (++state->randc)) & 0xffffffff; - for (m = mm, mend = m2 = m+(RK_ISAAC_STATE_LEN/2); m>6 , a, b, mm, m, m2, r, x); - rngstep( a<<2 , a, b, mm, m, m2, r, x); - rngstep( a>>16, a, b, mm, m, m2, r, x); - } - for (m2 = mm; m2>6 , a, b, mm, m, m2, r, x); - rngstep( a<<2 , a, b, mm, m, m2, r, x); - rngstep( a>>16, a, b, mm, m, m2, r, x); - } - state->randb = b; state->randa = a; - state->randcnt=RK_ISAAC_STATE_LEN-1; - } - return state->randrsl[state->randcnt] & 0xFFFFFFFF; -} - -#define mix(a,b,c,d,e,f,g,h) \ -{ \ - a^=b<<11; d+=a; b+=c; \ - b^=(c & 0xFFFFFFFF)>>2; e+=b; c+=d; \ - c^=d<<8; f+=c; d+=e; \ - d^=(e & 0xFFFFFFFF)>>16; g+=d; e+=f; \ - e^=f<<10; h+=e; f+=g; \ - f^=(g & 0xFFFFFFFF)>>4; a+=f; g+=h; \ - g^=h<<8; b+=g; h+=a; \ - h^=(a & 0xFFFFFFFF)>>9; c+=h; a+=b; \ -} - -/* if (flag==1), then use the contents of randrsl[] to initialize mm[]. */ -void isaac_init(rk_isaac_state *rk_isaac_state) -{ - int i; - unsigned long a,b,c,d,e,f,g,h; - unsigned long *m,*r; - rk_isaac_state->randa = rk_isaac_state->randb = rk_isaac_state->randc = 0; - m=rk_isaac_state->randmem; - r=rk_isaac_state->randrsl; - a=b=c=d=e=f=g=h=0x9e3779b9; /* the golden ratio */ - - for (i=0; i<4; ++i) /* scramble it */ - mix(a,b,c,d,e,f,g,h); - - /* initialize using the contents of r[] as the seed */ - for (i=0; irandcnt=0; - rk_isaac_state->has_gauss=0; -} - -long rk_isaac_long(rk_isaac_state *state) -{ - return rk_isaac_ulong(state) >> 1; -} - -unsigned long rk_isaac_ulong(rk_isaac_state *state) -{ -#if ULONG_MAX <= 0xffffffffUL - return rk_isaac_random(state); -#else - /* Assumes 64 bits */ - return (rk_isaac_random(state) << 32) | (rk_isaac_random(state)); -#endif -} - -unsigned long rk_isaac_interval(unsigned long max, rk_isaac_state *state) -{ - unsigned long mask = max, value; - - if (max == 0) return 0; - - /* Smallest bit mask >= max */ - mask |= mask >> 1; - mask |= mask >> 2; - mask |= mask >> 4; - mask |= mask >> 8; - mask |= mask >> 16; -#if ULONG_MAX > 0xffffffffUL - mask |= mask >> 32; -#endif - - /* Search a random value in [0..mask] <= max */ - while ((value = (rk_isaac_ulong(state) & mask)) > max); - - return value; -} - -double rk_isaac_double(rk_isaac_state *state) -{ - /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */ - long a = rk_isaac_random(state) >> 5, b = rk_isaac_random(state) >> 6; - return (a * 67108864.0 + b) / 9007199254740992.0; -} - -void rk_isaac_copy(rk_isaac_state *copy, rk_isaac_state *orig) -{ - memcpy(copy, orig, sizeof(rk_isaac_state)); -} - -double rk_isaac_gauss(rk_isaac_state *state) -{ - if (state->has_gauss) - { - state->has_gauss = 0; - return state->gauss; - } - else - { - double f, x1, x2, r2; - do - { - x1 = 2.0*rk_isaac_double(state) - 1.0; - x2 = 2.0*rk_isaac_double(state) - 1.0; - r2 = x1*x1 + x2*x2; - } - while (r2 >= 1.0 || r2 == 0.0); - - f = sqrt(-2.0*log(r2)/r2); /* Box-Muller transform */ - state->has_gauss = 1; - state->gauss = f*x1; /* Keep for next call */ - return f*x2; - } -} - -void rk_isaac_fill(void *buffer, size_t size, rk_isaac_state *state) -{ - unsigned long r; - unsigned char *buf = buffer; - rk_isaac_state tempstate; - - if (size > 0 && state == NULL) - { - rk_isaac_randomseed(&tempstate); - state = &tempstate; - } - - for (; size >= 4; size -= 4) - { - r = rk_isaac_random(state); - *(buf++) = r & 0xFF; - *(buf++) = (r >> 8) & 0xFF; - *(buf++) = (r >> 16) & 0xFF; - *(buf++) = (r >> 24) & 0xFF; - } - - if (!size) return; - - r = rk_isaac_random(state); - - for (; size; r >>= 8, size --) - *(buf++) = (unsigned char)(r & 0xFF); -} - -void rk_seed_isaac(rk_isaac_state *i_state, rk_state *state) -{ - rk_isaac_fill(state->key, sizeof(state->key), i_state); - state->pos = RK_STATE_LEN; - state->has_gauss = 0; -} diff --git a/tardis/montecarlo/src/randomkit/rk_isaac.h b/tardis/montecarlo/src/randomkit/rk_isaac.h deleted file mode 100644 index c8dfdcb451e..00000000000 --- a/tardis/montecarlo/src/randomkit/rk_isaac.h +++ /dev/null @@ -1,142 +0,0 @@ -/* Random kit 1.6 */ - -/* - * Copyright (c) 2004-2006, Jean-Sebastien Roy (js@jeannot.org) - * - * ISAAC RNG by By Bob Jenkins. Based on Bob Jenkins public domain code. - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -/* @(#) $Jeannot: rk_isaac.h,v 1.2 2006/02/19 14:40:26 js Exp $ */ - -/* - * Typical use: - * - * { - * rk_isaac_state state; - * unsigned long seed = 1, random_value; - * - * rk_isaac_seed(seed, &state); // Initialize the RNG - * ... - * random_value = rk_isaac_random(&state); // a random value in [0..RK_MAX] - * } - * - * Instead of rk_isaac_seed, you can use rk_isaac_randomseed which will get a - * random seed from /dev/random (or the clock, if /dev/random is unavailable): - * - * { - * rk_isaac_state state; - * unsigned long random_value; - * - * rk_isaac_randomseed(&state); // Initialize the RNG with a random seed - * ... - * random_value = rk_isaac_random(&state); // a random value in [0..RK_MAX] - * } - */ - -#ifndef _RK_ISAAC_ -#define _RK_ISAAC_ - -#include "rk_mt.h" - -#define RK_ISAAC_STATE_POW 8 -#define RK_ISAAC_STATE_LEN (1< -#include -#include -#include -#include -#include -#include -#include - -#ifdef _WIN32 -/* Windows */ -#include -#ifndef RK_NO_WINCRYPT -/* Windows crypto */ -#ifndef _WIN32_WINNT -#define _WIN32_WINNT 0x0400 -#endif -#include -#include -#endif -#else -/* Unix */ -#include -#include -#endif - -#include "randomkit.h" - -#ifndef RK_DEV_URANDOM -#define RK_DEV_URANDOM "/dev/urandom" -#endif - -#ifndef RK_DEV_RANDOM -#define RK_DEV_RANDOM "/dev/random" -#endif - -char *rk_strerror[RK_ERR_MAX] = -{ - "no error", - "random device unvavailable" -}; - -/* static functions */ -static unsigned long rk_hash(unsigned long key); - -void rk_knuth_fill(unsigned long seed, unsigned long *key, unsigned long len) -{ - int pos; - seed &= 0xffffffffUL; - - /* Knuth's PRNG as used in the Mersenne Twister reference implementation */ - for (pos=0; pos> 30)) + pos + 1) & 0xffffffffUL; - } -} - -void rk_seed(unsigned long seed, rk_state *state) -{ - rk_knuth_fill(seed, state->key, RK_STATE_LEN); - state->pos = RK_STATE_LEN; - state->has_gauss = 0; -} - -/* Thomas Wang 32 bits integer hash function */ -unsigned long rk_hash(unsigned long key) -{ - key += ~(key << 15); - key ^= (key >> 10); - key += (key << 3); - key ^= (key >> 6); - key += ~(key << 11); - key ^= (key >> 16); - return key; -} - -unsigned long rk_seedfromsystem() -{ -#ifndef _WIN32 - struct timeval tv; -#else - struct _timeb tv; -#endif - -#ifndef _WIN32 - gettimeofday(&tv, NULL); - return rk_hash(getpid()) ^ rk_hash(tv.tv_sec) ^ rk_hash(tv.tv_usec) - ^ rk_hash(clock()); -#else - _ftime(&tv); - return rk_hash(tv.time) ^ rk_hash(tv.millitm) ^ rk_hash(clock()); -#endif -} - -rk_error rk_randomseed(rk_state *state) -{ - if(rk_devfill(state->key, sizeof(state->key), 0) == RK_NOERR) - { - state->key[0] |= 0x80000000UL; /* ensures non-zero key */ - state->pos = RK_STATE_LEN; - state->has_gauss = 0; - return RK_NOERR; - } - - rk_seed(rk_seedfromsystem(), state); - - return RK_ENODEV; -} - -/* Magic Mersenne Twister constants */ -#define N 624 -#define M 397 -#define MATRIX_A 0x9908b0dfUL -#define UPPER_MASK 0x80000000UL -#define LOWER_MASK 0x7fffffffUL - -/* Slightly optimised reference implementation of the Mersenne Twister */ -unsigned long rk_random(rk_state *state) -{ - unsigned long y; - - if (state->pos == RK_STATE_LEN) - { - int i; - - for (i=0;ikey[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK); - state->key[i] = state->key[i+M] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); - } - for (;ikey[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK); - state->key[i] = state->key[i+(M-N)] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); - } - y = (state->key[N-1] & UPPER_MASK) | (state->key[0] & LOWER_MASK); - state->key[N-1] = state->key[M-1] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); - - state->pos = 0; - } - - y = state->key[state->pos++]; - - /* Tempering */ - y ^= (y >> 11); - y ^= (y << 7) & 0x9d2c5680UL; - y ^= (y << 15) & 0xefc60000UL; - y ^= (y >> 18); - - return y; -} - -long rk_long(rk_state *state) -{ - return rk_ulong(state) >> 1; -} - -unsigned long rk_ulong(rk_state *state) -{ -#if ULONG_MAX <= 0xffffffffUL - return rk_random(state); -#else - /* Assumes 64 bits */ - return (rk_random(state) << 32) | (rk_random(state)); -#endif -} - -unsigned long rk_interval(unsigned long max, rk_state *state) -{ - unsigned long mask = max, value; - - if (max == 0) return 0; - - /* Smallest bit mask >= max */ - mask |= mask >> 1; - mask |= mask >> 2; - mask |= mask >> 4; - mask |= mask >> 8; - mask |= mask >> 16; -#if ULONG_MAX > 0xffffffffUL - mask |= mask >> 32; -#endif - - /* Search a random value in [0..mask] <= max */ - while ((value = (rk_ulong(state) & mask)) > max); - - return value; -} - -double rk_double(rk_state *state) -{ - /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */ - long a = rk_random(state) >> 5, b = rk_random(state) >> 6; - return (a * 67108864.0 + b) / 9007199254740992.0; -} - -void rk_copy(rk_state *copy, rk_state *orig) -{ - memcpy(copy, orig, sizeof(rk_state)); -} - -void rk_fill(void *buffer, size_t size, rk_state *state) -{ - unsigned long r; - unsigned char *buf = buffer; - rk_state tempstate; - - if (size > 0 && state == NULL) - { - rk_randomseed(&tempstate); - state = &tempstate; - } - - for (; size >= 4; size -= 4) - { - r = rk_random(state); - *(buf++) = r & 0xFF; - *(buf++) = (r >> 8) & 0xFF; - *(buf++) = (r >> 16) & 0xFF; - *(buf++) = (r >> 24) & 0xFF; - } - - if (!size) return; - - r = rk_random(state); - - for (; size; r >>= 8, size --) - *(buf++) = (unsigned char)(r & 0xFF); -} - -rk_error rk_devfill(void *buffer, size_t size, int strong) -{ -#ifndef _WIN32 - FILE *rfile; - int done; - - if (strong) - rfile = fopen(RK_DEV_RANDOM, "rb"); - else - rfile = fopen(RK_DEV_URANDOM, "rb"); - if (rfile == NULL) - return RK_ENODEV; - done = fread(buffer, size, 1, rfile); - fclose(rfile); - if (done) - return RK_NOERR; -#else - -#ifndef RK_NO_WINCRYPT - HCRYPTPROV hCryptProv; - BOOL done; - - if (!CryptAcquireContext(&hCryptProv, NULL, NULL, PROV_RSA_FULL, - CRYPT_VERIFYCONTEXT) || !hCryptProv) - return RK_ENODEV; - done = CryptGenRandom(hCryptProv, size, (unsigned char *)buffer); - CryptReleaseContext(hCryptProv, 0); - if (done) - return RK_NOERR; -#endif - -#endif - - return RK_ENODEV; -} - -rk_error rk_altfill(void *buffer, size_t size, int strong, rk_state *state) -{ - rk_error err; - - err = rk_devfill(buffer, size, strong); - if (err) - rk_fill(buffer, size, state); - - return err; -} - -double rk_gauss(rk_state *state) -{ - if (state->has_gauss) - { - state->has_gauss = 0; - return state->gauss; - } - else - { - double f, x1, x2, r2; - do - { - x1 = 2.0*rk_double(state) - 1.0; - x2 = 2.0*rk_double(state) - 1.0; - r2 = x1*x1 + x2*x2; - } - while (r2 >= 1.0 || r2 == 0.0); - - f = sqrt(-2.0*log(r2)/r2); /* Box-Muller transform */ - state->has_gauss = 1; - state->gauss = f*x1; /* Keep for next call */ - return f*x2; - } -} diff --git a/tardis/montecarlo/src/randomkit/rk_mt.h b/tardis/montecarlo/src/randomkit/rk_mt.h deleted file mode 100644 index ff0e22a0879..00000000000 --- a/tardis/montecarlo/src/randomkit/rk_mt.h +++ /dev/null @@ -1,194 +0,0 @@ -/* Random kit 1.6 */ - -/* - * Anyone who attempts to generate random numbers by deterministic means is, - * of course, living in a state of sin. - * -- John von Neumann - */ - -/* - * Copyright (c) 2003-2006, Jean-Sebastien Roy (js@jeannot.org) - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -/* @(#) $Jeannot: rk_mt.h,v 1.6 2006/02/20 19:13:37 js Exp $ */ - -/* - * Typical use: - * - * { - * rk_state state; - * unsigned long seed = 1, random_value; - * - * rk_seed(seed, &state); // Initialize the RNG - * ... - * random_value = rk_random(&state); // a random value in [0..RK_MAX] - * } - * - * Instead of rk_seed, you can use rk_randomseed which will get a random seed - * from /dev/urandom (or the clock, if /dev/urandom is unavailable): - * - * { - * rk_state state; - * unsigned long random_value; - * - * rk_randomseed(&state); // Initialize the RNG with a random seed - * ... - * random_value = rk_random(&state); // a random value in [0..RK_MAX] - * } - */ - -/* - * Useful macro: - * RK_DEV_RANDOM: the device used for random seeding. - * defaults to "/dev/urandom" - */ - -#include - -#ifndef _RK_MT_ -#define _RK_MT_ - -#define RK_STATE_LEN 624 - -typedef struct rk_state_ -{ - unsigned long key[RK_STATE_LEN]; - int pos; - int has_gauss; /* !=0: gauss contains a gaussian deviate */ - double gauss; -} -rk_state; - -typedef enum { - RK_NOERR = 0, /* no error */ - RK_ENODEV = 1, /* no RK_DEV_RANDOM device */ - RK_ERR_MAX = 2 -} rk_error; - -/* error strings */ -extern char *rk_strerror[RK_ERR_MAX]; - -/* Maximum generated random value */ -#define RK_MAX 0xFFFFFFFFUL - -#ifdef __cplusplus -extern "C" { -#endif - -/* - * Initialize the RNG state using the given seed. - */ -extern void rk_seed(unsigned long seed, rk_state *state); - -/* - * Initialize the RNG state using a random seed. - * Uses /dev/urandom or, when unavailable, the clock (see rk_mt.c). - * Returns RK_NOERR when no errors occurs. - * Returns RK_ENODEV when the use of RK_DEV_RANDOM failed (for example because - * there is no such device). In this case, the RNG was initialized using the - * clock. - */ -extern rk_error rk_randomseed(rk_state *state); - -/* - * Returns a random unsigned long between 0 and RK_MAX inclusive - */ -extern unsigned long rk_random(rk_state *state); - -/* - * Returns a random long between 0 and LONG_MAX inclusive - */ -extern long rk_long(rk_state *state); - -/* - * Returns a random unsigned long between 0 and ULONG_MAX inclusive - */ -extern unsigned long rk_ulong(rk_state *state); - -/* - * Returns a random unsigned long between 0 and max inclusive. - */ -extern unsigned long rk_interval(unsigned long max, rk_state *state); - -/* - * Returns a random double between 0.0 and 1.0, 1.0 excluded. - */ -extern double rk_double(rk_state *state); - -/* - * Copy a random generator state. - */ -extern void rk_copy(rk_state *copy, rk_state *orig); - -/* - * fill the buffer with size random bytes. - * If state == NULL, the random generator is inilialized using rk_randomseed. - * Calling multiple times rk_randomseed should be avoided therefore calling - * multiple times rk_fill with state == NULL should be avoided. - */ -extern void rk_fill(void *buffer, size_t size, rk_state *state); - -/* - * fill the buffer with randombytes from the random device - * Returns RK_ENODEV if the device is unavailable, or RK_NOERR if it is - * On Unix, if strong is defined, RK_DEV_RANDOM is used. If not, RK_DEV_URANDOM - * is used instead. This parameter has no effect on Windows. - * Warning: on most unixes RK_DEV_RANDOM will wait for enough entropy to answer - * which can take a very long time on quiet systems. - */ -extern rk_error rk_devfill(void *buffer, size_t size, int strong); - -/* - * fill the buffer using rk_devfill if the random device is available and using - * rk_fill if is is not - * parameters have the same meaning as rk_fill and rk_devfill - * Returns RK_ENODEV if the device is unavailable, or RK_NOERR if it is - */ -extern rk_error rk_altfill(void *buffer, size_t size, int strong, - rk_state *state); - -/* - * return a random gaussian deviate with variance unity and zero mean. - */ -extern double rk_gauss(rk_state *state); - -/* Utility functions */ - -/* - * fill the key vector using Knuth RNG as used in MT reference implementation - * using the provided seed. The key vector length is len. - */ -extern void rk_knuth_fill(unsigned long seed, unsigned long *key, - unsigned long len); - -/* - * return a random unsigned long based upon the system state (clock, pid) - * used as seed when there is no random device. - */ -extern unsigned long rk_seedfromsystem(void); - - -#ifdef __cplusplus -} -#endif - -#endif /* _RK_MT_ */ diff --git a/tardis/montecarlo/src/randomkit/rk_primitive.c b/tardis/montecarlo/src/randomkit/rk_primitive.c deleted file mode 100644 index 25e70c6c7f5..00000000000 --- a/tardis/montecarlo/src/randomkit/rk_primitive.c +++ /dev/null @@ -1,520 +0,0 @@ -/* Random kit 1.6 */ -/* Primitivity test for binary polynomials of low degree */ - -/* - * Copyright (c) 2005-2006, Jean-Sebastien Roy (js@jeannot.org) - * - * Methodology inspired by scott duplichan's ppsearch code. - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -static char const rcsid[] = - "@(#) $Jeannot: rk_primitive.c,v 1.7 2006/02/19 13:48:34 js Exp $"; - -#include -#include -#include "rk_primitive.h" - -#ifndef LONG_BIT -#if ULONG_MAX <= 0xffffffffUL -#define LONG_BIT 32 -#else -#define LONG_BIT 64 -#endif -#endif - -static unsigned long modmul(unsigned long poly1, unsigned long poly2, - unsigned long modulo, unsigned long mask); -static unsigned long modpow(unsigned long polynomial, unsigned long power, - unsigned long modulo, int degree); - -/* - * For all powers i of two up to 64, list all the number of the form - * (2^i-1)/p for p a prime factor of 2^i-1. - */ -static const unsigned long divisors[][12]={ -/* 2^0-1 */ - {1UL, - 0UL}, -/* 2^1-1 */ - {1UL, - 0UL}, -/* 2^2-1 */ - {1UL, - 0UL}, -/* 2^3-1 */ - {1UL, - 0UL}, -/* 2^4-1 */ - {5UL, - 3UL, - 0UL}, -/* 2^5-1 */ - {1UL, - 0UL}, -/* 2^6-1 */ - {21UL, - 9UL, - 0UL}, -/* 2^7-1 */ - {1UL, - 0UL}, -/* 2^8-1 */ - {85UL, - 51UL, - 15UL, - 0UL}, -/* 2^9-1 */ - {73UL, - 7UL, - 0UL}, -/* 2^10-1 */ - {341UL, - 93UL, - 33UL, - 0UL}, -/* 2^11-1 */ - {89UL, - 23UL, - 0UL}, -/* 2^12-1 */ - {1365UL, - 819UL, - 585UL, - 315UL, - 0UL}, -/* 2^13-1 */ - {1UL, - 0UL}, -/* 2^14-1 */ - {5461UL, - 381UL, - 129UL, - 0UL}, -/* 2^15-1 */ - {4681UL, - 1057UL, - 217UL, - 0UL}, -/* 2^16-1 */ - {21845UL, - 13107UL, - 3855UL, - 255UL, - 0UL}, -/* 2^17-1 */ - {1UL, - 0UL}, -/* 2^18-1 */ - {87381UL, - 37449UL, - 13797UL, - 3591UL, - 0UL}, -/* 2^19-1 */ - {1UL, - 0UL}, -/* 2^20-1 */ - {349525UL, - 209715UL, - 95325UL, - 33825UL, - 25575UL, - 0UL}, -/* 2^21-1 */ - {299593UL, - 16513UL, - 6223UL, - 0UL}, -/* 2^22-1 */ - {1398101UL, - 182361UL, - 47127UL, - 6141UL, - 0UL}, -/* 2^23-1 */ - {178481UL, - 47UL, - 0UL}, -/* 2^24-1 */ - {5592405UL, - 3355443UL, - 2396745UL, - 1290555UL, - 986895UL, - 69615UL, - 0UL}, -/* 2^25-1 */ - {1082401UL, - 55831UL, - 18631UL, - 0UL}, -/* 2^26-1 */ - {22369621UL, - 24573UL, - 8193UL, - 0UL}, -/* 2^27-1 */ - {19173961UL, - 1838599UL, - 511UL, - 0UL}, -/* 2^28-1 */ - {89478485UL, - 53687091UL, - 9256395UL, - 6242685UL, - 2375535UL, - 2113665UL, - 0UL}, -/* 2^29-1 */ - {2304167UL, - 486737UL, - 256999UL, - 0UL}, -/* 2^30-1 */ - {357913941UL, - 153391689UL, - 97612893UL, - 34636833UL, - 7110873UL, - 3243933UL, - 0UL}, -/* 2^31-1 */ - {1UL, - 0UL}, -/* 2^32-1 */ - {1431655765UL, - 858993459UL, - 252645135UL, - 16711935UL, - 65535UL, - 0UL}, -#if LONG_BIT > 32 -/* 2^33-1 */ - {1227133513UL, - 373475417UL, - 96516119UL, - 14329UL, - 0UL}, -/* 2^34-1 */ - {5726623061UL, - 393213UL, - 131073UL, - 0UL}, -/* 2^35-1 */ - {1108378657UL, - 483939977UL, - 270549121UL, - 279527UL, - 0UL}, -/* 2^36-1 */ - {22906492245UL, - 13743895347UL, - 9817068105UL, - 5286113595UL, - 3616814565UL, - 1857283155UL, - 941362695UL, - 630453915UL, - 0UL}, -/* 2^37-1 */ - {616318177UL, - 223UL, - 0UL}, -/* 2^38-1 */ - {91625968981UL, - 1572861UL, - 524289UL, - 0UL}, -/* 2^39-1 */ - {78536544841UL, - 6958934353UL, - 67117057UL, - 4529623UL, - 0UL}, -/* 2^40-1 */ - {366503875925UL, - 219902325555UL, - 99955602525UL, - 64677154575UL, - 35468117025UL, - 26817356775UL, - 17825775UL, - 0UL}, -/* 2^41-1 */ - {164511353UL, - 13367UL, - 0UL}, -/* 2^42-1 */ - {1466015503701UL, - 628292358729UL, - 102280151421UL, - 34630287489UL, - 13050583119UL, - 811597437UL, - 0UL}, -/* 2^43-1 */ - {20408568497UL, - 905040953UL, - 4188889UL, - 0UL}, -/* 2^44-1 */ - {5864062014805UL, - 3518437208883UL, - 764877654105UL, - 197665011735UL, - 44312811195UL, - 25757227005UL, - 8325691455UL, - 0UL}, -/* 2^45-1 */ - {5026338869833UL, - 1134979744801UL, - 481977699847UL, - 233009086681UL, - 55759702201UL, - 1509346321UL, - 0UL}, -/* 2^46-1 */ - {23456248059221UL, - 1497207322929UL, - 394264623UL, - 25165821UL, - 0UL}, -/* 2^47-1 */ - {59862819377UL, - 31184907679UL, - 10610063UL, - 0UL}, -/* 2^48-1 */ - {93824992236885UL, - 56294995342131UL, - 40210710958665UL, - 21651921285435UL, - 16557351571215UL, - 2901803883615UL, - 1167945961455UL, - 1095233372415UL, - 418239192735UL, - 0UL}, -/* 2^49-1 */ - {4432676798593UL, - 127UL, - 0UL}, -/* 2^50-1 */ - {375299968947541UL, - 102354536985693UL, - 36319351833633UL, - 4485656999373UL, - 1873377548823UL, - 625152641223UL, - 277931351973UL, - 0UL}, -/* 2^51-1 */ - {321685687669321UL, - 21862134113449UL, - 1050769861729UL, - 202518195313UL, - 17180000257UL, - 0UL}, -/* 2^52-1 */ - {1501199875790165UL, - 900719925474099UL, - 84973577874915UL, - 28685347945035UL, - 2792064245115UL, - 1649066139645UL, - 549822930945UL, - 0UL}, -/* 2^53-1 */ - {1416003655831UL, - 129728784761UL, - 441650591UL, - 0UL}, -/* 2^54-1 */ - {6004799503160661UL, - 2573485501354569UL, - 948126237341157UL, - 246772582321671UL, - 206561081853UL, - 68585259519UL, - 0UL}, -/* 2^55-1 */ - {1566469435607129UL, - 1162219258676257UL, - 404817944033303UL, - 40895342813807UL, - 11290754314937UL, - 178394823847UL, - 0UL}, -/* 2^56-1 */ - {24019198012642645UL, - 14411518807585587UL, - 4238682002231055UL, - 2484744621997515UL, - 1675758000882045UL, - 637677823344495UL, - 567382630219905UL, - 4563402735UL, - 0UL}, -/* 2^57-1 */ - {20587884010836553UL, - 4451159405623UL, - 274878431233UL, - 118823881393UL, - 0UL}, -/* 2^58-1 */ - {96076792050570581UL, - 4885260612740877UL, - 1237040240994471UL, - 261314937580881UL, - 137975287770087UL, - 95026151247UL, - 0UL}, -/* 2^59-1 */ - {3203431780337UL, - 179951UL, - 0UL}, -/* 2^60-1 */ - {384307168202282325UL, - 230584300921369395UL, - 164703072086692425UL, - 104811045873349725UL, - 88686269585142075UL, - 37191016277640225UL, - 28120036697727975UL, - 18900352534538475UL, - 7635241752363225UL, - 3483146539597725UL, - 872764197279975UL, - 0UL}, -/* 2^61-1 */ - {1UL, - 0UL}, -/* 2^62-1 */ - {1537228672809129301UL, - 6442450941UL, - 2147483649UL, - 0UL}, -/* 2^63-1 */ - {1317624576693539401UL, - 126347562148695559UL, - 72624976668147841UL, - 27369056489183311UL, - 99457304386111UL, - 14197294936951UL, - 0UL}, -/* 2^64-1 */ - {6148914691236517205UL, - 3689348814741910323UL, - 1085102592571150095UL, - 71777214294589695UL, - 28778071877862015UL, - 281470681808895UL, - 2753074036095UL, - 0UL}, -#if LONG_BIT > 64 -#error Factorization of numbers up to 2^LONG_BIT required -#endif -#endif -}; - -/* - * Modular multiply for two binary polynomial - * mask is 1UL << the degree of the modulus. - */ -unsigned long modmul(unsigned long poly1, unsigned long poly2, - unsigned long modulo, unsigned long mask) -{ - unsigned long result = 0; - - for (; poly1; poly1 >>= 1) - { - if (poly1 & 1) - result ^= poly2; - - poly2 <<= 1; - if (poly2 & mask) - poly2 ^= modulo; - } - return result; -} - -/* - * Modular exponentiation for a binary polynomial - * degree is the degree of the modulus. - */ -unsigned long modpow(unsigned long polynomial, unsigned long power, - unsigned long modulo, int degree) -{ - unsigned long result = 1, mask = 1UL << degree; - - for (; power; power >>= 1) - { - if (power & 1) - result = modmul(result, polynomial, modulo, mask); - polynomial = modmul(polynomial, polynomial, modulo, mask); - } - return result; -} - -/* - * Test the primitivity of a polynomial - */ -int rk_isprimitive(unsigned long polynomial) -{ - unsigned long pelement = 2, temp = polynomial >> 1; - int k, degree = 0, weight = 1; - - /* Special case for polynomials of degree < 2 */ - if (polynomial < 4) - return (polynomial == 3) || (polynomial == 1); - - /* A binary primitive polynomial has a constant term */ - if (!(polynomial & 1)) - return 0; - - /* - * A binary primitive polynomial of degree > 1 has an odd number of terms. - * temp ^= temp >> 16; temp ^= temp >> 8; ... would be sligthly faster. - * Compute the degree at the same time. - */ - for (; temp; degree++, temp >>= 1) - weight += temp & 1; - if (!(weight & 1)) - return 0; - - /* - * Check if the period is 2^degree-1. - * Sufficient if 2^degree-1 is prime. - */ - if (modpow(pelement, 1UL << degree, polynomial, degree) != pelement) - return 0; - - if (divisors[degree][0] != 1) - /* Primitivity test */ - for (k = 0; divisors[degree][k]; k++) - if (modpow(pelement, divisors[degree][k], polynomial, degree) == 1) - return 0; - - return 1; -} diff --git a/tardis/montecarlo/src/randomkit/rk_primitive.h b/tardis/montecarlo/src/randomkit/rk_primitive.h deleted file mode 100644 index f7ab374cef5..00000000000 --- a/tardis/montecarlo/src/randomkit/rk_primitive.h +++ /dev/null @@ -1,54 +0,0 @@ -/* Random kit 1.6 */ -/* Primitivity test for binary polynomials of low degree */ - -/* - * Copyright (c) 2005-2006, Jean-Sebastien Roy (js@jeannot.org) - * - * Methodology inspired by scott duplichan's ppsearch code. - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -/* @(#) $Jeannot: rk_primitive.h,v 1.6 2006/02/19 13:48:34 js Exp $ */ - -#ifndef _RK_PRIMITIVE_ -#define _RK_PRIMITIVE_ - -#ifdef __cplusplus -extern "C" { -#endif - -/* - * Return 1 if the binary polynomial is primitive. - * - * Note that if p is primitive, the the polynomial obtained by reversing the - * bits of p is also primitive. (see list_primitive.c for an example) - * - * Typical use: - * int test; - * test = rk_isprimitive(3, &divisors); - */ -extern int rk_isprimitive(unsigned long polynomial); - -#ifdef __cplusplus -} -#endif - -#endif /* _RK_PRIMITIVE_ */ diff --git a/tardis/montecarlo/src/randomkit/rk_sobol.c b/tardis/montecarlo/src/randomkit/rk_sobol.c deleted file mode 100644 index c6a016bb033..00000000000 --- a/tardis/montecarlo/src/randomkit/rk_sobol.c +++ /dev/null @@ -1,991 +0,0 @@ -/* Random kit 1.6 */ - -/* - * Copyright (c) 2004-2006, Jean-Sebastien Roy (js@jeannot.org) - * - * Original algorithm from Numerical Recipes, 2nd edition, by Press et al. - * The inverse normal cdf formulas are from Peter J. Acklam. - * The initialization directions were found in Ferdinando Ametrano's - * implementation in QuantLib. - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -static char const rcsid[] = - "@(#) $Jeannot: rk_sobol.c,v 1.9 2006/02/19 13:48:34 js Exp $"; - -#include -#include -#include -#include "rk_sobol.h" -#include "rk_mt.h" -#include "rk_primitive.h" - -#ifndef M_PI -#define M_PI 3.14159265358979323846 -#endif -#ifndef M_SQRT1_2 -#define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */ -#endif -#define RK_SOBOL_M_SQRT2PI 2.506628274631000502415 /* sqrt(2*pi) */ - -#ifndef LONG_BIT -#if ULONG_MAX <= 0xffffffffUL -#define LONG_BIT 32 -#else -#define LONG_BIT 64 -#endif -#endif - -char *rk_sobol_strerror[] = -{ - "no error", - "invalid dimension", - "too many numbers generated", - "not enough memory" -}; - -static double inverse_normal(double p); - -/* - * Sobol/Levitan coefficients of the free direction integers as given - * by Bratley, P., Fox, B.L. (1988) - */ - -const unsigned long rk_sobol_SLdirections[] = { - 1, - 1, 1, - 1, 3, 7, - 1, 1, 5, - 1, 3, 1, 1, - 1, 1, 3, 7, - 1, 3, 3, 9, 9, - 1, 3, 7, 13, 3, - 1, 1, 5, 11, 27, - 1, 3, 5, 1, 15, - 1, 1, 7, 3, 29, - 1, 3, 7, 7, 21, - 1, 1, 1, 9, 23, 37, - 1, 3, 3, 5, 19, 33, - 1, 1, 3, 13, 11, 7, - 1, 1, 7, 13, 25, 5, - 1, 3, 5, 11, 7, 11, - 1, 1, 1, 3, 13, 39, - 1, 3, 1, 15, 17, 63, 13, - 1, 1, 5, 5, 1, 27, 33, - 1, 3, 3, 3, 25, 17, 115, - 1, 1, 3, 15, 29, 15, 41, - 1, 3, 1, 7, 3, 23, 79, - 1, 3, 7, 9, 31, 29, 17, - 1, 1, 5, 13, 11, 3, 29, - 1, 3, 1, 9, 5, 21, 119, - 1, 1, 3, 1, 23, 13, 75, - 1, 3, 3, 11, 27, 31, 73, - 1, 1, 7, 7, 19, 25, 105, - 1, 3, 5, 5, 21, 9, 7, - 1, 1, 1, 15, 5, 49, 59, - 1, 1, 1, 1, 1, 33, 65, - 1, 3, 5, 15, 17, 19, 21, - 1, 1, 7, 11, 13, 29, 3, - 1, 3, 7, 5, 7, 11, 113, - 1, 1, 5, 3, 15, 19, 61, - 1, 3, 1, 1, 9, 27, 89, 7, - 1, 1, 3, 7, 31, 15, 45, 23, - 1, 3, 3, 9, 9, 25, 107, 39, - 0 -}; - -/* - * Lemieux coefficients of the free direction integers as given - * in QuantLib by Christiane Lemieux, private communication, September 2004 - */ - -const unsigned long rk_sobol_Ldirections[] = { - 1, - 1, 1, - 1, 3, 7, - 1, 1, 5, - 1, 3, 1, 1, - 1, 1, 3, 7, - 1, 3, 3, 9, 9, - 1, 3, 7, 13, 3, - 1, 1, 5, 11, 27, - 1, 3, 5, 1, 15, - 1, 1, 7, 3, 29, - 1, 3, 7, 7, 21, - 1, 1, 1, 9, 23, 37, - 1, 3, 3, 5, 19, 33, - 1, 1, 3, 13, 11, 7, - 1, 1, 7, 13, 25, 5, - 1, 3, 5, 11, 7, 11, - 1, 1, 1, 3, 13, 39, - 1, 3, 1, 15, 17, 63, 13, - 1, 1, 5, 5, 1, 27, 33, - 1, 3, 3, 3, 25, 17, 115, - 1, 1, 3, 15, 29, 15, 41, - 1, 3, 1, 7, 3, 23, 79, - 1, 3, 7, 9, 31, 29, 17, - 1, 1, 5, 13, 11, 3, 29, - 1, 3, 1, 9, 5, 21, 119, - 1, 1, 3, 1, 23, 13, 75, - 1, 3, 3, 11, 27, 31, 73, - 1, 1, 7, 7, 19, 25, 105, - 1, 3, 5, 5, 21, 9, 7, - 1, 1, 1, 15, 5, 49, 59, - 1, 1, 1, 1, 1, 33, 65, - 1, 3, 5, 15, 17, 19, 21, - 1, 1, 7, 11, 13, 29, 3, - 1, 3, 7, 5, 7, 11, 113, - 1, 1, 5, 3, 15, 19, 61, - 1, 3, 1, 1, 9, 27, 89, 7, - 1, 1, 3, 7, 31, 15, 45, 23, - 1, 3, 3, 9, 9, 25, 107, 39, - 1, 1, 3, 13, 7, 35, 61, 91, - 1, 1, 7, 11, 5, 35, 55, 75, - 1, 3, 5, 5, 11, 23, 29, 139, - 1, 1, 1, 7, 11, 15, 17, 81, - 1, 1, 7, 9, 5, 57, 79, 103, - 1, 1, 7, 13, 19, 5, 5, 185, - 1, 3, 1, 3, 13, 57, 97, 131, - 1, 1, 5, 5, 21, 25, 125, 197, - 1, 3, 3, 9, 31, 11, 103, 201, - 1, 1, 5, 3, 7, 25, 51, 121, - 1, 3, 7, 15, 19, 53, 73, 189, - 1, 1, 1, 15, 19, 55, 27, 183, - 1, 1, 7, 13, 3, 29, 109, 69, - 1, 1, 5, 15, 15, 23, 15, 1, 57, - 1, 3, 1, 3, 23, 55, 43, 143, 397, - 1, 1, 3, 11, 29, 9, 35, 131, 411, - 1, 3, 1, 7, 27, 39, 103, 199, 277, - 1, 3, 7, 3, 19, 55, 127, 67, 449, - 1, 3, 7, 3, 5, 29, 45, 85, 3, - 1, 3, 5, 5, 13, 23, 75, 245, 453, - 1, 3, 1, 15, 21, 47, 3, 77, 165, - 1, 1, 7, 9, 15, 5, 117, 73, 473, - 1, 3, 1, 9, 1, 21, 13, 173, 313, - 1, 1, 7, 3, 11, 45, 63, 77, 49, - 1, 1, 1, 1, 1, 25, 123, 39, 259, - 1, 1, 1, 5, 23, 11, 59, 11, 203, - 1, 3, 3, 15, 21, 1, 73, 71, 421, - 1, 1, 5, 11, 15, 31, 115, 95, 217, - 1, 1, 3, 3, 7, 53, 37, 43, 439, - 1, 1, 1, 1, 27, 53, 69, 159, 321, - 1, 1, 5, 15, 29, 17, 19, 43, 449, - 1, 1, 3, 9, 1, 55, 121, 205, 255, - 1, 1, 3, 11, 9, 47, 107, 11, 417, - 1, 1, 1, 5, 17, 25, 21, 83, 95, - 1, 3, 5, 13, 31, 25, 61, 157, 407, - 1, 1, 7, 9, 25, 33, 41, 35, 17, - 1, 3, 7, 15, 13, 39, 61, 187, 461, - 1, 3, 7, 13, 5, 57, 23, 177, 435, - 1, 1, 3, 15, 11, 27, 115, 5, 337, - 1, 3, 7, 3, 15, 63, 61, 171, 339, - 1, 3, 3, 13, 15, 61, 59, 47, 1, - 1, 1, 5, 15, 13, 5, 39, 83, 329, - 1, 1, 5, 5, 5, 27, 25, 39, 301, - 1, 1, 5, 11, 31, 41, 35, 233, 27, - 1, 3, 5, 15, 7, 37, 119, 171, 419, - 1, 3, 5, 5, 3, 29, 21, 189, 417, - 1, 1, 1, 1, 21, 41, 117, 119, 351, - 1, 1, 3, 1, 7, 27, 87, 19, 213, - 1, 1, 1, 1, 17, 7, 97, 217, 477, - 1, 1, 7, 1, 29, 61, 103, 231, 269, - 1, 1, 7, 13, 9, 27, 107, 207, 311, - 1, 1, 7, 5, 25, 21, 107, 179, 423, - 1, 3, 5, 11, 7, 1, 17, 245, 281, - 1, 3, 5, 9, 1, 5, 53, 59, 125, - 1, 1, 7, 1, 31, 57, 71, 245, 125, - 1, 1, 7, 5, 5, 57, 53, 253, 441, - 1, 3, 1, 13, 19, 35, 119, 235, 381, - 1, 3, 1, 7, 19, 59, 115, 33, 361, - 1, 1, 3, 5, 13, 1, 49, 143, 501, - 1, 1, 3, 5, 1, 63, 101, 85, 189, - 1, 1, 5, 11, 27, 63, 13, 131, 5, - 1, 1, 5, 7, 15, 45, 75, 59, 455, 585, - 1, 3, 1, 3, 7, 7, 111, 23, 119, 959, - 1, 3, 3, 9, 11, 41, 109, 163, 161, 879, - 1, 3, 5, 1, 21, 41, 121, 183, 315, 219, - 1, 1, 3, 9, 15, 3, 9, 223, 441, 929, - 1, 1, 7, 9, 3, 5, 93, 57, 253, 457, - 1, 1, 7, 13, 15, 29, 83, 21, 35, 45, - 1, 1, 3, 7, 13, 61, 119, 219, 85, 505, - 1, 1, 3, 3, 17, 13, 35, 197, 291, 109, - 1, 1, 3, 3, 5, 1, 113, 103, 217, 253, - 1, 1, 7, 1, 15, 39, 63, 223, 17, 9, - 1, 3, 7, 1, 17, 29, 67, 103, 495, 383, - 1, 3, 3, 15, 31, 59, 75, 165, 51, 913, - 1, 3, 7, 9, 5, 27, 79, 219, 233, 37, - 1, 3, 5, 15, 1, 11, 15, 211, 417, 811, - 1, 3, 5, 3, 29, 27, 39, 137, 407, 231, - 1, 1, 3, 5, 29, 43, 125, 135, 109, 67, - 1, 1, 1, 5, 11, 39, 107, 159, 323, 381, - 1, 1, 1, 1, 9, 11, 33, 55, 169, 253, - 1, 3, 5, 5, 11, 53, 63, 101, 251, 897, - 1, 3, 7, 1, 25, 15, 83, 119, 53, 157, - 1, 3, 5, 13, 5, 5, 3, 195, 111, 451, - 1, 3, 1, 15, 11, 1, 19, 11, 307, 777, - 1, 3, 7, 11, 5, 5, 17, 231, 345, 981, - 1, 1, 3, 3, 1, 33, 83, 201, 57, 475, - 1, 3, 7, 7, 17, 13, 35, 175, 499, 809, - 1, 1, 5, 3, 3, 17, 103, 119, 499, 865, - 1, 1, 1, 11, 27, 25, 37, 121, 401, 11, - 1, 1, 1, 11, 9, 25, 25, 241, 403, 3, - 1, 1, 1, 1, 11, 1, 39, 163, 231, 573, - 1, 1, 1, 13, 13, 21, 75, 185, 99, 545, - 1, 1, 1, 15, 3, 63, 69, 11, 173, 315, - 1, 3, 5, 15, 11, 3, 95, 49, 123, 765, - 1, 1, 1, 15, 3, 63, 77, 31, 425, 711, - 1, 1, 7, 15, 1, 37, 119, 145, 489, 583, - 1, 3, 5, 15, 3, 49, 117, 211, 165, 323, - 1, 3, 7, 1, 27, 63, 77, 201, 225, 803, - 1, 1, 1, 11, 23, 35, 67, 21, 469, 357, - 1, 1, 7, 7, 9, 7, 25, 237, 237, 571, - 1, 1, 3, 15, 29, 5, 107, 109, 241, 47, - 1, 3, 5, 11, 27, 63, 29, 13, 203, 675, - 1, 1, 3, 9, 9, 11, 103, 179, 449, 263, - 1, 3, 5, 11, 29, 63, 53, 151, 259, 223, - 1, 1, 3, 7, 9, 25, 5, 197, 237, 163, - 1, 3, 7, 13, 5, 57, 67, 193, 147, 241, - 1, 1, 5, 15, 15, 33, 17, 67, 161, 341, - 1, 1, 3, 13, 17, 43, 21, 197, 441, 985, - 1, 3, 1, 5, 15, 33, 33, 193, 305, 829, - 1, 1, 1, 13, 19, 27, 71, 187, 477, 239, - 1, 1, 1, 9, 9, 17, 41, 177, 229, 983, - 1, 3, 5, 9, 15, 45, 97, 205, 43, 767, - 1, 1, 1, 9, 31, 31, 77, 159, 395, 809, - 1, 3, 3, 3, 29, 19, 73, 123, 165, 307, - 1, 3, 1, 7, 5, 11, 77, 227, 355, 403, - 1, 3, 5, 5, 25, 31, 1, 215, 451, 195, - 1, 3, 7, 15, 29, 37, 101, 241, 17, 633, - 1, 1, 5, 1, 11, 3, 107, 137, 489, 5, - 1, 1, 1, 7, 19, 19, 75, 85, 471, 355, - 1, 1, 3, 3, 9, 13, 113, 167, 13, 27, - 1, 3, 5, 11, 21, 3, 89, 205, 377, 307, - 1, 1, 1, 9, 31, 61, 65, 9, 391, 141, 867, - 1, 1, 1, 9, 19, 19, 61, 227, 241, 55, 161, - 1, 1, 1, 11, 1, 19, 7, 233, 463, 171, 1941, - 1, 1, 5, 7, 25, 13, 103, 75, 19, 1021, 1063, - 1, 1, 1, 15, 17, 17, 79, 63, 391, 403, 1221, - 1, 3, 3, 11, 29, 25, 29, 107, 335, 475, 963, - 1, 3, 5, 1, 31, 33, 49, 43, 155, 9, 1285, - 1, 1, 5, 5, 15, 47, 39, 161, 357, 863, 1039, - 1, 3, 7, 15, 1, 39, 47, 109, 427, 393, 1103, - 1, 1, 1, 9, 9, 29, 121, 233, 157, 99, 701, - 1, 1, 1, 7, 1, 29, 75, 121, 439, 109, 993, - 1, 1, 1, 9, 5, 1, 39, 59, 89, 157, 1865, - 1, 1, 5, 1, 3, 37, 89, 93, 143, 533, 175, - 1, 1, 3, 5, 7, 33, 35, 173, 159, 135, 241, - 1, 1, 1, 15, 17, 37, 79, 131, 43, 891, 229, - 1, 1, 1, 1, 1, 35, 121, 177, 397, 1017, 583, - 1, 1, 3, 15, 31, 21, 43, 67, 467, 923, 1473, - 1, 1, 1, 7, 1, 33, 77, 111, 125, 771, 1975, - 1, 3, 7, 13, 1, 51, 113, 139, 245, 573, 503, - 1, 3, 1, 9, 21, 49, 15, 157, 49, 483, 291, - 1, 1, 1, 1, 29, 35, 17, 65, 403, 485, 1603, - 1, 1, 1, 7, 19, 1, 37, 129, 203, 321, 1809, - 1, 3, 7, 15, 15, 9, 5, 77, 29, 485, 581, - 1, 1, 3, 5, 15, 49, 97, 105, 309, 875, 1581, - 1, 3, 5, 1, 5, 19, 63, 35, 165, 399, 1489, - 1, 3, 5, 3, 23, 5, 79, 137, 115, 599, 1127, - 1, 1, 7, 5, 3, 61, 27, 177, 257, 91, 841, - 1, 1, 3, 5, 9, 31, 91, 209, 409, 661, 159, - 1, 3, 1, 15, 23, 39, 23, 195, 245, 203, 947, - 1, 1, 3, 1, 15, 59, 67, 95, 155, 461, 147, - 1, 3, 7, 5, 23, 25, 87, 11, 51, 449, 1631, - 1, 1, 1, 1, 17, 57, 7, 197, 409, 609, 135, - 1, 1, 1, 9, 1, 61, 115, 113, 495, 895, 1595, - 1, 3, 7, 15, 9, 47, 121, 211, 379, 985, 1755, - 1, 3, 1, 3, 7, 57, 27, 231, 339, 325, 1023, - 1, 1, 1, 1, 19, 63, 63, 239, 31, 643, 373, - 1, 3, 1, 11, 19, 9, 7, 171, 21, 691, 215, - 1, 1, 5, 13, 11, 57, 39, 211, 241, 893, 555, - 1, 1, 7, 5, 29, 21, 45, 59, 509, 223, 491, - 1, 1, 7, 9, 15, 61, 97, 75, 127, 779, 839, - 1, 1, 7, 15, 17, 33, 75, 237, 191, 925, 681, - 1, 3, 5, 7, 27, 57, 123, 111, 101, 371, 1129, - 1, 3, 5, 5, 29, 45, 59, 127, 229, 967, 2027, - 1, 1, 1, 1, 17, 7, 23, 199, 241, 455, 135, - 1, 1, 7, 15, 27, 29, 105, 171, 337, 503, 1817, - 1, 1, 3, 7, 21, 35, 61, 71, 405, 647, 2045, - 1, 1, 1, 1, 1, 15, 65, 167, 501, 79, 737, - 1, 1, 5, 1, 3, 49, 27, 189, 341, 615, 1287, - 1, 1, 1, 9, 1, 7, 31, 159, 503, 327, 1613, - 1, 3, 3, 3, 3, 23, 99, 115, 323, 997, 987, - 1, 1, 1, 9, 19, 33, 93, 247, 509, 453, 891, - 1, 1, 3, 1, 13, 19, 35, 153, 161, 633, 445, - 1, 3, 5, 15, 31, 5, 87, 197, 183, 783, 1823, - 1, 1, 7, 5, 19, 63, 69, 221, 129, 231, 1195, - 1, 1, 5, 5, 13, 23, 19, 231, 245, 917, 379, - 1, 3, 1, 15, 19, 43, 27, 223, 171, 413, 125, - 1, 1, 1, 9, 1, 59, 21, 15, 509, 207, 589, - 1, 3, 5, 3, 19, 31, 113, 19, 23, 733, 499, - 1, 1, 7, 1, 19, 51, 101, 165, 47, 925, 1093, - 1, 3, 3, 9, 15, 21, 43, 243, 237, 461, 1361, - 1, 1, 1, 9, 17, 15, 75, 75, 113, 715, 1419, - 1, 1, 7, 13, 17, 1, 99, 15, 347, 721, 1405, - 1, 1, 7, 15, 7, 27, 23, 183, 39, 59, 571, - 1, 3, 5, 9, 7, 43, 35, 165, 463, 567, 859, - 1, 3, 3, 11, 15, 19, 17, 129, 311, 343, 15, - 1, 1, 1, 15, 31, 59, 63, 39, 347, 359, 105, - 1, 1, 1, 15, 5, 43, 87, 241, 109, 61, 685, - 1, 1, 7, 7, 9, 39, 121, 127, 369, 579, 853, - 1, 1, 1, 1, 17, 15, 15, 95, 325, 627, 299, - 1, 1, 3, 13, 31, 53, 85, 111, 289, 811, 1635, - 1, 3, 7, 1, 19, 29, 75, 185, 153, 573, 653, - 1, 3, 7, 1, 29, 31, 55, 91, 249, 247, 1015, - 1, 3, 5, 7, 1, 49, 113, 139, 257, 127, 307, - 1, 3, 5, 9, 15, 15, 123, 105, 105, 225, 1893, - 1, 3, 3, 1, 15, 5, 105, 249, 73, 709, 1557, - 1, 1, 1, 9, 17, 31, 113, 73, 65, 701, 1439, - 1, 3, 5, 15, 13, 21, 117, 131, 243, 859, 323, - 1, 1, 1, 9, 19, 15, 69, 149, 89, 681, 515, - 1, 1, 1, 5, 29, 13, 21, 97, 301, 27, 967, - 1, 1, 3, 3, 15, 45, 107, 227, 495, 769, 1935, - 1, 1, 1, 11, 5, 27, 41, 173, 261, 703, 1349, - 1, 3, 3, 3, 11, 35, 97, 43, 501, 563, 1331, - 1, 1, 1, 7, 1, 17, 87, 17, 429, 245, 1941, - 1, 1, 7, 15, 29, 13, 1, 175, 425, 233, 797, - 1, 1, 3, 11, 21, 57, 49, 49, 163, 685, 701, - 1, 3, 3, 7, 11, 45, 107, 111, 379, 703, 1403, - 1, 1, 7, 3, 21, 7, 117, 49, 469, 37, 775, - 1, 1, 5, 15, 31, 63, 101, 77, 507, 489, 1955, - 1, 3, 3, 11, 19, 21, 101, 255, 203, 673, 665, - 1, 3, 3, 15, 17, 47, 125, 187, 271, 899, 2003, - 1, 1, 7, 7, 1, 35, 13, 235, 5, 337, 905, - 1, 3, 1, 15, 1, 43, 1, 27, 37, 695, 1429, - 1, 3, 1, 11, 21, 27, 93, 161, 299, 665, 495, - 1, 3, 3, 15, 3, 1, 81, 111, 105, 547, 897, - 1, 3, 5, 1, 3, 53, 97, 253, 401, 827, 1467, - 1, 1, 1, 5, 19, 59, 105, 125, 271, 351, 719, - 1, 3, 5, 13, 7, 11, 91, 41, 441, 759, 1827, - 1, 3, 7, 11, 29, 61, 61, 23, 307, 863, 363, - 1, 1, 7, 1, 15, 35, 29, 133, 415, 473, 1737, - 1, 1, 1, 13, 7, 33, 35, 225, 117, 681, 1545, - 1, 1, 1, 3, 5, 41, 83, 247, 13, 373, 1091, - 1, 3, 1, 13, 25, 61, 71, 217, 233, 313, 547, - 1, 3, 1, 7, 3, 29, 3, 49, 93, 465, 15, - 1, 1, 1, 9, 17, 61, 99, 163, 129, 485, 1087, - 1, 1, 1, 9, 9, 33, 31, 163, 145, 649, 253, - 1, 1, 1, 1, 17, 63, 43, 235, 287, 111, 567, - 1, 3, 5, 13, 29, 7, 11, 69, 153, 127, 449, - 1, 1, 5, 9, 11, 21, 15, 189, 431, 493, 1219, - 1, 1, 1, 15, 19, 5, 47, 91, 399, 293, 1743, - 1, 3, 3, 11, 29, 53, 53, 225, 409, 303, 333, - 1, 1, 1, 15, 31, 31, 21, 81, 147, 287, 1753, - 1, 3, 5, 5, 5, 63, 35, 125, 41, 687, 1793, - 1, 1, 1, 9, 19, 59, 107, 219, 455, 971, 297, - 1, 1, 3, 5, 3, 51, 121, 31, 245, 105, 1311, - 1, 3, 1, 5, 5, 57, 75, 107, 161, 431, 1693, - 1, 3, 1, 3, 19, 53, 27, 31, 191, 565, 1015, - 1, 3, 5, 13, 9, 41, 35, 249, 287, 49, 123, - 1, 1, 5, 7, 27, 17, 21, 3, 151, 885, 1165, - 1, 1, 7, 1, 15, 17, 65, 139, 427, 339, 1171, - 1, 1, 1, 5, 23, 5, 9, 89, 321, 907, 391, - 1, 1, 7, 9, 15, 1, 77, 71, 87, 701, 917, - 1, 1, 7, 1, 17, 37, 115, 127, 469, 779, 1543, - 1, 3, 7, 3, 5, 61, 15, 37, 301, 951, 1437, - 1, 1, 1, 13, 9, 51, 127, 145, 229, 55, 1567, - 1, 3, 7, 15, 19, 47, 53, 153, 295, 47, 1337, - 1, 3, 3, 5, 11, 31, 29, 133, 327, 287, 507, - 1, 1, 7, 7, 25, 31, 37, 199, 25, 927, 1317, - 1, 1, 7, 9, 3, 39, 127, 167, 345, 467, 759, - 1, 1, 1, 1, 31, 21, 15, 101, 293, 787, 1025, - 1, 1, 5, 3, 11, 41, 105, 109, 149, 837, 1813, - 1, 1, 3, 5, 29, 13, 19, 97, 309, 901, 753, - 1, 1, 7, 1, 19, 17, 31, 39, 173, 361, 1177, - 1, 3, 3, 3, 3, 41, 81, 7, 341, 491, 43, - 1, 1, 7, 7, 31, 35, 29, 77, 11, 335, 1275, - 1, 3, 3, 15, 17, 45, 19, 63, 151, 849, 129, - 1, 1, 7, 5, 7, 13, 47, 73, 79, 31, 499, - 1, 3, 1, 11, 1, 41, 59, 151, 247, 115, 1295, - 1, 1, 1, 9, 31, 37, 73, 23, 295, 483, 179, - 1, 3, 1, 15, 13, 63, 81, 27, 169, 825, 2037, - 1, 3, 5, 15, 7, 11, 73, 1, 451, 101, 2039, - 1, 3, 5, 3, 13, 53, 31, 137, 173, 319, 1521, - 1, 3, 1, 3, 29, 1, 73, 227, 377, 337, 1189, - 1, 3, 3, 13, 27, 9, 31, 101, 229, 165, 1983, - 1, 3, 1, 13, 13, 19, 19, 111, 319, 421, 223, - 1, 1, 7, 15, 25, 37, 61, 55, 359, 255, 1955, - 1, 1, 5, 13, 17, 43, 49, 215, 383, 915, 51, - 1, 1, 3, 1, 3, 7, 13, 119, 155, 585, 967, - 1, 3, 1, 13, 1, 63, 125, 21, 103, 287, 457, - 1, 1, 7, 1, 31, 17, 125, 137, 345, 379, 1925, - 1, 1, 3, 5, 5, 25, 119, 153, 455, 271, 2023, - 1, 1, 7, 9, 9, 37, 115, 47, 5, 255, 917, - 1, 3, 5, 3, 31, 21, 75, 203, 489, 593, 1, - 1, 3, 7, 15, 19, 63, 123, 153, 135, 977, 1875, - 1, 1, 1, 1, 5, 59, 31, 25, 127, 209, 745, - 1, 1, 1, 1, 19, 45, 67, 159, 301, 199, 535, - 1, 1, 7, 1, 31, 17, 19, 225, 369, 125, 421, - 1, 3, 3, 11, 7, 59, 115, 197, 459, 469, 1055, - 1, 3, 1, 3, 27, 45, 35, 131, 349, 101, 411, - 1, 3, 7, 11, 9, 3, 67, 145, 299, 253, 1339, - 1, 3, 3, 11, 9, 37, 123, 229, 273, 269, 515, - 1, 3, 7, 15, 11, 25, 75, 5, 367, 217, 951, - 1, 1, 3, 7, 9, 23, 63, 237, 385, 159, 1273, - 1, 1, 5, 11, 23, 5, 55, 193, 109, 865, 663, - 1, 1, 7, 15, 1, 57, 17, 141, 51, 217, 1259, - 1, 1, 3, 3, 15, 7, 89, 233, 71, 329, 203, - 1, 3, 7, 11, 11, 1, 19, 155, 89, 437, 573, - 1, 3, 1, 9, 27, 61, 47, 109, 161, 913, 1681, - 1, 1, 7, 15, 1, 33, 19, 15, 23, 913, 989, - 1, 3, 1, 1, 25, 39, 119, 193, 13, 571, 157, - 1, 1, 7, 13, 9, 55, 59, 147, 361, 935, 515, - 1, 1, 1, 9, 7, 59, 67, 117, 71, 855, 1493, - 1, 3, 1, 3, 13, 19, 57, 141, 305, 275, 1079, - 1, 1, 1, 9, 17, 61, 33, 7, 43, 931, 781, - 1, 1, 3, 1, 11, 17, 21, 97, 295, 277, 1721, - 1, 3, 1, 13, 15, 43, 11, 241, 147, 391, 1641, - 1, 1, 1, 1, 1, 19, 37, 21, 255, 263, 1571, - 1, 1, 3, 3, 23, 59, 89, 17, 475, 303, 757, 543, - 1, 3, 3, 9, 11, 55, 35, 159, 139, 203, 1531, 1825, - 1, 1, 5, 3, 17, 53, 51, 241, 269, 949, 1373, 325, - 1, 3, 7, 7, 5, 29, 91, 149, 239, 193, 1951, 2675, - 1, 3, 5, 1, 27, 33, 69, 11, 51, 371, 833, 2685, - 1, 1, 1, 15, 1, 17, 35, 57, 171, 1007, 449, 367, - 1, 1, 1, 7, 25, 61, 73, 219, 379, 53, 589, 4065, - 1, 3, 5, 13, 21, 29, 45, 19, 163, 169, 147, 597, - 1, 1, 5, 11, 21, 27, 7, 17, 237, 591, 255, 1235, - 1, 1, 7, 7, 17, 41, 69, 237, 397, 173, 1229, 2341, - 1, 1, 3, 1, 1, 33, 125, 47, 11, 783, 1323, 2469, - 1, 3, 1, 11, 3, 39, 35, 133, 153, 55, 1171, 3165, - 1, 1, 5, 11, 27, 23, 103, 245, 375, 753, 477, 2165, - 1, 3, 1, 15, 15, 49, 127, 223, 387, 771, 1719, 1465, - 1, 1, 1, 9, 11, 9, 17, 185, 239, 899, 1273, 3961, - 1, 1, 3, 13, 11, 51, 73, 81, 389, 647, 1767, 1215, - 1, 3, 5, 15, 19, 9, 69, 35, 349, 977, 1603, 1435, - 1, 1, 1, 1, 19, 59, 123, 37, 41, 961, 181, 1275, - 1, 1, 1, 1, 31, 29, 37, 71, 205, 947, 115, 3017, - 1, 1, 7, 15, 5, 37, 101, 169, 221, 245, 687, 195, - 1, 1, 1, 1, 19, 9, 125, 157, 119, 283, 1721, 743, - 1, 1, 7, 3, 1, 7, 61, 71, 119, 257, 1227, 2893, - 1, 3, 3, 3, 25, 41, 25, 225, 31, 57, 925, 2139, - 0 -}; - - -/* - * coefficients of the free direction integers as given in - * "Monte Carlo Methods in Finance", by Peter J�ckel, section 8.3 - */ - -const unsigned long rk_sobol_Jdirections[] = { - 1, - 1, 1, - 1, 3, 7, - 1, 1, 5, - 1, 3, 1, 1, - 1, 1, 3, 7, - 1, 3, 3, 9, 9, - 1, 3, 7, 7, 21, - 1, 1, 5, 11, 27, - 1, 1, 7, 3, 29, - 1, 3, 7, 13, 3, - 1, 3, 5, 1, 15, - 1, 1, 1, 9, 23, 37, - 1, 1, 3, 13, 11, 7, - 1, 3, 3, 5, 19, 33, - 1, 1, 7, 13, 25, 5, - 1, 1, 1, 3, 13, 39, - 1, 3, 5, 11, 7, 11, - 1, 3, 1, 7, 3, 23, 79, - 1, 3, 1, 15, 17, 63, 13, - 1, 3, 3, 3, 25, 17, 115, - 1, 3, 7, 9, 31, 29, 17, - 1, 1, 3, 15, 29, 15, 41, - 1, 3, 1, 9, 5, 21, 119, - 1, 1, 5, 5, 1, 27, 33, - 1, 1, 3, 1, 23, 13, 75, - 1, 1, 7, 7, 19, 25, 105, - 1, 3, 5, 5, 21, 9, 7, - 1, 1, 1, 15, 5, 49, 59, - 1, 3, 5, 15, 17, 19, 21, - 1, 1, 7, 11, 13, 29, 3, - 0 -}; - -/* - * 0 terminated list of primitive polynomials to speed up initialization - * All polynomials up to degree 13 (ie. 1111 polynomials) - */ -static const unsigned long rk_sobol_primitive_polynomials[] = { - 0x1UL, 0x3UL, 0x7UL, 0xBUL, 0xDUL, 0x13UL, 0x19UL, 0x25UL, 0x29UL, - 0x2FUL, 0x37UL, 0x3BUL, 0x3DUL, 0x43UL, 0x5BUL, 0x61UL, 0x67UL, 0x6DUL, - 0x73UL, 0x83UL, 0x89UL, 0x8FUL, 0x91UL, 0x9DUL, 0xA7UL, 0xABUL, 0xB9UL, - 0xBFUL, 0xC1UL, 0xCBUL, 0xD3UL, 0xD5UL, 0xE5UL, 0xEFUL, 0xF1UL, 0xF7UL, - 0xFDUL, 0x11DUL, 0x12BUL, 0x12DUL, 0x14DUL, 0x15FUL, 0x163UL, 0x165UL, - 0x169UL, 0x171UL, 0x187UL, 0x18DUL, 0x1A9UL, 0x1C3UL, 0x1CFUL, 0x1E7UL, - 0x1F5UL, 0x211UL, 0x21BUL, 0x221UL, 0x22DUL, 0x233UL, 0x259UL, 0x25FUL, - 0x269UL, 0x26FUL, 0x277UL, 0x27DUL, 0x287UL, 0x295UL, 0x2A3UL, 0x2A5UL, - 0x2AFUL, 0x2B7UL, 0x2BDUL, 0x2CFUL, 0x2D1UL, 0x2DBUL, 0x2F5UL, 0x2F9UL, - 0x313UL, 0x315UL, 0x31FUL, 0x323UL, 0x331UL, 0x33BUL, 0x34FUL, 0x35BUL, - 0x361UL, 0x36BUL, 0x36DUL, 0x373UL, 0x37FUL, 0x385UL, 0x38FUL, 0x3B5UL, - 0x3B9UL, 0x3C7UL, 0x3CBUL, 0x3CDUL, 0x3D5UL, 0x3D9UL, 0x3E3UL, 0x3E9UL, - 0x3FBUL, 0x409UL, 0x41BUL, 0x427UL, 0x42DUL, 0x465UL, 0x46FUL, 0x481UL, - 0x48BUL, 0x4C5UL, 0x4D7UL, 0x4E7UL, 0x4F3UL, 0x4FFUL, 0x50DUL, 0x519UL, - 0x523UL, 0x531UL, 0x53DUL, 0x543UL, 0x557UL, 0x56BUL, 0x585UL, 0x58FUL, - 0x597UL, 0x5A1UL, 0x5C7UL, 0x5E5UL, 0x5F7UL, 0x5FBUL, 0x613UL, 0x615UL, - 0x625UL, 0x637UL, 0x643UL, 0x64FUL, 0x65BUL, 0x679UL, 0x67FUL, 0x689UL, - 0x6B5UL, 0x6C1UL, 0x6D3UL, 0x6DFUL, 0x6FDUL, 0x717UL, 0x71DUL, 0x721UL, - 0x739UL, 0x747UL, 0x74DUL, 0x755UL, 0x759UL, 0x763UL, 0x77DUL, 0x78DUL, - 0x793UL, 0x7B1UL, 0x7DBUL, 0x7F3UL, 0x7F9UL, 0x805UL, 0x817UL, 0x82BUL, - 0x82DUL, 0x847UL, 0x863UL, 0x865UL, 0x871UL, 0x87BUL, 0x88DUL, 0x895UL, - 0x89FUL, 0x8A9UL, 0x8B1UL, 0x8CFUL, 0x8D1UL, 0x8E1UL, 0x8E7UL, 0x8EBUL, - 0x8F5UL, 0x90DUL, 0x913UL, 0x925UL, 0x929UL, 0x93BUL, 0x93DUL, 0x945UL, - 0x949UL, 0x951UL, 0x95BUL, 0x973UL, 0x975UL, 0x97FUL, 0x983UL, 0x98FUL, - 0x9ABUL, 0x9ADUL, 0x9B9UL, 0x9C7UL, 0x9D9UL, 0x9E5UL, 0x9F7UL, 0xA01UL, - 0xA07UL, 0xA13UL, 0xA15UL, 0xA29UL, 0xA49UL, 0xA61UL, 0xA6DUL, 0xA79UL, - 0xA7FUL, 0xA85UL, 0xA91UL, 0xA9DUL, 0xAA7UL, 0xAABUL, 0xAB3UL, 0xAB5UL, - 0xAD5UL, 0xADFUL, 0xAE9UL, 0xAEFUL, 0xAF1UL, 0xAFBUL, 0xB03UL, 0xB09UL, - 0xB11UL, 0xB33UL, 0xB3FUL, 0xB41UL, 0xB4BUL, 0xB59UL, 0xB5FUL, 0xB65UL, - 0xB6FUL, 0xB7DUL, 0xB87UL, 0xB8BUL, 0xB93UL, 0xB95UL, 0xBAFUL, 0xBB7UL, - 0xBBDUL, 0xBC9UL, 0xBDBUL, 0xBDDUL, 0xBE7UL, 0xBEDUL, 0xC0BUL, 0xC0DUL, - 0xC19UL, 0xC1FUL, 0xC57UL, 0xC61UL, 0xC6BUL, 0xC73UL, 0xC85UL, 0xC89UL, - 0xC97UL, 0xC9BUL, 0xC9DUL, 0xCB3UL, 0xCBFUL, 0xCC7UL, 0xCCDUL, 0xCD3UL, - 0xCD5UL, 0xCE3UL, 0xCE9UL, 0xCF7UL, 0xD03UL, 0xD0FUL, 0xD1DUL, 0xD27UL, - 0xD2DUL, 0xD41UL, 0xD47UL, 0xD55UL, 0xD59UL, 0xD63UL, 0xD6FUL, 0xD71UL, - 0xD93UL, 0xD9FUL, 0xDA9UL, 0xDBBUL, 0xDBDUL, 0xDC9UL, 0xDD7UL, 0xDDBUL, - 0xDE1UL, 0xDE7UL, 0xDF5UL, 0xE05UL, 0xE1DUL, 0xE21UL, 0xE27UL, 0xE2BUL, - 0xE33UL, 0xE39UL, 0xE47UL, 0xE4BUL, 0xE55UL, 0xE5FUL, 0xE71UL, 0xE7BUL, - 0xE7DUL, 0xE81UL, 0xE93UL, 0xE9FUL, 0xEA3UL, 0xEBBUL, 0xECFUL, 0xEDDUL, - 0xEF3UL, 0xEF9UL, 0xF0BUL, 0xF19UL, 0xF31UL, 0xF37UL, 0xF5DUL, 0xF6BUL, - 0xF6DUL, 0xF75UL, 0xF83UL, 0xF91UL, 0xF97UL, 0xF9BUL, 0xFA7UL, 0xFADUL, - 0xFB5UL, 0xFCDUL, 0xFD3UL, 0xFE5UL, 0xFE9UL, 0x1053UL, 0x1069UL, - 0x107BUL, 0x107DUL, 0x1099UL, 0x10D1UL, 0x10EBUL, 0x1107UL, 0x111FUL, - 0x1123UL, 0x113BUL, 0x114FUL, 0x1157UL, 0x1161UL, 0x116BUL, 0x1185UL, - 0x11B3UL, 0x11D9UL, 0x11DFUL, 0x120DUL, 0x1237UL, 0x123DUL, 0x1267UL, - 0x1273UL, 0x127FUL, 0x12B9UL, 0x12C1UL, 0x12CBUL, 0x130FUL, 0x131DUL, - 0x1321UL, 0x1339UL, 0x133FUL, 0x134DUL, 0x1371UL, 0x1399UL, 0x13A3UL, - 0x13A9UL, 0x1407UL, 0x1431UL, 0x1437UL, 0x144FUL, 0x145DUL, 0x1467UL, - 0x1475UL, 0x14A7UL, 0x14ADUL, 0x14D3UL, 0x150FUL, 0x151DUL, 0x154DUL, - 0x1593UL, 0x15C5UL, 0x15D7UL, 0x15DDUL, 0x15EBUL, 0x1609UL, 0x1647UL, - 0x1655UL, 0x1659UL, 0x16A5UL, 0x16BDUL, 0x1715UL, 0x1719UL, 0x1743UL, - 0x1745UL, 0x1775UL, 0x1789UL, 0x17ADUL, 0x17B3UL, 0x17BFUL, 0x17C1UL, - 0x1857UL, 0x185DUL, 0x1891UL, 0x1897UL, 0x18B9UL, 0x18EFUL, 0x191BUL, - 0x1935UL, 0x1941UL, 0x1965UL, 0x197BUL, 0x198BUL, 0x19B1UL, 0x19BDUL, - 0x19C9UL, 0x19CFUL, 0x19E7UL, 0x1A1BUL, 0x1A2BUL, 0x1A33UL, 0x1A69UL, - 0x1A8BUL, 0x1AD1UL, 0x1AE1UL, 0x1AF5UL, 0x1B0BUL, 0x1B13UL, 0x1B1FUL, - 0x1B57UL, 0x1B91UL, 0x1BA7UL, 0x1BBFUL, 0x1BC1UL, 0x1BD3UL, 0x1C05UL, - 0x1C11UL, 0x1C17UL, 0x1C27UL, 0x1C4DUL, 0x1C87UL, 0x1C9FUL, 0x1CA5UL, - 0x1CBBUL, 0x1CC5UL, 0x1CC9UL, 0x1CCFUL, 0x1CF3UL, 0x1D07UL, 0x1D23UL, - 0x1D43UL, 0x1D51UL, 0x1D5BUL, 0x1D75UL, 0x1D85UL, 0x1D89UL, 0x1E15UL, - 0x1E19UL, 0x1E2FUL, 0x1E45UL, 0x1E51UL, 0x1E67UL, 0x1E73UL, 0x1E8FUL, - 0x1EE3UL, 0x1F11UL, 0x1F1BUL, 0x1F27UL, 0x1F71UL, 0x1F99UL, 0x1FBBUL, - 0x1FBDUL, 0x1FC9UL, 0x201BUL, 0x2027UL, 0x2035UL, 0x2053UL, 0x2065UL, - 0x206FUL, 0x208BUL, 0x208DUL, 0x209FUL, 0x20A5UL, 0x20AFUL, 0x20BBUL, - 0x20BDUL, 0x20C3UL, 0x20C9UL, 0x20E1UL, 0x20F3UL, 0x210DUL, 0x2115UL, - 0x2129UL, 0x212FUL, 0x213BUL, 0x2143UL, 0x2167UL, 0x216BUL, 0x2179UL, - 0x2189UL, 0x2197UL, 0x219DUL, 0x21BFUL, 0x21C1UL, 0x21C7UL, 0x21CDUL, - 0x21DFUL, 0x21E3UL, 0x21F1UL, 0x21FBUL, 0x2219UL, 0x2225UL, 0x2237UL, - 0x223DUL, 0x2243UL, 0x225BUL, 0x225DUL, 0x2279UL, 0x227FUL, 0x2289UL, - 0x2297UL, 0x229BUL, 0x22B3UL, 0x22BFUL, 0x22CDUL, 0x22EFUL, 0x22F7UL, - 0x22FBUL, 0x2305UL, 0x2327UL, 0x232BUL, 0x2347UL, 0x2355UL, 0x2359UL, - 0x236FUL, 0x2371UL, 0x237DUL, 0x2387UL, 0x238DUL, 0x2395UL, 0x23A3UL, - 0x23A9UL, 0x23B1UL, 0x23B7UL, 0x23BBUL, 0x23E1UL, 0x23EDUL, 0x23F9UL, - 0x240BUL, 0x2413UL, 0x241FUL, 0x2425UL, 0x2429UL, 0x243DUL, 0x2451UL, - 0x2457UL, 0x2461UL, 0x246DUL, 0x247FUL, 0x2483UL, 0x249BUL, 0x249DUL, - 0x24B5UL, 0x24BFUL, 0x24C1UL, 0x24C7UL, 0x24CBUL, 0x24E3UL, 0x2509UL, - 0x2517UL, 0x251DUL, 0x2521UL, 0x252DUL, 0x2539UL, 0x2553UL, 0x2555UL, - 0x2563UL, 0x2571UL, 0x2577UL, 0x2587UL, 0x258BUL, 0x2595UL, 0x2599UL, - 0x259FUL, 0x25AFUL, 0x25BDUL, 0x25C5UL, 0x25CFUL, 0x25D7UL, 0x25EBUL, - 0x2603UL, 0x2605UL, 0x2611UL, 0x262DUL, 0x263FUL, 0x264BUL, 0x2653UL, - 0x2659UL, 0x2669UL, 0x2677UL, 0x267BUL, 0x2687UL, 0x2693UL, 0x2699UL, - 0x26B1UL, 0x26B7UL, 0x26BDUL, 0x26C3UL, 0x26EBUL, 0x26F5UL, 0x2713UL, - 0x2729UL, 0x273BUL, 0x274FUL, 0x2757UL, 0x275DUL, 0x276BUL, 0x2773UL, - 0x2779UL, 0x2783UL, 0x2791UL, 0x27A1UL, 0x27B9UL, 0x27C7UL, 0x27CBUL, - 0x27DFUL, 0x27EFUL, 0x27F1UL, 0x2807UL, 0x2819UL, 0x281FUL, 0x2823UL, - 0x2831UL, 0x283BUL, 0x283DUL, 0x2845UL, 0x2867UL, 0x2875UL, 0x2885UL, - 0x28ABUL, 0x28ADUL, 0x28BFUL, 0x28CDUL, 0x28D5UL, 0x28DFUL, 0x28E3UL, - 0x28E9UL, 0x28FBUL, 0x2909UL, 0x290FUL, 0x2911UL, 0x291BUL, 0x292BUL, - 0x2935UL, 0x293FUL, 0x2941UL, 0x294BUL, 0x2955UL, 0x2977UL, 0x297DUL, - 0x2981UL, 0x2993UL, 0x299FUL, 0x29AFUL, 0x29B7UL, 0x29BDUL, 0x29C3UL, - 0x29D7UL, 0x29F3UL, 0x29F5UL, 0x2A03UL, 0x2A0FUL, 0x2A1DUL, 0x2A21UL, - 0x2A33UL, 0x2A35UL, 0x2A4DUL, 0x2A69UL, 0x2A6FUL, 0x2A71UL, 0x2A7BUL, - 0x2A7DUL, 0x2AA5UL, 0x2AA9UL, 0x2AB1UL, 0x2AC5UL, 0x2AD7UL, 0x2ADBUL, - 0x2AEBUL, 0x2AF3UL, 0x2B01UL, 0x2B15UL, 0x2B23UL, 0x2B25UL, 0x2B2FUL, - 0x2B37UL, 0x2B43UL, 0x2B49UL, 0x2B6DUL, 0x2B7FUL, 0x2B85UL, 0x2B97UL, - 0x2B9BUL, 0x2BADUL, 0x2BB3UL, 0x2BD9UL, 0x2BE5UL, 0x2BFDUL, 0x2C0FUL, - 0x2C21UL, 0x2C2BUL, 0x2C2DUL, 0x2C3FUL, 0x2C41UL, 0x2C4DUL, 0x2C71UL, - 0x2C8BUL, 0x2C8DUL, 0x2C95UL, 0x2CA3UL, 0x2CAFUL, 0x2CBDUL, 0x2CC5UL, - 0x2CD1UL, 0x2CD7UL, 0x2CE1UL, 0x2CE7UL, 0x2CEBUL, 0x2D0DUL, 0x2D19UL, - 0x2D29UL, 0x2D2FUL, 0x2D37UL, 0x2D3BUL, 0x2D45UL, 0x2D5BUL, 0x2D67UL, - 0x2D75UL, 0x2D89UL, 0x2D8FUL, 0x2DA7UL, 0x2DABUL, 0x2DB5UL, 0x2DE3UL, - 0x2DF1UL, 0x2DFDUL, 0x2E07UL, 0x2E13UL, 0x2E15UL, 0x2E29UL, 0x2E49UL, - 0x2E4FUL, 0x2E5BUL, 0x2E5DUL, 0x2E61UL, 0x2E6BUL, 0x2E8FUL, 0x2E91UL, - 0x2E97UL, 0x2E9DUL, 0x2EABUL, 0x2EB3UL, 0x2EB9UL, 0x2EDFUL, 0x2EFBUL, - 0x2EFDUL, 0x2F05UL, 0x2F09UL, 0x2F11UL, 0x2F17UL, 0x2F3FUL, 0x2F41UL, - 0x2F4BUL, 0x2F4DUL, 0x2F59UL, 0x2F5FUL, 0x2F65UL, 0x2F69UL, 0x2F95UL, - 0x2FA5UL, 0x2FAFUL, 0x2FB1UL, 0x2FCFUL, 0x2FDDUL, 0x2FE7UL, 0x2FEDUL, - 0x2FF5UL, 0x2FFFUL, 0x3007UL, 0x3015UL, 0x3019UL, 0x302FUL, 0x3049UL, - 0x304FUL, 0x3067UL, 0x3079UL, 0x307FUL, 0x3091UL, 0x30A1UL, 0x30B5UL, - 0x30BFUL, 0x30C1UL, 0x30D3UL, 0x30D9UL, 0x30E5UL, 0x30EFUL, 0x3105UL, - 0x310FUL, 0x3135UL, 0x3147UL, 0x314DUL, 0x315FUL, 0x3163UL, 0x3171UL, - 0x317BUL, 0x31A3UL, 0x31A9UL, 0x31B7UL, 0x31C5UL, 0x31C9UL, 0x31DBUL, - 0x31E1UL, 0x31EBUL, 0x31EDUL, 0x31F3UL, 0x31FFUL, 0x3209UL, 0x320FUL, - 0x321DUL, 0x3227UL, 0x3239UL, 0x324BUL, 0x3253UL, 0x3259UL, 0x3265UL, - 0x3281UL, 0x3293UL, 0x3299UL, 0x329FUL, 0x32A9UL, 0x32B7UL, 0x32BBUL, - 0x32C3UL, 0x32D7UL, 0x32DBUL, 0x32E7UL, 0x3307UL, 0x3315UL, 0x332FUL, - 0x3351UL, 0x335DUL, 0x3375UL, 0x3397UL, 0x339BUL, 0x33ABUL, 0x33B9UL, - 0x33C1UL, 0x33C7UL, 0x33D5UL, 0x33E3UL, 0x33E5UL, 0x33F7UL, 0x33FBUL, - 0x3409UL, 0x341BUL, 0x3427UL, 0x3441UL, 0x344DUL, 0x345FUL, 0x3469UL, - 0x3477UL, 0x347BUL, 0x3487UL, 0x3493UL, 0x3499UL, 0x34A5UL, 0x34BDUL, - 0x34C9UL, 0x34DBUL, 0x34E7UL, 0x34F9UL, 0x350DUL, 0x351FUL, 0x3525UL, - 0x3531UL, 0x3537UL, 0x3545UL, 0x354FUL, 0x355DUL, 0x356DUL, 0x3573UL, - 0x357FUL, 0x359DUL, 0x35A1UL, 0x35B9UL, 0x35CDUL, 0x35D5UL, 0x35D9UL, - 0x35E3UL, 0x35E9UL, 0x35EFUL, 0x3601UL, 0x360BUL, 0x361FUL, 0x3625UL, - 0x362FUL, 0x363BUL, 0x3649UL, 0x3651UL, 0x365BUL, 0x3673UL, 0x3675UL, - 0x3691UL, 0x369BUL, 0x369DUL, 0x36ADUL, 0x36CBUL, 0x36D3UL, 0x36D5UL, - 0x36E3UL, 0x36EFUL, 0x3705UL, 0x370FUL, 0x371BUL, 0x3721UL, 0x372DUL, - 0x3739UL, 0x3741UL, 0x3747UL, 0x3753UL, 0x3771UL, 0x3777UL, 0x378BUL, - 0x3795UL, 0x3799UL, 0x37A3UL, 0x37C5UL, 0x37CFUL, 0x37D1UL, 0x37D7UL, - 0x37DDUL, 0x37E1UL, 0x37F3UL, 0x3803UL, 0x3805UL, 0x3817UL, 0x381DUL, - 0x3827UL, 0x3833UL, 0x384BUL, 0x3859UL, 0x3869UL, 0x3871UL, 0x38A3UL, - 0x38B1UL, 0x38BBUL, 0x38C9UL, 0x38CFUL, 0x38E1UL, 0x38F3UL, 0x38F9UL, - 0x3901UL, 0x3907UL, 0x390BUL, 0x3913UL, 0x3931UL, 0x394FUL, 0x3967UL, - 0x396DUL, 0x3983UL, 0x3985UL, 0x3997UL, 0x39A1UL, 0x39A7UL, 0x39ADUL, - 0x39CBUL, 0x39CDUL, 0x39D3UL, 0x39EFUL, 0x39F7UL, 0x39FDUL, 0x3A07UL, - 0x3A29UL, 0x3A2FUL, 0x3A3DUL, 0x3A51UL, 0x3A5DUL, 0x3A61UL, 0x3A67UL, - 0x3A73UL, 0x3A75UL, 0x3A89UL, 0x3AB9UL, 0x3ABFUL, 0x3ACDUL, 0x3AD3UL, - 0x3AD5UL, 0x3ADFUL, 0x3AE5UL, 0x3AE9UL, 0x3AFBUL, 0x3B11UL, 0x3B2BUL, - 0x3B2DUL, 0x3B35UL, 0x3B3FUL, 0x3B53UL, 0x3B59UL, 0x3B63UL, 0x3B65UL, - 0x3B6FUL, 0x3B71UL, 0x3B77UL, 0x3B8BUL, 0x3B99UL, 0x3BA5UL, 0x3BA9UL, - 0x3BB7UL, 0x3BBBUL, 0x3BD1UL, 0x3BE7UL, 0x3BF3UL, 0x3BFFUL, 0x3C0DUL, - 0x3C13UL, 0x3C15UL, 0x3C1FUL, 0x3C23UL, 0x3C25UL, 0x3C3BUL, 0x3C4FUL, - 0x3C5DUL, 0x3C6DUL, 0x3C83UL, 0x3C8FUL, 0x3C9DUL, 0x3CA7UL, 0x3CABUL, - 0x3CB9UL, 0x3CC7UL, 0x3CE9UL, 0x3CFBUL, 0x3CFDUL, 0x3D03UL, 0x3D17UL, - 0x3D1BUL, 0x3D21UL, 0x3D2DUL, 0x3D33UL, 0x3D35UL, 0x3D41UL, 0x3D4DUL, - 0x3D65UL, 0x3D69UL, 0x3D7DUL, 0x3D81UL, 0x3D95UL, 0x3DB1UL, 0x3DB7UL, - 0x3DC3UL, 0x3DD1UL, 0x3DDBUL, 0x3DE7UL, 0x3DEBUL, 0x3DF9UL, 0x3E05UL, - 0x3E09UL, 0x3E0FUL, 0x3E1BUL, 0x3E2BUL, 0x3E3FUL, 0x3E41UL, 0x3E53UL, - 0x3E65UL, 0x3E69UL, 0x3E8BUL, 0x3EA3UL, 0x3EBDUL, 0x3EC5UL, 0x3ED7UL, - 0x3EDDUL, 0x3EE1UL, 0x3EF9UL, 0x3F0DUL, 0x3F19UL, 0x3F1FUL, 0x3F25UL, - 0x3F37UL, 0x3F3DUL, 0x3F43UL, 0x3F45UL, 0x3F49UL, 0x3F51UL, 0x3F57UL, - 0x3F61UL, 0x3F83UL, 0x3F89UL, 0x3F91UL, 0x3FABUL, 0x3FB5UL, 0x3FE3UL, - 0x3FF7UL, 0x3FFDUL, - 0UL -}; - -rk_sobol_error rk_sobol_init(size_t dimension, rk_sobol_state *s, - rk_state *rs_dir, const unsigned long *directions, - const unsigned long *polynomials) -{ - rk_state rs_dir_temp; - int j, l, degree = 0, last_degree = 0, ooord = 0; - size_t k, cdir = 0, cpol = 0; - unsigned long polynomial = 1, rev = 0, last = 0; - - if (dimension == 0) - return RK_SOBOL_EINVAL; - - if (polynomials == NULL) - polynomials = rk_sobol_primitive_polynomials; - - /* Allocate the structure */ - s->direction = NULL; s->numerator = NULL; - s->direction = malloc(sizeof(*(s->direction))*dimension*LONG_BIT); - s->numerator = malloc(sizeof(*(s->numerator))*dimension); - if (!s->direction | !s->numerator) - { - if (!s->direction) free(s->direction); - if (!s->numerator) free(s->numerator); - return RK_SOBOL_ENOMEM; - } - - /* Initialize directions */ - /* Degree 0 */ - for (j = degree; j < LONG_BIT; j++) - s->direction[j*dimension] = 1UL << (LONG_BIT-j-1); - - /* Skip unused first polynomial */ - if (polynomials[cpol]) - cpol++; - - /* Degree >0 */ - for (k = 1; k < dimension; k++) - { - unsigned long temp; - - /* Find a new primitive polynomial */ - if (polynomials[cpol]) - polynomial = polynomials[cpol++]; - else if (rev) - { - /* We are generating polynomials out of order: - use the reverse of the previous polynomial */ - last = polynomial; - polynomial = rev; - rev = 0; - } - else - { - if (last) - { - polynomial = last; - last = 0; - } - - /* Find a new primitive polynomial */ - while(1) - { - if (polynomial == ULONG_MAX) - { - /* Not enough polynomials */ - free(s->direction); - free(s->numerator); - return RK_SOBOL_EINVAL; - } - - polynomial += 2; - - if (ooord) - { - unsigned long copy = polynomial; - /* We are generating polynomials out of order: - check if the reverse was already checked */ - for (rev = 0; copy; copy >>= 1) - rev = (rev << 1) | (copy & 1); - if (ooord && rev < polynomial) - continue; - } - - if (rk_isprimitive(polynomial)) - break; - } - - if (rev == polynomial) - /* We are generating polynomials out of order: - the reverse is not different, discard it */ - rev = 0; - } - - /* Compute the degree */ - for (temp = polynomial >> 1, degree = 0; temp; degree++, temp >>= 1); - - for (j=0; jdirection[j*dimension+k] = m << (LONG_BIT-j-1); - } - - /* Scaled recursion for directions */ - for (j = degree; j < LONG_BIT; j++) - { - unsigned long effdir = s->direction[(j-degree)*dimension+k], - ptemp = polynomial >> 1; - effdir ^= (effdir >> degree); - for (l = degree-1; l >= 1; l--, ptemp >>= 1) - if (ptemp & 1) - effdir ^= s->direction[(j-l)*dimension+k]; - s->direction[j*dimension+k] = effdir; - } - - /* Can we generate polynomials out of order ? */ - if (!ooord && polynomials[cpol] == 0 && degree > last_degree - && (directions == NULL || directions[cdir] == 0)) - ooord = 0; - else - last_degree = degree; - } - - /* Initialize numerator */ - for (k=0; knumerator[k] = 0; - - s->dimension = dimension; - s->gcount = 0; - s->count = 0; - return RK_SOBOL_OK; -} - -void rk_sobol_reinit(rk_sobol_state *s) -{ - size_t k; - - /* Initialize numerator */ - for (k=0; kdimension; k++) - s->numerator[k] = 0; - - s->count = 0; - s->gcount = 0; -} - -void rk_sobol_randomshift(rk_sobol_state *s, rk_state *rs_num) -{ - rk_state rs_num_temp; - size_t k; - - if (rs_num == NULL) - { - rs_num = &rs_num_temp; - rk_randomseed(rs_num); - } - - /* Initialize numerator */ - for (k=0; kdimension; k++) - s->numerator[k] = rk_ulong(rs_num); -} - -rk_sobol_error rk_sobol_copy(rk_sobol_state *copy, rk_sobol_state *orig) -{ - size_t k; - - /* Allocate the structure */ - copy->direction = NULL; copy->numerator = NULL; - copy->direction = malloc(sizeof(*(copy->direction))*orig->dimension*LONG_BIT); - copy->numerator = malloc(sizeof(*(copy->numerator))*orig->dimension); - if (!copy->direction | !copy->numerator) - { - if (!copy->direction) free(copy->direction); - if (!copy->numerator) free(copy->numerator); - return RK_SOBOL_ENOMEM; - } - - /* Initialize numerator */ - for (k=0; kdimension; k++) - copy->numerator[k] = orig->numerator[k]; - for (k=0; k<(orig->dimension*LONG_BIT); k++) - copy->direction[k] = orig->direction[k]; - - copy->count = orig->count; - copy->gcount = orig->gcount; - copy->dimension = orig->dimension; - - return RK_SOBOL_OK; -} - -rk_sobol_error rk_sobol_double(rk_sobol_state *s, double *x) -{ - int j; - size_t k; - unsigned long im; - const double inverse_denominator=1.0/(ULONG_MAX+1.0); - - if (s->count == ULONG_MAX) - j = 0; - else - for (im = s->count, j=0; im & 1; j++, im >>= 1); - s->count++; - - for (k=0; kdimension; k++) - { - s->numerator[k] ^= s->direction[j*s->dimension+k]; - x[k] = s->numerator[k]*inverse_denominator; - } - - if ((s->gcount++) == ULONG_MAX) return RK_SOBOL_EXHAUST; - return RK_SOBOL_OK; -} - -void rk_sobol_setcount(rk_sobol_state *s, unsigned long count) -{ - s->count = count; -} - -void rk_sobol_free(rk_sobol_state *s) -{ - free(s->direction); - free(s->numerator); -} - -double inverse_normal(double p) -{ - double q, t, x; - - /* Peter J. Acklam constants for the rational approximation */ - const double a[6] = - { - -3.969683028665376e+01, 2.209460984245205e+02, - -2.759285104469687e+02, 1.383577518672690e+02, - -3.066479806614716e+01, 2.506628277459239e+00 - }; - const double b[5] = - { - -5.447609879822406e+01, 1.615858368580409e+02, - -1.556989798598866e+02, 6.680131188771972e+01, - -1.328068155288572e+01 - }; - const double c[6] = - { - -7.784894002430293e-03, -3.223964580411365e-01, - -2.400758277161838e+00, -2.549732539343734e+00, - 4.374664141464968e+00, 2.938163982698783e+00 - }; - const double d[4] = - { - 7.784695709041462e-03, 3.224671290700398e-01, - 2.445134137142996e+00, 3.754408661907416e+00 - }; - - if (p <= 0) - return -HUGE_VAL; - else if (p >= 1) - return HUGE_VAL; - - q = p<0.5 ? p : 1-p; - if (q > 0.02425) - { - /* Rational approximation for central region */ - x = q-0.5; - t = x*x; - x = x*(((((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4])*t+a[5]) - /(((((b[0]*t+b[1])*t+b[2])*t+b[3])*t+b[4])*t+1); - } - else - { - /* Rational approximation for tail region */ - t = sqrt(-2*log(q)); - x = (((((c[0]*t+c[1])*t+c[2])*t+c[3])*t+c[4])*t+c[5]) - /((((d[0]*t+d[1])*t+d[2])*t+d[3])*t+1); - } - - /* If we have erfc, improve the precision */ -#ifndef WIN32 - /* Halley's rational method */ - t = (erfc(-x*M_SQRT1_2)/2 - q) * RK_SOBOL_M_SQRT2PI * exp(x*x/2); - x -= t/(1 + x*t/2); -#endif - - return p>0.5 ? -x : x; -} - -rk_sobol_error rk_sobol_gauss(rk_sobol_state *s, double *x) -{ - size_t k; - rk_sobol_error rc = rk_sobol_double(s, x); - - for (k=0; kdimension; k++) - x[k] = inverse_normal(x[k]); - - return rc; -} diff --git a/tardis/montecarlo/src/randomkit/rk_sobol.h b/tardis/montecarlo/src/randomkit/rk_sobol.h deleted file mode 100644 index aca71bebe1e..00000000000 --- a/tardis/montecarlo/src/randomkit/rk_sobol.h +++ /dev/null @@ -1,173 +0,0 @@ -/* Random kit 1.6 */ - -/* - * Copyright (c) 2004-2006, Jean-Sebastien Roy (js@jeannot.org) - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -/* @(#) $Jeannot: rk_sobol.h,v 1.7 2006/02/19 13:48:34 js Exp $ */ - -/* - * Typical use: - * - * int dimension = 2; - * rk_sobol_state s; - * rk_sobol_error rc; - * double x[dimension], y[dimension]; - * - * // Init - * if (rc = rk_sobol_init(dimension, &s, NULL, NULL, NULL)) - * { - * fprintf(stderr, "%s\n", rk_sobol_strerror[rc]); - * abort(); - * } - * - * // Draw uniform quasirandom doubles - * if (rc = rk_sobol_double(&s, x)) - * { - * fprintf(stderr, "%s\n", rk_sobol_strerror[rc]); - * abort(); - * } - * - * // Draw gaussian quasirandom doubles - * if (rc = rk_sobol_gauss(&s, y)) - * { - * fprintf(stderr, "%s\n", rk_sobol_strerror[rc]); - * abort(); - * } - * - * // Free allocated memory - * rk_sobol_free(&s); - */ - - -#ifndef _RK_SOBOL_ -#define _RK_SOBOL_ - -#include "rk_mt.h" - -typedef enum { - RK_SOBOL_OK = 0, /* No error */ - RK_SOBOL_EINVAL = 1, /* Invalid dimension (<= 0 or too large) */ - RK_SOBOL_EXHAUST = 2, /* Too many number generated */ - RK_SOBOL_ENOMEM = 3, /* Not enough memory */ - RK_SOBOL_ERR_MAX = 4 -} rk_sobol_error; - -/* error strings */ -extern char *rk_sobol_strerror[]; - -typedef struct -{ - size_t dimension; - unsigned long *direction; - unsigned long *numerator; - unsigned long count; - unsigned long gcount; -} rk_sobol_state; - -#ifdef __cplusplus -extern "C" { -#endif - -/* Sobol directions initializations (zero terminated lists) */ - -/* - * Sobol/Levitan coefficients of the free direction integers as given - * by Bratley, P., Fox, B.L. (1988) - * Defined up to dimension 40. - */ -extern const unsigned long rk_sobol_SLdirections[]; - -/* - * Lemieux coefficients of the free direction integers as given - * in QuantLib by Christiane Lemieux, private communication, September 2004 - * Defined up to dimension 360. - */ -extern const unsigned long rk_sobol_Ldirections[]; - -/* - * Peter J�ckel coefficients of the free direction integers as given - * in "Monte Carlo Methods in Finance", by Peter J�ckel, section 8.3 - * Defined up to dimension 32. - */ -extern const unsigned long rk_sobol_Jdirections[]; - -/* - * Initialize a sobol quasirandom number generator. - * 1 <= dimension <= the number of primitive polylonimals of degree < LONG_BIT - * If directions == NULL (or more directions than provided are required), - * the directions are picked at random using rs_dir. - * If rs_dir == NULL, it is initialized using rk_randomseed. - * polynomials is a zero terminated list of primitive polynomials to use if - * it is != NULL to speed up initialization for dimension > 1024. - */ -extern rk_sobol_error rk_sobol_init(size_t dimension, rk_sobol_state *s, - rk_state *rs_dir, const unsigned long *directions, - const unsigned long *polynomials); - -/* - * Reinitialize the random generator with same directions. - */ -extern void rk_sobol_reinit(rk_sobol_state *s); - -/* - * You can change the starting rank in the sequence by changing s->count. - */ -extern void rk_sobol_setcount(rk_sobol_state *s, unsigned long count); - -/* - * XOR the numerators at random using rs_num. - * To be used once, after (re-)initialization. - * Useful for randomized quasi monte carlo. - * If rs_num == NULL, it is initialized using rk_randomseed. - */ -extern void rk_sobol_randomshift(rk_sobol_state *s, rk_state *rs_num); - -/* - * Copy a sobol generator. - * Can be used to avoid the time consuming initialization. - */ -extern rk_sobol_error rk_sobol_copy(rk_sobol_state *copy, rk_sobol_state *orig); - -/* - * Free the memory allocated by rk_sobol_init - */ -extern void rk_sobol_free(rk_sobol_state *s); - -/* - * return a vector of dimension quasirandom uniform deviates between 0 and 1 - */ -extern rk_sobol_error rk_sobol_double(rk_sobol_state *s, double *x); - -/* - * return a vector of dimension quasirandom gaussian deviates - * with variance unity and zero mean. - * On Windows, the standard function erfc is missing, which results in - * lower precision (9 digits instead of full precision). - */ -extern rk_sobol_error rk_sobol_gauss(rk_sobol_state *s, double *x); - -#ifdef __cplusplus -} -#endif - -#endif /* _RK_SOBOL_ */ diff --git a/tardis/montecarlo/src/rpacket.c b/tardis/montecarlo/src/rpacket.c deleted file mode 100644 index cb52cddc698..00000000000 --- a/tardis/montecarlo/src/rpacket.c +++ /dev/null @@ -1,55 +0,0 @@ -#include -#include -#include "rpacket.h" -#include "storage.h" - -extern tardis_error_t line_search (const double *nu, double nu_insert, - int64_t number_of_lines, int64_t * result); - -tardis_error_t -rpacket_init (rpacket_t * packet, storage_model_t * storage, int packet_index, - int virtual_packet_flag, double * chi_bf_tmp_partial) -{ - int64_t current_line_id; - tardis_error_t ret_val = TARDIS_ERROR_OK; - double current_nu = storage->packet_nus[packet_index]; - double current_energy = storage->packet_energies[packet_index]; - double current_mu = storage->packet_mus[packet_index]; - double comov_current_nu = current_nu; - int current_shell_id = 0; - double current_r = storage->r_inner[0]; - double beta = current_r * storage->inverse_time_explosion * INVERSE_C; - - if (storage->full_relativity) - { - current_nu = current_nu * (1 + beta * current_mu) / sqrt(1 - beta * beta); - current_energy = current_energy * (1 + beta * current_mu) / sqrt(1 - beta * beta); - current_mu = (current_mu + beta) / (1 + beta * current_mu); - } - else - { - current_nu = current_nu / (1 - beta * current_mu); - current_energy = current_energy / (1 - beta * current_mu); - } - if ((ret_val = - line_search (storage->line_list_nu, comov_current_nu, - storage->no_of_lines, - ¤t_line_id)) != TARDIS_ERROR_OK) - { - return ret_val; - } - bool last_line = (current_line_id == storage->no_of_lines); - rpacket_set_nu (packet, current_nu); - rpacket_set_mu (packet, current_mu); - rpacket_set_energy (packet, current_energy); - rpacket_set_r (packet, current_r); - rpacket_set_current_shell_id (packet, current_shell_id); - rpacket_set_next_line_id (packet, current_line_id); - rpacket_set_last_line (packet, last_line); - rpacket_set_close_line (packet, false); - rpacket_set_virtual_packet_flag (packet, virtual_packet_flag); - packet->chi_bf_tmp_partial = chi_bf_tmp_partial; - packet->compute_chi_bf = true; - packet->vpacket_weight = 1.0; - return ret_val; -} diff --git a/tardis/montecarlo/src/rpacket.h b/tardis/montecarlo/src/rpacket.h deleted file mode 100644 index 2a30546a837..00000000000 --- a/tardis/montecarlo/src/rpacket.h +++ /dev/null @@ -1,330 +0,0 @@ -#ifndef TARDIS_RPACKET_H -#define TARDIS_RPACKET_H - -#include -#include -#include -#include "randomkit/randomkit.h" -#include "status.h" -#include "storage.h" - -#define MISS_DISTANCE 1e99 -#define C 29979245800.0 -#define INVERSE_C 3.33564095198152e-11 -#define H 6.6260755e-27 // erg * s, converted to CGS units from the NIST Constant Index -#define KB 1.3806488e-16 // erg / K converted to CGS units from the NIST Constant Index - -/** - * @brief A photon packet. - */ -typedef struct RPacket -{ - double nu; /**< Frequency of the packet in Hz. */ - double mu; /**< Cosine of the angle of the packet. */ - double energy; /**< Energy of the packet in erg. */ - double r; /**< Distance from center in cm. */ - double tau_event; - double nu_line; - int64_t current_shell_id; /**< ID of the current shell. */ - int64_t next_line_id; /**< The index of the next line that the packet will encounter. */ - /** - * @brief The packet has a nu red-ward of the last line. - * It will not encounter any lines anymore. - */ - int64_t last_line; - /** - * @brief The packet just encountered a line that is very close to the next line. - * The next iteration will automatically make an interaction with the next line - * (avoiding numerical problems). - */ - int64_t close_line; - /** - * @brief packet is a virtual packet and will ignore any d_line or d_electron checks. - * It now whenever a d_line is calculated only adds the tau_line to an - * internal float. - */ - int64_t current_continuum_id; /* Packet can interact with bf-continua with an index equal or bigger than this */ - int64_t virtual_packet_flag; - int64_t virtual_packet; - double d_line; /**< Distance to electron event. */ - double d_electron; /**< Distance to line event. */ - double d_boundary; /**< Distance to shell boundary. */ - double d_cont; /**< Distance to continuum event */ - int64_t next_shell_id; /**< ID of the next shell packet visits. */ - rpacket_status_t status; /**< Packet status (in process, emitted or reabsorbed). */ - int64_t id; - double chi_th; /**< Opacity due to electron scattering */ - double chi_cont; /**< Opacity due to continuum processes */ - double chi_ff; /**< Opacity due to free-free processes */ - double chi_bf; /**< Opacity due to bound-free processes */ - double *chi_bf_tmp_partial; - int64_t macro_atom_activation_level; - bool compute_chi_bf; - double vpacket_weight; -} rpacket_t; - -static inline double rpacket_get_nu (const rpacket_t * packet) -{ - return packet->nu; -} - -static inline void rpacket_set_nu (rpacket_t * packet, double nu) -{ - packet->nu = nu; -} - -static inline double rpacket_get_mu (const rpacket_t * packet) -{ - return packet->mu; -} - -static inline void rpacket_set_mu (rpacket_t * packet, double mu) -{ - packet->mu = mu; -} - -static inline double rpacket_get_energy (const rpacket_t * packet) -{ - return packet->energy; -} - -static inline void rpacket_set_energy (rpacket_t * packet, double energy) -{ - packet->energy = energy; -} - -static inline double rpacket_get_r (const rpacket_t * packet) -{ - return packet->r; -} - -static inline void rpacket_set_r (rpacket_t * packet, double r) -{ - packet->r = r; -} - -static inline double rpacket_get_tau_event (const rpacket_t * packet) -{ - return packet->tau_event; -} - -static inline void rpacket_set_tau_event (rpacket_t * packet, double tau_event) -{ - packet->tau_event = tau_event; -} - -static inline double rpacket_get_nu_line (const rpacket_t * packet) -{ - return packet->nu_line; -} - -static inline void rpacket_set_nu_line (rpacket_t * packet, double nu_line) -{ - packet->nu_line = nu_line; -} - -static inline unsigned int rpacket_get_current_shell_id (const rpacket_t * packet) -{ - return packet->current_shell_id; -} - -static inline void rpacket_set_current_shell_id (rpacket_t * packet, - unsigned int current_shell_id) -{ - packet->current_shell_id = current_shell_id; -} - -static inline unsigned int rpacket_get_next_line_id (const rpacket_t * packet) -{ - return packet->next_line_id; -} - -static inline void rpacket_set_next_line_id (rpacket_t * packet, - unsigned int next_line_id) -{ - packet->next_line_id = next_line_id; -} - -static inline bool rpacket_get_last_line (const rpacket_t * packet) -{ - return packet->last_line; -} - -static inline void rpacket_set_last_line (rpacket_t * packet, bool last_line) -{ - packet->last_line = last_line; -} - -static inline bool rpacket_get_close_line (const rpacket_t * packet) -{ - return packet->close_line; -} - -static inline void rpacket_set_close_line (rpacket_t * packet, bool close_line) -{ - packet->close_line = close_line; -} - -static inline int rpacket_get_virtual_packet_flag (const rpacket_t * packet) -{ - return packet->virtual_packet_flag; -} - -static inline void rpacket_set_virtual_packet_flag (rpacket_t * packet, - int virtual_packet_flag) -{ - packet->virtual_packet_flag = virtual_packet_flag; -} - -static inline int rpacket_get_virtual_packet (const rpacket_t * packet) -{ - return packet->virtual_packet; -} - -static inline void rpacket_set_virtual_packet (rpacket_t * packet, - int virtual_packet) -{ - packet->virtual_packet = virtual_packet; -} - -static inline double rpacket_get_d_boundary (const rpacket_t * packet) -{ - return packet->d_boundary; -} - -static inline void rpacket_set_d_boundary (rpacket_t * packet, double d_boundary) -{ - packet->d_boundary = d_boundary; -} - -static inline double rpacket_get_d_electron (const rpacket_t * packet) -{ - return packet->d_electron; -} - -static inline void rpacket_set_d_electron (rpacket_t * packet, double d_electron) -{ - packet->d_electron = d_electron; -} - -static inline double rpacket_get_d_line (const rpacket_t * packet) -{ - return packet->d_line; -} - -static inline void rpacket_set_d_line (rpacket_t * packet, double d_line) -{ - packet->d_line = d_line; -} - -static inline int rpacket_get_next_shell_id (const rpacket_t * packet) -{ - return packet->next_shell_id; -} - -static inline void rpacket_set_next_shell_id (rpacket_t * packet, int next_shell_id) -{ - packet->next_shell_id = next_shell_id; -} - -static inline rpacket_status_t rpacket_get_status (const rpacket_t * packet) -{ - return packet->status; -} - -static inline void rpacket_set_status (rpacket_t * packet, rpacket_status_t status) -{ - packet->status = status; -} - -static inline int rpacket_get_id (const rpacket_t * packet) -{ - return packet->id; -} - -static inline void rpacket_set_id (rpacket_t * packet, int id) -{ - packet->id = id; -} - -static inline void rpacket_reset_tau_event (rpacket_t * packet, rk_state *mt_state) -{ - rpacket_set_tau_event (packet, -log (rk_double (mt_state))); -} - -tardis_error_t rpacket_init (rpacket_t * packet, storage_model_t * storage, - int packet_index, int virtual_packet_flag, double * chi_bf_tmp_partial); - -/* New getter and setter methods for continuum implementation */ - -static inline void rpacket_set_d_continuum (rpacket_t * packet, double d_continuum) -{ - packet->d_cont = d_continuum; -} - -static inline double rpacket_get_d_continuum (const rpacket_t * packet) -{ - return packet->d_cont; -} - -static inline void rpacket_set_chi_electron (rpacket_t * packet, double chi_electron) -{ - packet->chi_th = chi_electron; -} - -static inline double rpacket_get_chi_electron (const rpacket_t * packet) -{ - return packet->chi_th; -} - -static inline void rpacket_set_chi_continuum (rpacket_t * packet, double chi_continuum) -{ - packet->chi_cont = chi_continuum; -} - -static inline double rpacket_get_chi_continuum (const rpacket_t * packet) -{ - return packet->chi_cont; -} - -static inline void rpacket_set_chi_freefree (rpacket_t * packet, double chi_freefree) -{ - packet->chi_ff = chi_freefree; -} - -static inline double rpacket_get_chi_freefree (const rpacket_t * packet) -{ - return packet->chi_ff; -} - -static inline void rpacket_set_chi_boundfree (rpacket_t * packet, double chi_boundfree) -{ - packet->chi_bf = chi_boundfree; -} - -static inline double rpacket_get_chi_boundfree (const rpacket_t * packet) -{ - return packet->chi_bf; -} - -static inline unsigned int rpacket_get_current_continuum_id (const rpacket_t * packet) -{ - return packet->current_continuum_id; -} - -static inline void rpacket_set_current_continuum_id (rpacket_t * packet, unsigned int current_continuum_id) -{ - packet->current_continuum_id = current_continuum_id; -} - -static inline void rpacket_set_macro_atom_activation_level (rpacket_t * packet, unsigned int activation_level) -{ - packet->macro_atom_activation_level = activation_level; -} - -static inline unsigned int rpacket_get_macro_atom_activation_level (const rpacket_t * packet) -{ - return packet->macro_atom_activation_level; -} - -#endif // TARDIS_RPACKET_H diff --git a/tardis/montecarlo/src/status.h b/tardis/montecarlo/src/status.h deleted file mode 100644 index 17a01720b79..00000000000 --- a/tardis/montecarlo/src/status.h +++ /dev/null @@ -1,38 +0,0 @@ -#ifndef TARDIS_STATUS_H -#define TARDIS_STATUS_H - -typedef enum -{ - TARDIS_ERROR_OK = 0, - TARDIS_ERROR_BOUNDS_ERROR = 1, - TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE = 2 -} tardis_error_t; - -typedef enum -{ - TARDIS_PACKET_STATUS_IN_PROCESS = 0, - TARDIS_PACKET_STATUS_EMITTED = 1, - TARDIS_PACKET_STATUS_REABSORBED = 2 -} rpacket_status_t; - -typedef enum -{ - CONTINUUM_OFF = 0, - CONTINUUM_ON = 1 -} ContinuumProcessesStatus; - -typedef enum -{ - BB_EMISSION = -1, - BF_EMISSION = -2, - FF_EMISSION = -3, - ADIABATIC_COOLING = -4 -} emission_type; - -typedef enum -{ - LIN_INTERPOLATION = 0, - HYDROGENIC = 1 -} bound_free_treatment; - -#endif // TARDIS_STATUS_H diff --git a/tardis/montecarlo/src/storage.h b/tardis/montecarlo/src/storage.h deleted file mode 100644 index a408a9f7bc2..00000000000 --- a/tardis/montecarlo/src/storage.h +++ /dev/null @@ -1,102 +0,0 @@ -#ifndef TARDIS_STORAGE_H -#define TARDIS_STORAGE_H - -#include - -#include "status.h" - -typedef struct photo_xsect_1level -{ - double * nu; - double * x_sect; - int64_t no_of_points; -} photo_xsect_1level; - -typedef struct StorageModel -{ - double *packet_nus; - double *packet_mus; - double *packet_energies; - double *output_nus; - double *output_energies; - double *last_interaction_in_nu; - int64_t *last_line_interaction_in_id; - int64_t *last_line_interaction_out_id; - int64_t *last_line_interaction_shell_id; - int64_t *last_interaction_type; - int64_t *last_interaction_out_type; - int64_t no_of_packets; - int64_t no_of_shells; - int64_t no_of_shells_i; - double *r_inner; - double *r_outer; - double *r_inner_i; - double *r_outer_i; - double *v_inner; - double time_explosion; - double inverse_time_explosion; - double *electron_densities; - double *electron_densities_i; - double *inverse_electron_densities; - double *line_list_nu; - double *continuum_list_nu; - double *line_lists_tau_sobolevs; - double *line_lists_tau_sobolevs_i; - int64_t line_lists_tau_sobolevs_nd; - double *line_lists_j_blues; - int64_t line_lists_j_blues_nd; - double *line_lists_Edotlu; - int64_t no_of_lines; - int64_t no_of_edges; - int64_t line_interaction_id; - double *transition_probabilities; - int64_t transition_probabilities_nd; - int64_t *line2macro_level_upper; - int64_t *macro_block_references; - int64_t *transition_type; - int64_t *destination_level_id; - int64_t *transition_line_id; - double *js; - double *nubars; - double spectrum_start_nu; - double spectrum_delta_nu; - double spectrum_end_nu; - double spectrum_virt_start_nu; - double spectrum_virt_end_nu; - double *spectrum_virt_nu; - double sigma_thomson; - double inverse_sigma_thomson; - double inner_boundary_albedo; - int64_t reflective_inner_boundary; - int64_t current_packet_id; - photo_xsect_1level **photo_xsect; - double *chi_ff_factor; - double *t_electrons; - double *l_pop; - double *l_pop_r; - ContinuumProcessesStatus cont_status; - bound_free_treatment bf_treatment; - double *virt_packet_nus; - double *virt_packet_energies; - double *virt_packet_last_interaction_in_nu; - int64_t *virt_packet_last_interaction_type; - int64_t *virt_packet_last_line_interaction_in_id; - int64_t *virt_packet_last_line_interaction_out_id; - int64_t virt_packet_count; - int64_t virt_array_size; - int64_t kpacket2macro_level; - int64_t *cont_edge2macro_level; - double *photo_ion_estimator; - double *stim_recomb_estimator; - int64_t *photo_ion_estimator_statistics; - double *bf_heating_estimator; - double *ff_heating_estimator; - double *stim_recomb_cooling_estimator; - int full_relativity; - double survival_probability; - double tau_russian; - double *tau_bias; - int enable_biasing; -} storage_model_t; - -#endif // TARDIS_STORAGE_H diff --git a/tardis/montecarlo/tests/conftest.py b/tardis/montecarlo/tests/conftest.py index 13cca00085e..054f6b455c4 100644 --- a/tardis/montecarlo/tests/conftest.py +++ b/tardis/montecarlo/tests/conftest.py @@ -9,7 +9,6 @@ c_ulong, ) -from tardis.montecarlo import montecarlo from tardis.montecarlo.struct import ( RPacket, StorageModel, RKState, TARDIS_PACKET_STATUS_IN_PROCESS, @@ -18,12 +17,6 @@ ) -# Wrap the shared object containing C methods, which are tested here. -@pytest.fixture(scope='session') -def clib(): - return CDLL(os.path.join(montecarlo.__file__)) - - @pytest.fixture(scope="function") def packet(): """Fixture to return `RPacket` object with default params initialized.""" diff --git a/tardis/montecarlo/tests/test_formal_integral.py b/tardis/montecarlo/tests/test_formal_integral.py index 616de0d2d7d..ccf52665e5c 100644 --- a/tardis/montecarlo/tests/test_formal_integral.py +++ b/tardis/montecarlo/tests/test_formal_integral.py @@ -2,23 +2,16 @@ import numpy as np from tardis import constants as c -import ctypes -from ctypes import ( - c_double, - c_int - ) - from numpy.ctypeslib import ( as_array, as_ctypes, - ndpointer ) import numpy.testing as ntest from tardis.util.base import intensity_black_body -from tardis.montecarlo.struct import StorageModel +import tardis.montecarlo.formal_integral as formal_integral @pytest.mark.parametrize( @@ -29,11 +22,10 @@ (1, 1), ] ) -def test_intensity_black_body(clib, nu, T): - func = clib.intensity_black_body - func.restype = c_double - func.argtypes = [c_double, c_double] +def test_intensity_black_body(nu, T): + func = formal_integral.intensity_black_body actual = func(nu, T) + print(actual, type(actual)) expected = intensity_black_body(nu, T) ntest.assert_almost_equal( actual, @@ -45,16 +37,10 @@ def test_intensity_black_body(clib, nu, T): 'N', (1e2, 1e3, 1e4, 1e5) ) -def test_trapezoid_integration(clib, N): - func = clib.trapezoid_integration - func.restype = c_double +def test_trapezoid_integration(N): + func = formal_integral.trapezoid_integration h = 1. N = int(N) - func.argtypes = [ - ndpointer(c_double), - c_double, - c_int - ] data = np.random.random(N) actual = func(data, h, int(N)) @@ -70,7 +56,7 @@ def test_trapezoid_integration(clib, N): True, reason='static inline functions are not inside the library' ) -def test_calculate_z(clib): +def test_calculate_z(): pass @@ -97,7 +83,7 @@ def calculate_z(r, p): def formal_integral_model(request, model): r = request.param['r'] model.no_of_shells_i = r.shape[0] - 1 - model.inverse_time_explosion = c.c.cgs.value + model.time_explosion = 1/c.c.cgs.value model.r_outer_i.contents = as_ctypes(r[1:]) model.r_inner_i.contents = as_ctypes(r[:-1]) return model @@ -107,20 +93,16 @@ def formal_integral_model(request, model): 'p', [0, 0.5, 1] ) -def test_populate_z_photosphere(clib, formal_integral_model, p): +def test_populate_z_photosphere(formal_integral_model, p): ''' Test the case where p < r[0] That means we 'hit' all shells from inside to outside. ''' - func = clib.populate_z - func.restype = ctypes.c_int64 - func.argtypes = [ - ctypes.POINTER(StorageModel), # storage - c_double, # p - ndpointer(dtype=np.float64), # oz - ndpointer(dtype=np.int64) # oshell_id - ] - + integrator = formal_integral.FormalIntegrator( + formal_integral_model, + None, + None) + func = integrator.populate_z size = formal_integral_model.no_of_shells_i r_inner = as_array(formal_integral_model.r_inner_i, (size,)) r_outer = as_array(formal_integral_model.r_outer_i, (size,)) @@ -130,7 +112,6 @@ def test_populate_z_photosphere(clib, formal_integral_model, p): oshell_id = np.zeros_like(oz, dtype=np.int64) N = func( - formal_integral_model, p, oz, oshell_id @@ -153,18 +134,15 @@ def test_populate_z_photosphere(clib, formal_integral_model, p): 'p', [1e-5, 0.5, 0.99, 1] ) -def test_populate_z_shells(clib, formal_integral_model, p): +def test_populate_z_shells(formal_integral_model, p): ''' Test the case where p > r[0] ''' - func = clib.populate_z - func.restype = ctypes.c_int64 - func.argtypes = [ - ctypes.POINTER(StorageModel), # storage - c_double, # p - ndpointer(dtype=np.float64), # oz - ndpointer(dtype=np.int64) # oshell_id - ] + integrator = formal_integral.FormalIntegrator( + formal_integral_model, + None, + None) + func = integrator.populate_z size = formal_integral_model.no_of_shells_i r_inner = as_array(formal_integral_model.r_inner_i, (size,)) @@ -194,7 +172,6 @@ def test_populate_z_shells(clib, formal_integral_model, p): p) N = func( - formal_integral_model, p, oz, oshell_id @@ -221,14 +198,10 @@ def test_populate_z_shells(clib, formal_integral_model, p): 1000, 10000, ]) -def test_calculate_p_values(clib, N): +def test_calculate_p_values(N): r = 1. - func = clib.calculate_p_values - func.argtypes = [ - c_double, - c_int, - ndpointer(dtype=np.float64) - ] + func = formal_integral.calculate_p_values + expected = r/(N - 1) * np.arange(0, N, dtype=np.float64) actual = np.zeros_like(expected, dtype=np.float64) diff --git a/tardis/montecarlo/tests/test_montecarlo.py b/tardis/montecarlo/tests/test_montecarlo.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tardis/scripts/debug/run_numba_single.py b/tardis/scripts/debug/run_numba_single.py new file mode 100644 index 00000000000..86d0b22f331 --- /dev/null +++ b/tardis/scripts/debug/run_numba_single.py @@ -0,0 +1,22 @@ +from tardis import run_tardis +import numpy as np +from tardis.montecarlo.montecarlo_numba.base import montecarlo_main_loop +import os +import numba +import sys +import yaml + + +SEED = eval(sys.argv[1].split('=')[1])[0] + +yaml_file, params = 'tardis_example_single.yml', None + +with open(yaml_file) as f: + params = yaml.safe_load(f) + +params['montecarlo']['single_packet_seed'] = SEED + +with open(yaml_file, 'w') as f: + yaml.safe_dump(params, f) + +mdl = run_tardis(yaml_file) diff --git a/tardis/scripts/debug/run_numba_single.run.xml b/tardis/scripts/debug/run_numba_single.run.xml new file mode 100644 index 00000000000..987a1e7808a --- /dev/null +++ b/tardis/scripts/debug/run_numba_single.run.xml @@ -0,0 +1,25 @@ + + + + + \ No newline at end of file diff --git a/tardis/scripts/debug/tardis_example_single.yml b/tardis/scripts/debug/tardis_example_single.yml new file mode 100644 index 00000000000..35c17071b18 --- /dev/null +++ b/tardis/scripts/debug/tardis_example_single.yml @@ -0,0 +1,45 @@ +atom_data: kurucz_cd23_chianti_H_He.h5 +model: + abundances: + Si: 1.0 + type: uniform + structure: + density: + type: branch85_w7 + type: specific + velocity: + num: 20 + start: 1.1e4 km/s + stop: 20000 km/s +montecarlo: + convergence_strategy: + damping_constant: 1.0 + fraction: 0.8 + hold_iterations: 3 + t_inner: + damping_constant: 1.0 + threshold: 0.05 + type: damped + debug_packets: true + iterations: 2 + last_no_of_packets: 1000000.0 + logger_buffer: 1 + no_of_packets: 40000.0 + no_of_virtual_packets: 0 + nthreads: 6 + seed: 23111963 + single_packet_seed: 46 +plasma: + disable_electron_scattering: false + excitation: lte + ionization: lte + line_interaction_type: macroatom + radiative_rates_type: dilute-blackbody +spectrum: + num: 10000 + start: 500 angstrom + stop: 20000 angstrom +supernova: + luminosity_requested: 9.44 log_lsun + time_explosion: 10 day +tardis_config_version: v1.0 diff --git a/tardis/simulation/base.py b/tardis/simulation/base.py index 8d9d0d30159..fd3161701bf 100644 --- a/tardis/simulation/base.py +++ b/tardis/simulation/base.py @@ -257,7 +257,7 @@ def advance_state(self): # case it needs some extra kwargs. if 'j_blue_estimator' in self.plasma.outputs_dict: update_properties.update(t_inner=next_t_inner, - j_blue_estimator=self.runner.j_b_lu_estimator) + j_blue_estimator=self.runner.j_blue_estimator) self.plasma.update(**update_properties) @@ -299,7 +299,7 @@ def run(self): self.model.t_rad, self.plasma.electron_densities, self.model.t_inner) - self.iterate(self.last_no_of_packets, self.no_of_virtual_packets, True) + self.iterate(self.last_no_of_packets, self.no_of_virtual_packets, last_run=True) self.reshape_plasma_state_store(self.iterations_executed) diff --git a/tardis/tests/test_tardis_full.py b/tardis/tests/test_tardis_full.py index 4e98249b685..682e4645728 100644 --- a/tardis/tests/test_tardis_full.py +++ b/tardis/tests/test_tardis_full.py @@ -52,7 +52,7 @@ def test_j_blue_estimators(self, runner, refdata): j_blue_estimator = refdata('j_blue_estimator').values npt.assert_allclose( - runner.j_b_lu_estimator, + runner.j_blue_estimator, j_blue_estimator) def test_spectrum(self, runner, refdata): diff --git a/tardis/tests/test_tardis_full_formal_integral.py b/tardis/tests/test_tardis_full_formal_integral.py index 3be437b8113..d3de68501c4 100644 --- a/tardis/tests/test_tardis_full_formal_integral.py +++ b/tardis/tests/test_tardis_full_formal_integral.py @@ -79,7 +79,7 @@ def test_j_blue_estimators(self, runner, refdata): j_blue_estimator = refdata('j_blue_estimator').values npt.assert_allclose( - runner.j_b_lu_estimator, + runner.j_blue_estimator, j_blue_estimator) def test_spectrum(self, runner, refdata): diff --git a/tardis_env3.yml b/tardis_env3.yml index 8f0539b7893..d9e7678fd6e 100644 --- a/tardis_env3.yml +++ b/tardis_env3.yml @@ -11,7 +11,7 @@ dependencies: - pandas=0.24 - astropy=3.2.1 - - numba=0.43 + - numba=0.49.1 - numexpr - Cython=0.29