From 328ec77df73f151f2f0c51dc9a851a3d96c1c37d Mon Sep 17 00:00:00 2001 From: Israel Roldan Date: Mon, 6 May 2024 08:28:31 -0600 Subject: [PATCH] GSoC 2024 | Benchmark | First objective | Israel Roldan (#2565) * Added the development environment to make easy the test of the changes. * Improved the ASV framework; Added benchmarks taking the unit tests. * Added the script for the run the ASV benchmarks by tag release in Git. * Added my conact information in the mail map. * Removed the Docker files. * Updated the release script which runs the benchmarks. * Updated the runtime timeout for the benchmarks. * Updated the function get relative paths for benchmarks. * Added comments for the deprecated function 'delim_whitespace'. * Fixed the usage of the regression data class. * Added proposal to divide the main functionality of the Benchmark Base class. * Update the benchmark for Gamma Ray. * Renamed the class name for NLTE in benchmarks. * Applied the 'black' formatter to main_gamma_ray_loop. * Revert "Applied the 'black' formatter to main_gamma_ray_loop." This reverts commit abe3b0d5b206e45156049788df3ac75e6f970fdf. * Added documentation for the benchmarks. * Fixed the format error warning by black. * Fixed again the format error warning by black. * Deleted all benchmarks except montecarlo numba and tardis; --------- Co-authored-by: Israel Roldan --- .gitignore | 9 + .mailmap | 4 + asv.conf.json | 2 +- benchmarks/asv_by_release.bash | 99 ++++++ benchmarks/benchmark_base.py | 329 ++++++++++++++++++ benchmarks/benchmark_run_tardis.py | 20 -- benchmarks/data/tardis_configv1_benchmark.yml | 2 +- ...montecarlo_montecarlo_numba_interaction.py | 103 ++++++ ...ontecarlo_numba_numba_formal_integral_p.py | 178 ++++++++++ ...ecarlo_montecarlo_numba_numba_interface.py | 77 ++++ .../montecarlo_montecarlo_numba_opacities.py | 91 +++++ .../montecarlo_montecarlo_numba_packet.py | 307 ++++++++++++++++ .../montecarlo_montecarlo_numba_r_packet.py | 67 ++++ .../montecarlo_montecarlo_numba_vpacket.py | 119 +++++++ benchmarks/run_tardis.py | 25 ++ benchmarks/util/__init__.py | 0 benchmarks/util/base.py | 20 ++ benchmarks/util/nlte.py | 78 +++++ docs/contributing/development/benchmarks.rst | 89 +++++ docs/contributing/development/index.rst | 1 + .../tardis_configv1_verysimple_logger.yml | 1 + tardis/io/model/readers/artis.py | 6 +- tardis/io/model/readers/blondin_toymodel.py | 14 +- tardis/io/model/readers/stella.py | 6 +- .../tardis_configv1_ascii_density_abund.yml | 1 - .../data/tardis_configv1_isotope_iabund.yml | 64 ++-- .../data/tardis_configv1_isotope_uniabund.yml | 2 +- .../tests/data/plasma_base_test_config.yml | 4 - tardis/util/base.py | 6 +- 29 files changed, 1657 insertions(+), 67 deletions(-) create mode 100755 benchmarks/asv_by_release.bash create mode 100644 benchmarks/benchmark_base.py delete mode 100644 benchmarks/benchmark_run_tardis.py create mode 100644 benchmarks/montecarlo_montecarlo_numba_interaction.py create mode 100644 benchmarks/montecarlo_montecarlo_numba_numba_formal_integral_p.py create mode 100644 benchmarks/montecarlo_montecarlo_numba_numba_interface.py create mode 100644 benchmarks/montecarlo_montecarlo_numba_opacities.py create mode 100644 benchmarks/montecarlo_montecarlo_numba_packet.py create mode 100644 benchmarks/montecarlo_montecarlo_numba_r_packet.py create mode 100644 benchmarks/montecarlo_montecarlo_numba_vpacket.py create mode 100644 benchmarks/run_tardis.py create mode 100644 benchmarks/util/__init__.py create mode 100644 benchmarks/util/base.py create mode 100644 benchmarks/util/nlte.py create mode 100644 docs/contributing/development/benchmarks.rst diff --git a/.gitignore b/.gitignore index 7b0c1a23a4a..4f02f9da6f8 100644 --- a/.gitignore +++ b/.gitignore @@ -15,6 +15,7 @@ __pycache__ */cython_version.py htmlcov .coverage +coverage.xml MANIFEST .ipynb_checkpoints @@ -77,3 +78,11 @@ pip-wheel-metadata/ # Random files .hypothesis/unicode_data/11.0.0/charmap.json.gz + +# Data files +benchmarks/data/*.h5 + +# ASV +.asv/ +pkgs/ +release_hashes.txt diff --git a/.mailmap b/.mailmap index a96f4a29f90..3aa26341c9c 100644 --- a/.mailmap +++ b/.mailmap @@ -270,3 +270,7 @@ Kim Lingemann kimsina Kim Lingemann kim Sumit Gupta + +Israel Roldan Israel Roldan +Israel Roldan AirvZxf +Israel Roldan airv_zxf diff --git a/asv.conf.json b/asv.conf.json index ab19c67ff26..351c2a78502 100644 --- a/asv.conf.json +++ b/asv.conf.json @@ -11,7 +11,7 @@ ], "branches": ["master"], "environment_type": "mamba", - "show_commit_url": "https://github.com/tardis-sn/tardis/commit", + "show_commit_url": "https://github.com/tardis-sn/tardis/commit/", "conda_environment_file": "tardis_env3.yml", "benchmark_dir": "benchmarks", "env_dir": ".asv/env", diff --git a/benchmarks/asv_by_release.bash b/benchmarks/asv_by_release.bash new file mode 100755 index 00000000000..815f1ebb636 --- /dev/null +++ b/benchmarks/asv_by_release.bash @@ -0,0 +1,99 @@ +#!/usr/bin/env bash + +RELEASE_LIST=$(git tag -l "release-202[34]*" | sort -r) + +readarray -t RELEASE_TAGS <<<"${RELEASE_LIST[@]}" +RELEASE_HASHES=() +for release_tag in "${RELEASE_TAGS[@]}"; do + echo "Tag: ${release_tag}" + HASH_COMMIT=$(git show-ref -s "${release_tag}") + RELEASE_HASHES+=("${HASH_COMMIT}") +done +echo "RELEASE_HASHES: ${#RELEASE_HASHES[*]}" + +ASV_CONFIG_PATH="/app/asv" +cd "${ASV_CONFIG_PATH}" || exit + +rm -f release_hashes.txt +touch release_hashes.txt +for release_hash in "${RELEASE_HASHES[@]}"; do + echo "${release_hash}" >>release_hashes.txt +done + +function show_timed_time { + local time=${1} + local milliseconds="${time: -3}" + local seconds=$((time / 1000)) + local minutes=0 + local minutes_display="" + local hours=0 + local hours_display="" + local days=0 + local days_display="" + + if [[ "${seconds}" -gt 59 ]]; then + minutes=$((seconds / 60)) + seconds=$((seconds % 60)) + minutes_display="${minutes}m " + fi + + if [[ "${minutes}" -gt 59 ]]; then + hours=$((minutes / 60)) + minutes=$((minutes % 60)) + minutes_display="${minutes}m " + hours_display="${hours}h " + fi + + if [[ "${hours}" -gt 23 ]]; then + days=$((hours / 24)) + hours=$((hours % 24)) + hours_display="${hours}h " + days_display="${days}d " + fi + + echo "${days_display}${hours_display}${minutes_display}${seconds}.${milliseconds}s" +} + +start=$(date +%s%N | cut -b1-13) + +# ASV has an argument called “bench”, which filters the benchmarks. +# I had two problems with ASV regarding benchmark filtering. +# 1. If we want to run only the `time_read_stella_model` benchmark. +# ASV will run 3 benchmarks using this argument +# `--bench time_read_stella_model`: +# - `time_read_stella_model`. +# - `time_read_stella_model_meta`. +# - `time_read_stella_model_data`. +# It is because the ASV uses regular expressions to match the name. +# To run only the benchmark `time_read_stella_model` we need to run: +# - `--bench "time_read_stella_model$"` +# - The `$` means that it matches the end of a string without +# consuming any characters. +# Note: +# - If the benchmark doesn't have parameters (`@parameterize`), +# then the name in ASV is the benchmark name without parameters. +# 2. The second problem is when I want to search for some benchmark that +# has parameters (`@parameterize`) because the name of the benchmark +# includes parenthesis and the parameters and their values. +# One example is `time_get_inverse_doppler_factor` which has 3 benchmarks +# with this prefix, and has parameters. +# To prevent it, we need to add the argument with this syntax: +# `--bench time_benchmark_name([^_A-Za-z]|$)` +# The regular expression match with all the parameters generated by ASV. +#time asv run \ +# --bench "time_read_stella_model$" +# --bench "time_get_inverse_doppler_factor([^_A-Za-z]|$)" \ +# --bench "time_get_inverse_doppler_factor_full_relativity([^_A-Za-z]|$)" \ +# release-2023.01.11..master + +# This command runs all benchmarks for all commits that have not been run. +time asv run \ + --skip-existing-commits \ + ALL + +end=$(date +%s%N | cut -b1-13) +runtime=$((end - start)) +display_time="$(show_timed_time ${runtime})" +echo "" +echo "Time: ${display_time}" +echo "" diff --git a/benchmarks/benchmark_base.py b/benchmarks/benchmark_base.py new file mode 100644 index 00000000000..012d33935bb --- /dev/null +++ b/benchmarks/benchmark_base.py @@ -0,0 +1,329 @@ +import re +from copy import deepcopy +from os.path import dirname, realpath, join +from pathlib import Path +from tempfile import mkstemp + +import astropy.units as u +import numpy as np +import pandas as pd +from numba import njit + +from benchmarks.util.nlte import NLTE +from tardis.io.atom_data import AtomData +from tardis.io.configuration import config_reader +from tardis.io.configuration.config_reader import Configuration +from tardis.io.util import yaml_load_file, YAMLLoader, HDFWriterMixin +from tardis.model import SimulationState +from tardis.montecarlo import NumbaModel, opacity_state_initialize +from tardis.montecarlo.montecarlo_numba import RPacket +from tardis.montecarlo.montecarlo_numba.packet_collections import ( + VPacketCollection, +) +from tardis.simulation import Simulation +from tardis.tests.fixtures.atom_data import DEFAULT_ATOM_DATA_UUID +from tardis.tests.fixtures.regression_data import RegressionData + + +class BenchmarkBase: + # It allows 10 minutes of runtime for each benchmark and includes + # the total time for all the repetitions for each benchmark. + timeout = 600 + + def __init__(self): + self.nlte = NLTE() + + @staticmethod + def get_relative_path(partial_path: str): + path = dirname(realpath(__file__)) + targets = Path(partial_path).parts + + for target in targets: + path = join(path, target) + + return path + + def get_absolute_path(self, partial_path): + partial_path = "../" + partial_path + + return self.get_relative_path(partial_path) + + @property + def tardis_config_verysimple(self): + filename = self.get_absolute_path( + "tardis/io/configuration/tests/data/tardis_configv1_verysimple.yml" + ) + return yaml_load_file( + filename, + YAMLLoader, + ) + + @property + def tardis_ref_path(self): + # TODO: This route is fixed but needs to get from the arguments given in the command line. + # /app/tardis-refdata + return "/app/tardis-refdata" + + @property + def atomic_dataset(self) -> AtomData: + atomic_data = AtomData.from_hdf(self.atomic_data_fname) + + if atomic_data.md5 != DEFAULT_ATOM_DATA_UUID: + message = f'Need default Kurucz atomic dataset (md5="{DEFAULT_ATOM_DATA_UUID}")' + raise Exception(message) + else: + return atomic_data + + @property + def atomic_data_fname(self): + atomic_data_fname = ( + f"{self.tardis_ref_path}/atom_data/kurucz_cd23_chianti_H_He.h5" + ) + + if not Path(atomic_data_fname).exists(): + atom_data_missing_str = ( + f"{atomic_data_fname} atomic datafiles " + f"does not seem to exist" + ) + raise Exception(atom_data_missing_str) + + return atomic_data_fname + + @property + def example_configuration_dir(self): + return self.get_absolute_path("tardis/io/configuration/tests/data") + + @property + def hdf_file_path(self): + # TODO: Delete this files after ASV runs the benchmarks. + # ASV create a temporal directory in runtime per test: `tmpiuxngvlv`. + # The ASV and ASV_Runner, not has some way to get this temporal directory. + # The idea is use this temporal folders to storage this created temporal file. + _, path = mkstemp("-tardis-benchmark-hdf_buffer-test.hdf") + return path + + def create_temporal_file(self, suffix=None): + # TODO: Delete this files after ASV runs the benchmarks. + # ASV create a temporal directory in runtime per test: `tmpiuxngvlv`. + # The ASV and ASV_Runner, not has some way to get this temporal directory. + # The idea is use this temporal folders to storage this created temporal file. + suffix_str = "" if suffix is None else f"-{suffix}" + _, path = mkstemp(suffix_str) + return path + + @property + def gamma_ray_simulation_state(self): + self.gamma_ray_config.model.structure.velocity.start = 1.0 * u.km / u.s + self.gamma_ray_config.model.structure.density.rho_0 = ( + 5.0e2 * u.g / u.cm**3 + ) + self.gamma_ray_config.supernova.time_explosion = 150 * u.d + + return SimulationState.from_config( + self.gamma_ray_config, atom_data=self.atomic_dataset + ) + + @property + def gamma_ray_config(self): + yml_path = f"{self.example_configuration_dir}/tardis_configv1_density_exponential_nebular_multi_isotope.yml" + + return config_reader.Configuration.from_yaml(yml_path) + + @property + def example_model_file_dir(self): + return self.get_absolute_path("tardis/io/model/readers/tests/data") + + @property + def kurucz_atomic_data(self) -> AtomData: + return deepcopy(self.atomic_dataset) + + @property + def example_csvy_file_dir(self): + return self.get_absolute_path("tardis/model/tests/data/") + + @property + def simulation_verysimple(self): + atomic_data = deepcopy(self.atomic_dataset) + sim = Simulation.from_config( + self.config_verysimple, atom_data=atomic_data + ) + sim.iterate(4000) + return sim + + @property + def config_verysimple(self): + return Configuration.from_yaml( + f"{self.example_configuration_dir}/tardis_configv1_verysimple.yml" + ) + + class CustomPyTestRequest: + def __init__( + self, + tardis_regression_data_path: str, + node_name: str, + node_module_name: str, + regression_data_dir: str, + ): + self.tardis_regression_data_path = tardis_regression_data_path + self.node_name = node_name + self.node_module_name = node_module_name + self.regression_data_dir = regression_data_dir + + @property + def config(self): + class SubClass: + @staticmethod + def getoption(option): + if option == "--tardis-regression-data": + return self.tardis_regression_data_path + return None + + return SubClass() + + @property + def node(self): + class SubClass: + def __init__(self, parent): + self.parent = parent + + @property + def name(self): + return self.parent.node_name + + @property + def module(self): + class SubSubClass: + def __init__(self, parent): + self.parent = parent + + @property + def __name__(self): + return self.parent.node_module_name + + return SubSubClass(self.parent) + + return SubClass(self) + + @property + def cls(self): + return None + + @property + def relative_regression_data_dir(self): + return self.regression_data_dir + + @staticmethod + def regression_data(request: CustomPyTestRequest): + return RegressionData(request) + + @property + def packet(self): + return RPacket( + r=7.5e14, + nu=self.verysimple_packet_collection.initial_nus[0], + mu=self.verysimple_packet_collection.initial_mus[0], + energy=self.verysimple_packet_collection.initial_energies[0], + seed=1963, + index=0, + ) + + @property + def verysimple_packet_collection(self): + return ( + self.nb_simulation_verysimple.transport.transport_state.packet_collection + ) + + @property + def nb_simulation_verysimple(self): + atomic_data = deepcopy(self.atomic_dataset) + sim = Simulation.from_config( + self.config_verysimple, atom_data=atomic_data + ) + sim.iterate(10) + return sim + + @property + def verysimple_numba_model(self): + model = self.nb_simulation_verysimple.simulation_state + return NumbaModel( + model.time_explosion.to("s").value, + ) + + @property + def verysimple_opacity_state(self): + return opacity_state_initialize( + self.nb_simulation_verysimple.plasma, + line_interaction_type="macroatom", + ) + + @property + def static_packet(self): + return RPacket( + r=7.5e14, + nu=0.4, + mu=0.3, + energy=0.9, + seed=1963, + index=0, + ) + + @property + def set_seed_fixture(self): + def set_seed(value): + np.random.seed(value) + + return njit(set_seed) + + @property + def verysimple_3vpacket_collection(self): + spectrum_frequency = ( + self.nb_simulation_verysimple.transport.spectrum_frequency.value + ) + return VPacketCollection( + source_rpacket_index=0, + spectrum_frequency=spectrum_frequency, + number_of_vpackets=3, + v_packet_spawn_start_frequency=0, + v_packet_spawn_end_frequency=np.inf, + temporary_v_packet_bins=0, + ) + + @property + def verysimple_numba_radial_1d_geometry(self): + return ( + self.nb_simulation_verysimple.simulation_state.geometry.to_numba() + ) + + @property + def simulation_verysimple_vpacket_tracking(self): + atomic_data = deepcopy(self.atomic_dataset) + sim = Simulation.from_config( + self.config_verysimple, + atom_data=atomic_data, + virtual_packet_logging=True, + ) + sim.last_no_of_packets = 4000 + sim.run_final() + return sim + + @property + def generate_reference(self): + # TODO: Investigate how to get the `--generate-reference` parameter passed in the command line. + # `request.config.getoption("--generate-reference")` + option = None + if option is None: + return False + else: + return option + + @property + def tardis_ref_data(self): + # TODO: This function is not working in the benchmarks. + if self.generate_reference: + mode = "w" + else: + mode = "r" + with pd.HDFStore( + f"{self.tardis_ref_path}/unit_test_data.h5", mode=mode + ) as store: + yield store diff --git a/benchmarks/benchmark_run_tardis.py b/benchmarks/benchmark_run_tardis.py deleted file mode 100644 index 7a62dc5822a..00000000000 --- a/benchmarks/benchmark_run_tardis.py +++ /dev/null @@ -1,20 +0,0 @@ -"""Basic TARDIS Benchmark.""" -import os -from tardis.io.configuration.config_reader import Configuration -from tardis import run_tardis - -class Benchmarkruntardis: - """Class to benchmark the run_tardis function. - """ - timeout = 200 - - def setup(self): - filename = "tardis_configv1_benchmark.yml" - dir_path = os.path.dirname(os.path.realpath(__file__)) - path = os.path.join(dir_path, "data", filename) - config = Configuration.from_yaml(path) - config.atom_data = "kurucz_cd23_chianti_H_He.h5" - self.config = config - - def time_run_tardis(self): - sim = run_tardis(self.config, log_level="ERROR", show_progress_bars=False) diff --git a/benchmarks/data/tardis_configv1_benchmark.yml b/benchmarks/data/tardis_configv1_benchmark.yml index 22e1c9ac30b..578cd76b1f7 100644 --- a/benchmarks/data/tardis_configv1_benchmark.yml +++ b/benchmarks/data/tardis_configv1_benchmark.yml @@ -4,7 +4,7 @@ supernova: luminosity_requested: 2.8e9 solLum time_explosion: 13 day -atom_data: kurucz_atom_pure_simple.h5 +atom_data: kurucz_cd23_chianti_H_He.h5 model: structure: diff --git a/benchmarks/montecarlo_montecarlo_numba_interaction.py b/benchmarks/montecarlo_montecarlo_numba_interaction.py new file mode 100644 index 00000000000..fa7a2f34552 --- /dev/null +++ b/benchmarks/montecarlo_montecarlo_numba_interaction.py @@ -0,0 +1,103 @@ +""" +Basic TARDIS Benchmark. +""" + +import numpy as np +from asv_runner.benchmarks.mark import parameterize, skip_benchmark + +import tardis.montecarlo.montecarlo_numba.interaction as interaction +from benchmarks.benchmark_base import BenchmarkBase +from tardis.montecarlo.montecarlo_numba.numba_interface import ( + LineInteractionType, +) + + +@skip_benchmark +class BenchmarkMontecarloMontecarloNumbaInteraction(BenchmarkBase): + """ + Class to benchmark the numba interaction function. + """ + + def time_thomson_scatter(self): + packet = self.packet + init_mu = packet.mu + init_nu = packet.nu + init_energy = packet.energy + time_explosion = self.verysimple_numba_model.time_explosion + + interaction.thomson_scatter(packet, time_explosion) + + assert np.abs(packet.mu - init_mu) > 1e-7 + assert np.abs(packet.nu - init_nu) > 1e-7 + assert np.abs(packet.energy - init_energy) > 1e-7 + + @parameterize( + { + "Line interaction type": [ + LineInteractionType.SCATTER, + LineInteractionType.DOWNBRANCH, + LineInteractionType.MACROATOM, + ], + } + ) + def time_line_scatter(self, line_interaction_type): + packet = self.packet + init_mu = packet.mu + init_nu = packet.nu + init_energy = packet.energy + packet.initialize_line_id( + self.verysimple_opacity_state, self.verysimple_numba_model + ) + time_explosion = self.verysimple_numba_model.time_explosion + + interaction.line_scatter( + packet, + time_explosion, + line_interaction_type, + self.verysimple_opacity_state, + ) + + assert np.abs(packet.mu - init_mu) > 1e-7 + assert np.abs(packet.nu - init_nu) > 1e-7 + assert np.abs(packet.energy - init_energy) > 1e-7 + + @parameterize( + { + "Test packet": [ + { + "mu": 0.8599443103322428, + "emission_line_id": 1000, + "energy": 0.9114437898710559, + }, + { + "mu": -0.6975116557422458, + "emission_line_id": 2000, + "energy": 0.8803098648913266, + }, + { + "mu": -0.7115661419975774, + "emission_line_id": 0, + "energy": 0.8800385929341252, + }, + ] + } + ) + def time_line_emission(self, test_packet): + emission_line_id = test_packet["emission_line_id"] + packet = self.packet + packet.mu = test_packet["mu"] + packet.energy = test_packet["energy"] + packet.initialize_line_id( + self.verysimple_opacity_state, self.verysimple_numba_model + ) + + time_explosion = self.verysimple_numba_model.time_explosion + + interaction.line_emission( + packet, + emission_line_id, + time_explosion, + self.verysimple_opacity_state, + ) + + assert packet.next_line_id == emission_line_id + 1 diff --git a/benchmarks/montecarlo_montecarlo_numba_numba_formal_integral_p.py b/benchmarks/montecarlo_montecarlo_numba_numba_formal_integral_p.py new file mode 100644 index 00000000000..0f2632f34d6 --- /dev/null +++ b/benchmarks/montecarlo_montecarlo_numba_numba_formal_integral_p.py @@ -0,0 +1,178 @@ +""" +Basic TARDIS Benchmark. +""" + +import numpy as np +from asv_runner.benchmarks.mark import parameterize, skip_benchmark + +import tardis.montecarlo.montecarlo_numba.formal_integral as formal_integral +from benchmarks.benchmark_base import BenchmarkBase +from tardis import constants as c +from tardis.model.geometry.radial1d import NumbaRadial1DGeometry +from tardis.montecarlo.montecarlo_numba.numba_interface import NumbaModel +from tardis.util.base import intensity_black_body + + +class BenchmarkMontecarloMontecarloNumbaNumbaFormalIntegral(BenchmarkBase): + """ + Class to benchmark the numba formal integral function. + """ + + @parameterize( + { + "nu": [1e14, 0, 1], + "temperature": [1e4, 1, 1], + } + ) + def time_intensity_black_body(self, nu, temperature): + func = formal_integral.intensity_black_body + actual = func(nu, temperature) + print(actual, type(actual)) + intensity_black_body(nu, temperature) + + @parameterize({"N": (1e2, 1e3, 1e4, 1e5)}) + def time_trapezoid_integration(self, n): + func = formal_integral.trapezoid_integration + h = 1.0 + n = int(n) + data = np.random.random(n) + + func(data, h) + np.trapz(data) + + @staticmethod + def calculate_z(r, p): + return np.sqrt(r * r - p * p) + + TESTDATA = [ + np.linspace(1, 2, 3, dtype=np.float64), + np.linspace(0, 1, 3), + # np.linspace(1, 2, 10, dtype=np.float64), + ] + + def formal_integral_geometry(self, r): + # NOTE: PyTest is generating a full matrix with all the permutations. + # For the `time_calculate_z` function with values: [0.0, 0.5, 1.0] + # - p=0.0, formal_integral_geometry0-0.0, param["r"]: [1. 1.5 2. ] + # - p=0.5, formal_integral_geometry0-0.5, param["r"]: [1. 1.5 2. ] + # - p=1.0, formal_integral_geometry0-1.0, param["r"]: [1. 1.5 2. ] + # - p=0.0, formal_integral_geometry1-0.0, param["r"]: [0. 0.5 1. ] + # - p=1.0, formal_integral_geometry1-1.0, param["r"]: [0. 0.5 1. ] + # Same for `test_populate_z_photosphere` function + # And for `test_populate_z_shells` function + # - p=1e-05, formal_integral_geometry0-1e-05, param["r"]: [1. 1.5 2. ] + # - p=0.5, formal_integral_geometry0-0.5, param["r"]: [1. 1.5 2. ] + # - p=0.99, formal_integral_geometry0-0.99, param["r"]: [1. 1.5 2. ] + # - p=1, formal_integral_geometry0-1, param["r"]: [1. 1.5 2. ] + # - p=1e-05, formal_integral_geometry1-1e-05, param["r"]: [0. 0.5 1. ] + # - p=0.5, formal_integral_geometry1-0.5, param["r"]: [0. 0.5 1. ] + # - p=0.99, formal_integral_geometry1-0.99, param["r"]: [0. 0.5 1. ] + # - p=1, formal_integral_geometry1-1, param["r"]: [0. 0.5 1. ] + geometry = NumbaRadial1DGeometry( + r[:-1], + r[1:], + r[:-1] * c.c.cgs.value, + r[1:] * c.c.cgs.value, + ) + return geometry + + @property + def formal_integral_model(self): + model = NumbaModel( + 1 / c.c.cgs.value, + ) + return model + + @parameterize({"p": [0.0, 0.5, 1.0], "Test data": TESTDATA}) + def time_calculate_z(self, p, test_data): + func = formal_integral.calculate_z + inv_t = 1.0 / self.formal_integral_model.time_explosion + len(self.formal_integral_geometry(test_data).r_outer) + r_outer = self.formal_integral_geometry(test_data).r_outer + + for r in r_outer: + actual = func(r, p, inv_t) + if p >= r: + assert actual == 0 + else: + np.sqrt(r * r - p * p) * formal_integral.C_INV * inv_t + + @skip_benchmark + @parameterize({"p": [0, 0.5, 1], "Test data": TESTDATA}) + def time_populate_z_photosphere(self, p, test_data): + formal_integral.FormalIntegrator( + self.formal_integral_geometry(test_data), None, None + ) + func = formal_integral.populate_z + size = len(self.formal_integral_geometry(test_data).r_outer) + r_inner = self.formal_integral_geometry(test_data).r_inner + self.formal_integral_geometry(test_data).r_outer + + p = r_inner[0] * p + oz = np.zeros_like(r_inner) + oshell_id = np.zeros_like(oz, dtype=np.int64) + + n = func( + self.formal_integral_geometry(test_data), + self.formal_integral_geometry(test_data), + p, + oz, + oshell_id, + ) + assert n == size + + @skip_benchmark + @parameterize({"p": [1e-5, 0.5, 0.99, 1], "Test data": TESTDATA}) + def time_populate_z_shells(self, p, test_data): + formal_integral.FormalIntegrator( + self.formal_integral_geometry(test_data), None, None + ) + func = formal_integral.populate_z + + size = len(self.formal_integral_geometry(test_data).r_inner) + r_inner = self.formal_integral_geometry(test_data).r_inner + r_outer = self.formal_integral_geometry(test_data).r_outer + + p = r_inner[0] + (r_outer[-1] - r_inner[0]) * p + idx = np.searchsorted(r_outer, p, side="right") + + oz = np.zeros(size * 2) + oshell_id = np.zeros_like(oz, dtype=np.int64) + + offset = size - idx + + expected_n = (offset) * 2 + expected_oz = np.zeros_like(oz) + expected_oshell_id = np.zeros_like(oshell_id) + + # Calculated way to determine which shells get hit + expected_oshell_id[:expected_n] = ( + np.abs(np.arange(0.5, expected_n, 1) - offset) - 0.5 + idx + ) + + expected_oz[0:offset] = 1 + self.calculate_z( + r_outer[np.arange(size, idx, -1) - 1], p + ) + expected_oz[offset:expected_n] = 1 - self.calculate_z( + r_outer[np.arange(idx, size, 1)], p + ) + + n = func( + self.formal_integral_geometry(test_data), + self.formal_integral_geometry(test_data), + p, + oz, + oshell_id, + ) + + assert n == expected_n + + @parameterize({"N": [100, 1000, 10000]}) + def time_calculate_p_values(self, n): + r = 1.0 + func = formal_integral.calculate_p_values + + expected = r / (n - 1) * np.arange(0, n, dtype=np.float64) + np.zeros_like(expected, dtype=np.float64) + + func(r, n) diff --git a/benchmarks/montecarlo_montecarlo_numba_numba_interface.py b/benchmarks/montecarlo_montecarlo_numba_numba_interface.py new file mode 100644 index 00000000000..c921184a160 --- /dev/null +++ b/benchmarks/montecarlo_montecarlo_numba_numba_interface.py @@ -0,0 +1,77 @@ +""" +Basic TARDIS Benchmark. +""" + +import numpy as np +from asv_runner.benchmarks.mark import parameterize + +import tardis.montecarlo.montecarlo_numba.numba_interface as numba_interface +from benchmarks.benchmark_base import BenchmarkBase + + +class BenchmarkMontecarloMontecarloNumbaNumbaInterface(BenchmarkBase): + """ + Class to benchmark the numba interface function. + """ + + @parameterize({"Input params": ["scatter", "macroatom", "downbranch"]}) + def time_opacity_state_initialize(self, input_params): + line_interaction_type = input_params + plasma = self.nb_simulation_verysimple.plasma + numba_interface.opacity_state_initialize(plasma, line_interaction_type) + + if line_interaction_type == "scatter": + np.zeros(1, dtype=np.int64) + + def time_VPacketCollection_add_packet(self): + verysimple_3vpacket_collection = self.verysimple_3vpacket_collection + assert verysimple_3vpacket_collection.length == 0 + + nus = [3.0e15, 0.0, 1e15, 1e5] + energies = [0.4, 0.1, 0.6, 1e10] + initial_mus = [0.1, 0, 1, 0.9] + initial_rs = [3e42, 4.5e45, 0, 9.0e40] + last_interaction_in_nus = np.array( + [3.0e15, 0.0, 1e15, 1e5], dtype=np.float64 + ) + last_interaction_types = np.array([1, 1, 3, 2], dtype=np.int64) + last_interaction_in_ids = np.array([100, 0, 1, 1000], dtype=np.int64) + last_interaction_out_ids = np.array( + [1201, 123, 545, 1232], dtype=np.int64 + ) + last_interaction_shell_ids = np.array([2, -1, 6, 0], dtype=np.int64) + + for ( + nu, + energy, + initial_mu, + initial_r, + last_interaction_in_nu, + last_interaction_type, + last_interaction_in_id, + last_interaction_out_id, + last_interaction_shell_id, + ) in zip( + nus, + energies, + initial_mus, + initial_rs, + last_interaction_in_nus, + last_interaction_types, + last_interaction_in_ids, + last_interaction_out_ids, + last_interaction_shell_ids, + ): + verysimple_3vpacket_collection.add_packet( + nu, + energy, + initial_mu, + initial_r, + last_interaction_in_nu, + last_interaction_type, + last_interaction_in_id, + last_interaction_out_id, + last_interaction_shell_id, + ) + + assert verysimple_3vpacket_collection.length == 9 diff --git a/benchmarks/montecarlo_montecarlo_numba_opacities.py b/benchmarks/montecarlo_montecarlo_numba_opacities.py new file mode 100644 index 00000000000..c10092a129b --- /dev/null +++ b/benchmarks/montecarlo_montecarlo_numba_opacities.py @@ -0,0 +1,91 @@ +""" +Basic TARDIS Benchmark. +""" + +from asv_runner.benchmarks.mark import parameterize + +import tardis.montecarlo.montecarlo_numba.opacities as calculate_opacity +from benchmarks.benchmark_base import BenchmarkBase + + +class BenchmarkMontecarloMontecarloNumbaOpacities(BenchmarkBase): + """ + Class to benchmark the numba opacities function. + """ + + @parameterize( + { + "Electron number density": [ + 1.0e11, + 1e15, + 1e5, + ], + "Energy": [ + 511.0, + 255.5, + 511.0e7, + ], + } + ) + def time_compton_opacity_calculation(self, electron_number_density, energy): + calculate_opacity.compton_opacity_calculation( + energy, electron_number_density + ) + + @parameterize( + { + "Ejecta density": [ + 1.0, + 1e-2, + 1e-2, + 1e5, + ], + "Energy": [ + 511.0, + 255.5, + 255.5, + 511.0e7, + ], + "Iron group fraction": [ + 0.0, + 0.5, + 0.25, + 1.0, + ], + } + ) + def time_photoabsorption_opacity_calculation( + self, ejecta_density, energy, iron_group_fraction + ): + calculate_opacity.photoabsorption_opacity_calculation( + energy, ejecta_density, iron_group_fraction + ) + + @parameterize( + { + "Ejecta density": [ + 1.0, + 1e-2, + 1e-2, + 1e5, + ], + "Energy": [ + 511.0, + 1500, + 1200, + 511.0e7, + ], + "Iron group fraction": [ + 0.0, + 0.5, + 0.25, + 1.0, + ], + } + ) + def time_pair_creation_opacity_calculation( + self, ejecta_density, energy, iron_group_fraction + ): + calculate_opacity.pair_creation_opacity_calculation( + energy, ejecta_density, iron_group_fraction + ) diff --git a/benchmarks/montecarlo_montecarlo_numba_packet.py b/benchmarks/montecarlo_montecarlo_numba_packet.py new file mode 100644 index 00000000000..c16077596ac --- /dev/null +++ b/benchmarks/montecarlo_montecarlo_numba_packet.py @@ -0,0 +1,307 @@ +""" +Basic TARDIS Benchmark. +""" + +import numpy as np +from asv_runner.benchmarks.mark import parameterize, skip_benchmark + +import tardis.montecarlo.estimators.radfield_mc_estimators +import tardis.montecarlo.estimators.radfield_mc_estimators +import tardis.montecarlo.montecarlo_numba.numba_interface as numba_interface +import tardis.montecarlo.montecarlo_numba.opacities as opacities +import tardis.montecarlo.montecarlo_numba.r_packet as r_packet +import tardis.montecarlo.montecarlo_numba.utils as utils +import tardis.transport.frame_transformations as frame_transformations +import tardis.transport.geometry.calculate_distances as calculate_distances +import tardis.transport.r_packet_transport as r_packet_transport +from benchmarks.benchmark_base import BenchmarkBase +from tardis.model.geometry.radial1d import NumbaRadial1DGeometry +from tardis.montecarlo.estimators.radfield_mc_estimators import ( + update_line_estimators, +) + + +class BenchmarkMontecarloMontecarloNumbaPacket(BenchmarkBase): + """ + Class to benchmark the numba packet function. + """ + + @property + def geometry(self): + return NumbaRadial1DGeometry( + r_inner=np.array([6.912e14, 8.64e14], dtype=np.float64), + r_outer=np.array([8.64e14, 1.0368e15], dtype=np.float64), + v_inner=np.array([-1, -1], dtype=np.float64), + v_outer=np.array([-1, -1], dtype=np.float64), + ) + + @property + def model(self): + return numba_interface.NumbaModel( + time_explosion=5.2e7, + ) + + @property + def estimators(self): + return tardis.montecarlo.estimators.radfield_mc_estimators.RadiationFieldMCEstimators( + j_estimator=np.array([0.0, 0.0], dtype=np.float64), + nu_bar_estimator=np.array([0.0, 0.0], dtype=np.float64), + j_blue_estimator=np.array( + [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], dtype=np.float64 + ), + Edotlu_estimator=np.array( + [[0.0, 0.0, 1.0], [0.0, 0.0, 1.0]], dtype=np.float64 + ), + photo_ion_estimator=np.empty((0, 0), dtype=np.float64), + stim_recomb_estimator=np.empty((0, 0), dtype=np.float64), + bf_heating_estimator=np.empty((0, 0), dtype=np.float64), + stim_recomb_cooling_estimator=np.empty((0, 0), dtype=np.float64), + photo_ion_estimator_statistics=np.empty((0, 0), dtype=np.int64), + ) + + @parameterize( + { + "Packet params": [ + {"mu": 0.3, "r": 7.5e14}, + {"mu": -0.3, "r": 7.5e13}, + {"mu": -0.3, "r": 7.5e14}, + ] + } + ) + def time_calculate_distance_boundary(self, packet_params): + mu = packet_params["mu"] + r = packet_params["r"] + + calculate_distances.calculate_distance_boundary( + r, mu, self.geometry.r_inner[0], self.geometry.r_outer[0] + ) + + @parameterize( + { + "Parameters": [ + { + "packet": {"nu_line": 0.1, "is_last_line": True}, + "expected": None, + }, + { + "packet": {"nu_line": 0.2, "is_last_line": False}, + "expected": None, + }, + { + "packet": {"nu_line": 0.5, "is_last_line": False}, + "expected": utils.MonteCarloException, + }, + { + "packet": {"nu_line": 0.6, "is_last_line": False}, + "expected": utils.MonteCarloException, + }, + ] + } + ) + def time_calculate_distance_line(self, parameters): + packet_params = parameters["packet"] + expected_params = parameters["expected"] + nu_line = packet_params["nu_line"] + is_last_line = packet_params["is_last_line"] + + time_explosion = self.model.time_explosion + + doppler_factor = frame_transformations.get_doppler_factor( + self.static_packet.r, self.static_packet.mu, time_explosion + ) + comov_nu = self.static_packet.nu * doppler_factor + + obtained_tardis_error = None + try: + calculate_distances.calculate_distance_line( + self.static_packet, + comov_nu, + is_last_line, + nu_line, + time_explosion, + ) + except utils.MonteCarloException: + obtained_tardis_error = utils.MonteCarloException + + assert obtained_tardis_error == expected_params + + @parameterize( + { + "Parameters": [ + { + "electron_density": 1e-5, + "tua_event": 1e10, + }, + {"electron_density": 1.0, "tua_event": 1e10}, + ] + } + ) + def time_calculate_distance_electron(self, parameters): + electron_density = parameters["electron_density"] + tau_event = parameters["tua_event"] + calculate_distances.calculate_distance_electron( + electron_density, tau_event + ) + + @parameterize( + { + "Parameters": [ + { + "electron_density": 1e-5, + "distance": 1.0, + }, + { + "electron_density": 1e10, + "distance": 1e10, + }, + { + "electron_density": -1, + "distance": 0, + }, + { + "electron_density": -1e10, + "distance": -1e10, + }, + ] + } + ) + def time_calculate_tau_electron(self, parameters): + electron_density = parameters["electron_density"] + distance = parameters["distance"] + opacities.calculate_tau_electron(electron_density, distance) + + def time_get_random_mu(self): + self.set_seed_fixture(1963) + + output1 = utils.get_random_mu() + assert output1 == 0.9136407866175174 + + @parameterize( + { + "Parameters": [ + { + "cur_line_id": 0, + "distance_trace": 1e12, + "time_explosion": 5.2e7, + }, + { + "cur_line_id": 0, + "distance_trace": 0, + "time_explosion": 5.2e7, + }, + { + "cur_line_id": 1, + "distance_trace": 1e5, + "time_explosion": 1e10, + }, + ] + } + ) + def time_update_line_estimators(self, parameters): + cur_line_id = parameters["cur_line_id"] + distance_trace = parameters["distance_trace"] + time_explosion = parameters["time_explosion"] + update_line_estimators( + self.estimators, + self.static_packet, + cur_line_id, + distance_trace, + time_explosion, + ) + + @parameterize( + { + "Parameters": [ + { + "current_shell_id": 132, + "delta_shell": 11, + "no_of_shells": 132, + }, + { + "current_shell_id": 132, + "delta_shell": 1, + "no_of_shells": 133, + }, + { + "current_shell_id": 132, + "delta_shell": 2, + "no_of_shells": 133, + }, + ] + } + ) + def time_move_packet_across_shell_boundary_emitted(self, parameters): + current_shell_id = parameters["current_shell_id"] + delta_shell = parameters["delta_shell"] + no_of_shells = parameters["no_of_shells"] + packet = self.packet + packet.current_shell_id = current_shell_id + r_packet_transport.move_packet_across_shell_boundary( + packet, delta_shell, no_of_shells + ) + assert packet.status == r_packet.PacketStatus.EMITTED + + @skip_benchmark + @parameterize( + { + "Parameters": [ + { + "current_shell_id": 132, + "delta_shell": 132, + "no_of_shells": 132, + }, + { + "current_shell_id": -133, + "delta_shell": -133, + "no_of_shells": -1e9, + }, + { + "current_shell_id": 132, + "delta_shell": 133, + "no_of_shells": 133, + }, + ] + } + ) + def time_move_packet_across_shell_boundary_reabsorbed(self, parameters): + current_shell_id = parameters["current_shell_id"] + delta_shell = parameters["delta_shell"] + no_of_shells = parameters["no_of_shells"] + packet = self.packet + packet.current_shell_id = current_shell_id + r_packet_transport.move_packet_across_shell_boundary( + packet, delta_shell, no_of_shells + ) + assert packet.status == r_packet.PacketStatus.REABSORBED + + @parameterize( + { + "Parameters": [ + { + "current_shell_id": 132, + "delta_shell": -1, + "no_of_shells": 199, + }, + { + "current_shell_id": 132, + "delta_shell": 0, + "no_of_shells": 132, + }, + { + "current_shell_id": 132, + "delta_shell": 20, + "no_of_shells": 154, + }, + ] + } + ) + def time_move_packet_across_shell_boundary_increment(self, parameters): + current_shell_id = parameters["current_shell_id"] + delta_shell = parameters["delta_shell"] + no_of_shells = parameters["no_of_shells"] + packet = self.packet + packet.current_shell_id = current_shell_id + r_packet_transport.move_packet_across_shell_boundary( + packet, delta_shell, no_of_shells + ) + assert packet.current_shell_id == current_shell_id + delta_shell diff --git a/benchmarks/montecarlo_montecarlo_numba_r_packet.py b/benchmarks/montecarlo_montecarlo_numba_r_packet.py new file mode 100644 index 00000000000..c5f0128ad4e --- /dev/null +++ b/benchmarks/montecarlo_montecarlo_numba_r_packet.py @@ -0,0 +1,67 @@ +""" +Basic TARDIS Benchmark. +""" + +from copy import deepcopy + +from benchmarks.benchmark_base import BenchmarkBase +from tardis.base import run_tardis +from tardis.montecarlo.montecarlo_numba.r_packet import ( + rpacket_trackers_to_dataframe, +) + + +class BenchmarkMontecarloMontecarloNumbaRPacket(BenchmarkBase): + """ + Class to benchmark the numba R packet function. + """ + + @property + def simulation_rpacket_tracking_enabled(self): + config_verysimple = self.config_verysimple + config_verysimple.montecarlo.iterations = 3 + config_verysimple.montecarlo.no_of_packets = 4000 + config_verysimple.montecarlo.last_no_of_packets = -1 + config_verysimple.spectrum.virtual.virtual_packet_logging = True + config_verysimple.montecarlo.no_of_virtual_packets = 1 + config_verysimple.montecarlo.tracking.track_rpacket = True + config_verysimple.spectrum.num = 2000 + atomic_data = deepcopy(self.atomic_dataset) + sim = run_tardis( + config_verysimple, + atom_data=atomic_data, + show_convergence_plots=False, + ) + return sim + + def time_rpacket_trackers_to_dataframe(self): + sim = self.simulation_rpacket_tracking_enabled + transport_state = sim.transport.transport_state + rtracker_df = rpacket_trackers_to_dataframe( + transport_state.rpacket_tracker + ) + + # check df shape and column names + assert rtracker_df.shape == ( + sum( + [len(tracker.r) for tracker in transport_state.rpacket_tracker] + ), + 8, + ) + + # check all data with rpacket_tracker + expected_rtrackers = [] + for rpacket in transport_state.rpacket_tracker: + for rpacket_step_no in range(len(rpacket.r)): + expected_rtrackers.append( + [ + rpacket.status[rpacket_step_no], + rpacket.seed, + rpacket.r[rpacket_step_no], + rpacket.nu[rpacket_step_no], + rpacket.mu[rpacket_step_no], + rpacket.energy[rpacket_step_no], + rpacket.shell_id[rpacket_step_no], + rpacket.interaction_type[rpacket_step_no], + ] + ) diff --git a/benchmarks/montecarlo_montecarlo_numba_vpacket.py b/benchmarks/montecarlo_montecarlo_numba_vpacket.py new file mode 100644 index 00000000000..a1b711a2b22 --- /dev/null +++ b/benchmarks/montecarlo_montecarlo_numba_vpacket.py @@ -0,0 +1,119 @@ +""" +Basic TARDIS Benchmark. +""" + +import numpy as np + +import tardis.montecarlo.montecarlo_numba.vpacket as vpacket +from benchmarks.benchmark_base import BenchmarkBase +from tardis.transport.frame_transformations import ( + get_doppler_factor, +) + + +class BenchmarkMontecarloMontecarloNumbaVpacket(BenchmarkBase): + """ + Class to benchmark the single packet loop function. + """ + + @property + def v_packet(self): + return vpacket.VPacket( + r=7.5e14, + nu=4e15, + mu=0.3, + energy=0.9, + current_shell_id=0, + next_line_id=0, + index=0, + ) + + def v_packet_initialize_line_id(self, v_packet, opacity_state, numba_model): + inverse_line_list_nu = opacity_state.line_list_nu[::-1] + doppler_factor = get_doppler_factor( + v_packet.r, v_packet.mu, numba_model.time_explosion + ) + comov_nu = v_packet.nu * doppler_factor + next_line_id = len(opacity_state.line_list_nu) - np.searchsorted( + inverse_line_list_nu, comov_nu + ) + v_packet.next_line_id = next_line_id + + def time_trace_vpacket_within_shell(self): + v_packet = self.v_packet + verysimple_numba_radial_1d_geometry = ( + self.verysimple_numba_radial_1d_geometry + ) + verysimple_numba_model = self.verysimple_numba_model + verysimple_opacity_state = self.verysimple_opacity_state + + # Give the vpacket a reasonable line ID + self.v_packet_initialize_line_id( + v_packet, verysimple_opacity_state, verysimple_numba_model + ) + + ( + tau_trace_combined, + distance_boundary, + delta_shell, + ) = vpacket.trace_vpacket_within_shell( + v_packet, + verysimple_numba_radial_1d_geometry, + verysimple_numba_model, + verysimple_opacity_state, + ) + + assert delta_shell == 1 + + def time_trace_vpacket(self): + v_packet = self.v_packet + verysimple_numba_radial_1d_geometry = ( + self.verysimple_numba_radial_1d_geometry + ) + verysimple_numba_model = self.verysimple_numba_model + verysimple_opacity_state = self.verysimple_opacity_state + + # Set seed because of RNG in trace_vpacket + np.random.seed(1) + + # Give the vpacket a reasonable line ID + self.v_packet_initialize_line_id( + v_packet, verysimple_opacity_state, verysimple_numba_model + ) + + tau_trace_combined = vpacket.trace_vpacket( + v_packet, + verysimple_numba_radial_1d_geometry, + verysimple_numba_model, + verysimple_opacity_state, + ) + + assert v_packet.next_line_id == 2773 + assert v_packet.current_shell_id == 1 + + @property + def broken_packet(self): + return vpacket.VPacket( + r=1286064000000000.0, + nu=1660428912896553.2, + mu=0.4916053094346575, + energy=2.474533071386993e-07, + index=3, + current_shell_id=0, + next_line_id=5495, + ) + + def time_trace_bad_vpacket(self): + broken_packet = self.broken_packet + verysimple_numba_radial_1d_geometry = ( + self.verysimple_numba_radial_1d_geometry + ) + verysimple_numba_model = self.verysimple_numba_model + verysimple_opacity_state = self.verysimple_opacity_state + + vpacket.trace_vpacket( + broken_packet, + verysimple_numba_radial_1d_geometry, + verysimple_numba_model, + verysimple_opacity_state, + ) diff --git a/benchmarks/run_tardis.py b/benchmarks/run_tardis.py new file mode 100644 index 00000000000..8fca11f6030 --- /dev/null +++ b/benchmarks/run_tardis.py @@ -0,0 +1,25 @@ +""" +Basic TARDIS Benchmark. +""" + +from benchmarks.benchmark_base import BenchmarkBase +from tardis import run_tardis +from tardis.io.configuration.config_reader import Configuration + + +class BenchmarkRunTardis(BenchmarkBase): + """ + Class to benchmark the `run tardis` function. + """ + + def __init__(self): + super().__init__() + self.config = None + + def setup(self): + filename = "data/tardis_configv1_benchmark.yml" + path = self.get_relative_path(filename) + self.config = Configuration.from_yaml(path) + + def time_run_tardis(self): + run_tardis(self.config, log_level="ERROR", show_progress_bars=False) diff --git a/benchmarks/util/__init__.py b/benchmarks/util/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/benchmarks/util/base.py b/benchmarks/util/base.py new file mode 100644 index 00000000000..7dd61629e3e --- /dev/null +++ b/benchmarks/util/base.py @@ -0,0 +1,20 @@ +from os.path import dirname, realpath, join +from pathlib import Path, PosixPath + + +class Base: + @staticmethod + def get_path(partial_path: str) -> Path: + base_path = dirname(realpath(__file__)) + path = Path(base_path) / Path(partial_path) + return path + + @property + def tardis_ref_path(self) -> Path: + # TODO: This route is fixed but needs to get from the arguments given in the command line. + # /app/tardis-refdata + return Path("/app/tardis-refdata") + + @property + def example_configuration_dir(self) -> Path: + return self.get_path("../../tardis/io/configuration/tests/data") diff --git a/benchmarks/util/nlte.py b/benchmarks/util/nlte.py new file mode 100644 index 00000000000..7cfed4ef4dd --- /dev/null +++ b/benchmarks/util/nlte.py @@ -0,0 +1,78 @@ +from collections import OrderedDict +from copy import deepcopy +from pathlib import Path + +from benchmarks.util.base import Base +from tardis.io.atom_data import AtomData +from tardis.io.configuration.config_reader import Configuration +from tardis.io.util import yaml_load_file, YAMLLoader +from tardis.model import SimulationState + + +class NLTE: + def __init__(self): + self.__base = Base() + + @property + def tardis_config_verysimple_nlte(self) -> OrderedDict: + path: str = ( + "../../tardis/io/configuration/tests/data/tardis_configv1_nlte.yml" + ) + filename: Path = self.__base.get_path(path) + + return yaml_load_file( + filename, + YAMLLoader, + ) + + @property + def nlte_raw_model_root(self) -> SimulationState: + return SimulationState.from_config( + self.tardis_model_config_nlte_root, self.nlte_atom_data + ) + + @property + def nlte_raw_model_lu(self) -> SimulationState: + return SimulationState.from_config( + self.tardis_model_config_nlte_lu, self.nlte_atom_data + ) + + @property + def nlte_atom_data(self) -> AtomData: + atomic_data = deepcopy(self.nlte_atomic_dataset) + return atomic_data + + @property + def nlte_atomic_dataset(self) -> AtomData: + nlte_atomic_data = AtomData.from_hdf(self.nlte_atomic_data_fname) + return nlte_atomic_data + + @property + def nlte_atomic_data_fname(self) -> str: + atomic_data_fname = ( + f"{self.__base.tardis_ref_path}/nlte_atom_data/TestNLTE_He_Ti.h5" + ) + + if not Path(atomic_data_fname).exists(): + atom_data_missing_str = ( + f"Atomic datafiles {atomic_data_fname} does not seem to exist" + ) + raise Exception(atom_data_missing_str) + + return atomic_data_fname + + @property + def tardis_model_config_nlte_root(self) -> Configuration: + config = Configuration.from_yaml( + f"{self.__base.example_configuration_dir}/tardis_configv1_nlte.yml" + ) + config.plasma.nlte_solver = "root" + return config + + @property + def tardis_model_config_nlte_lu(self) -> Configuration: + config = Configuration.from_yaml( + f"{self.__base.example_configuration_dir}/tardis_configv1_nlte.yml" + ) + config.plasma.nlte_solver = "lu" + return config diff --git a/docs/contributing/development/benchmarks.rst b/docs/contributing/development/benchmarks.rst new file mode 100644 index 00000000000..f94d85301f3 --- /dev/null +++ b/docs/contributing/development/benchmarks.rst @@ -0,0 +1,89 @@ +.. _benchmarks: + +********** +Benchmarks +********** + +The objective of the benchmarking system is to detect regressions +that affect the performance of the TARDIS. This means we can visually +check whether there is a positive or negative spike in the TARDIS' performance. + + +AirSpeed Velocity (``ASV``) +=========================== + +TARDIS bases its benchmarking system on +`AirSpeed Velocity `_ (``ASV``). +Since it has a great advantage, which is that it is designed to run benchmarks on +`random servers `_ +such as those provided by +`GitHub Actions `_. +ASV eliminates the noise due to the technical differences between the servers +and produces graphs that indicate whether there is a regression or not. +It indicates the commit for the changes added or removed that affected performance +in some functions. + + +Installation +============ + +The complete guide is on the +`official ASV site `_, +however, here is detailed and summarized information to configure TARDIS with ASV. + +ASV needs `Conda `_ +(or `Miniconda `_) +and `Mamba `_. +To make configuration easier, you can use +`Mini-forge `_, +which includes the installers mentioned above. +In this step, Mamba installs Python, ASV, and Mamba; +however, this step does not configure Mamba with the TARDIS. + +.. code-block:: shell + + > export MAMBA_ENV_NAME="benchmark" + > mamba create --yes --name "${MAMBA_ENV_NAME}" python asv mamba + > mamba init + + +Set up +====== + +In this step, ASV configures TARDIS through Mamba. +Packages that use TARIDS are downloaded here. +These packages are mainly found in this ``tardis_env3.yml`` file. +The environment is also configured for ASV to execute benchmarks +and store the results through the ``asv.conf.json`` file. + +.. code-block:: shell + + > cd tardis + > export MAMBA_ENV_NAME="benchmark" + > mamba activate "${MAMBA_ENV_NAME}" + > asv setup + > asv machine --yes + + +Execution +========= + +ASV commands are used for execution. The first ``run`` command execute +the benchmarks found in the Python files that are in the ``benchmarks/`` +folder. Subsequently, the data and information are stored in the ``.asv/`` folder. + +.. code-block:: shell + + > cd tardis + > export MAMBA_ENV_NAME="benchmark" + > mamba activate "${MAMBA_ENV_NAME}" + > asv run + > asv publish + + +Visualization +============= + +There are two ways to view the data. The simplest thing is +to execute the ``asv preview`` command, creating a local web server. +The second is to run a local web server of your choice. diff --git a/docs/contributing/development/index.rst b/docs/contributing/development/index.rst index fa1d4666c8f..53fde2cc824 100644 --- a/docs/contributing/development/index.rst +++ b/docs/contributing/development/index.rst @@ -17,6 +17,7 @@ to the Astropy team for designing it. git_workflow documentation_guidelines running_tests + benchmarks code_quality developer_faq diff --git a/tardis/io/configuration/tests/data/tardis_configv1_verysimple_logger.yml b/tardis/io/configuration/tests/data/tardis_configv1_verysimple_logger.yml index 054deb5dd62..1c4110dac0a 100644 --- a/tardis/io/configuration/tests/data/tardis_configv1_verysimple_logger.yml +++ b/tardis/io/configuration/tests/data/tardis_configv1_verysimple_logger.yml @@ -5,6 +5,7 @@ supernova: time_explosion: 13 day atom_data: kurucz_atom_pure_simple.h5 + model: structure: type: specific diff --git a/tardis/io/model/readers/artis.py b/tardis/io/model/readers/artis.py index d013867cfde..21a353d0d94 100644 --- a/tardis/io/model/readers/artis.py +++ b/tardis/io/model/readers/artis.py @@ -50,7 +50,11 @@ def read_artis_density(fname): usecols=(0, 1, 2, 4, 5, 6, 7), dtype={item: np.float64 for item in artis_model_columns}, names=artis_model_columns, - delim_whitespace=True, + # The argument `delim_whitespace` was changed to `sep` + # because the first one is deprecated since version 2.2.0. + # The regular expression means: the separation is one or + # more spaces together (simple space, tabs, new lines). + sep=r"\s+", ).to_records(index=False) velocity = u.Quantity(artis_model["velocities"], "km/s").to("cm/s") diff --git a/tardis/io/model/readers/blondin_toymodel.py b/tardis/io/model/readers/blondin_toymodel.py index be2437e9b9e..d58ba5d2250 100644 --- a/tardis/io/model/readers/blondin_toymodel.py +++ b/tardis/io/model/readers/blondin_toymodel.py @@ -1,10 +1,8 @@ import re -import yaml - import numpy as np import pandas as pd - +import yaml from astropy import units as u from tardis.util.base import parse_quantity @@ -44,7 +42,15 @@ def read_blondin_toymodel(fname): ] raw_blondin_csv = pd.read_csv( - fname, delim_whitespace=True, comment="#", header=None, names=columns + fname, + # The argument `delim_whitespace` was changed to `sep` + # because the first one is deprecated since version 2.2.0. + # The regular expression means: the separation is one or + # more spaces together (simple space, tabs, new lines). + sep=r"\s+", + comment="#", + header=None, + names=columns, ) raw_blondin_csv.set_index("idx", inplace=True) diff --git a/tardis/io/model/readers/stella.py b/tardis/io/model/readers/stella.py index d198b614677..e314e716e3b 100644 --- a/tardis/io/model/readers/stella.py +++ b/tardis/io/model/readers/stella.py @@ -74,7 +74,11 @@ def read_stella_model(fname): # and the actual data data = pd.read_csv( fname, - delim_whitespace=True, + # The argument `delim_whitespace` was changed to `sep` + # because the first one is deprecated since version 2.2.0. + # The regular expression means: the separation is one or + # more spaces together (simple space, tabs, new lines). + sep=r"\s+", skiprows=DATA_START_ROW + 1, header=None, index_col=0, diff --git a/tardis/io/model/readers/tests/data/tardis_configv1_ascii_density_abund.yml b/tardis/io/model/readers/tests/data/tardis_configv1_ascii_density_abund.yml index 235546bbf41..99d8f1eb7b9 100644 --- a/tardis/io/model/readers/tests/data/tardis_configv1_ascii_density_abund.yml +++ b/tardis/io/model/readers/tests/data/tardis_configv1_ascii_density_abund.yml @@ -11,7 +11,6 @@ supernova: atom_data: kurucz_atom_pure_simple.h5 model: - structure: type: file filename: density.dat diff --git a/tardis/io/model/readers/tests/data/tardis_configv1_isotope_iabund.yml b/tardis/io/model/readers/tests/data/tardis_configv1_isotope_iabund.yml index f70ee7fa60e..7040021ecd9 100644 --- a/tardis/io/model/readers/tests/data/tardis_configv1_isotope_iabund.yml +++ b/tardis/io/model/readers/tests/data/tardis_configv1_isotope_iabund.yml @@ -1,45 +1,45 @@ tardis_config_version: v1.0 supernova: - luminosity_requested: 2.8e9 solLum - time_explosion: 13 day + luminosity_requested: 2.8e9 solLum + time_explosion: 13 day atom_data: kurucz_atom_pure_simple.h5 model: - structure: - type: specific - velocity: - start: 1.1e4 km/s - stop: 2.0e4 km/s - num: 2 - density: - type: branch85_w7 - abundances: - type: file - filename: non_uniform_isotope_abundance.dat - filetype: custom_composition + structure: + type: specific + velocity: + start: 1.1e4 km/s + stop: 2.0e4 km/s + num: 2 + density: + type: branch85_w7 + abundances: + type: file + filename: non_uniform_isotope_abundance.dat + filetype: custom_composition plasma: - ionization: lte - excitation: lte - radiative_rates_type: dilute-blackbody - line_interaction_type: macroatom + ionization: lte + excitation: lte + radiative_rates_type: dilute-blackbody + line_interaction_type: macroatom montecarlo: - seed: 23111963 - no_of_packets: 2.0e+5 - iterations: 5 - last_no_of_packets: 5.0e+5 - no_of_virtual_packets: 5 - convergence_strategy: - type: damped - damping_constant: 0.5 - threshold: 0.05 - lock_t_inner_cycles: 1 - t_inner_update_exponent: -0.5 + seed: 23111963 + no_of_packets: 2.0e+5 + iterations: 5 + last_no_of_packets: 5.0e+5 + no_of_virtual_packets: 5 + convergence_strategy: + type: damped + damping_constant: 0.5 + threshold: 0.05 + lock_t_inner_cycles: 1 + t_inner_update_exponent: -0.5 spectrum: - start: 500 angstrom - stop: 20000 angstrom - num: 10000 + start: 500 angstrom + stop: 20000 angstrom + num: 10000 diff --git a/tardis/io/model/readers/tests/data/tardis_configv1_isotope_uniabund.yml b/tardis/io/model/readers/tests/data/tardis_configv1_isotope_uniabund.yml index 02d2f7f9a60..7f1c583a2f9 100755 --- a/tardis/io/model/readers/tests/data/tardis_configv1_isotope_uniabund.yml +++ b/tardis/io/model/readers/tests/data/tardis_configv1_isotope_uniabund.yml @@ -5,8 +5,8 @@ supernova: time_explosion: 13 day atom_data: kurucz_atom_pure_simple.h5 -model: +model: structure: type: specific velocity: diff --git a/tardis/plasma/tests/data/plasma_base_test_config.yml b/tardis/plasma/tests/data/plasma_base_test_config.yml index f17612618da..b8e06f580e8 100644 --- a/tardis/plasma/tests/data/plasma_base_test_config.yml +++ b/tardis/plasma/tests/data/plasma_base_test_config.yml @@ -7,18 +7,14 @@ supernova: atom_data: kurucz_atom_pure_simple.h5 model: - structure: type: specific - velocity: start: 1.1e4 km/s stop: 2.0e4 km/s num: 20 - density: type: branch85_w7 - abundances: type: uniform He: 1 diff --git a/tardis/util/base.py b/tardis/util/base.py index 68693c3fa6b..8e21da9d667 100644 --- a/tardis/util/base.py +++ b/tardis/util/base.py @@ -32,7 +32,11 @@ ATOMIC_SYMBOLS_DATA = ( pd.read_csv( get_internal_data_path("atomic_symbols.dat"), - delim_whitespace=True, + # The argument `delim_whitespace` was changed to `sep` + # because the first one is deprecated since version 2.2.0. + # The regular expression means: the separation is one or + # more spaces together (simple space, tabs, new lines). + sep=r"\s+", names=["atomic_number", "symbol"], ) .set_index("atomic_number")