diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 0b23826..28da8bf 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -20,10 +20,8 @@ jobs: python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] os: [ubuntu-22.04] include: - - python-version: "3.11" - os: macos-12 - python-version: "3.12" - os: macos-12 + os: macos-13 steps: - uses: actions/checkout@v4 - name: Set up Python ${{ matrix.python-version }} @@ -41,7 +39,14 @@ jobs: - name: Run Tests env: SIMWEIGHTS_TESTDATA: . - run: python3 -m pytest --cov-report=xml + run: python3 -m pytest --cov-report=xml --junit-xml=test-results-${{matrix.os}}-${{matrix.python-version}}.junit.xml + - name: Upload Test Results + uses: actions/upload-artifact@v4 + if: always() + with: + if-no-files-found: error + name: test-results-${{matrix.os}}-${{matrix.python-version}} + path: test-results-${{matrix.os}}-${{matrix.python-version}}.junit.xml - name: Upload Coverage to Codecov if: ${{ !github.event.act }} uses: codecov/codecov-action@v4 @@ -50,3 +55,24 @@ jobs: with: fail_ci_if_error: false verbose: true + publish-test-results: + name: "Publish Tests Results" + needs: Tests + runs-on: ubuntu-latest + permissions: + checks: write + pull-requests: write + contents: read + if: always() + steps: + - name: Download Artifacts + uses: actions/download-artifact@v4 + with: + path: . + pattern: test-results-* + merge-multiple: true + - name: Publish Test Results + uses: EnricoMi/publish-unit-test-result-action@v2 + with: + files: "*.xml" + deduplicate_classes_by_file_name: true diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index de7e6ee..6e70819 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -44,7 +44,7 @@ repos: exclude: ^contrib/ additional_dependencies: [numpy, pandas] - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.1.14 + rev: v0.2.1 hooks: - id: ruff args: [--fix, --show-fixes] diff --git a/.reuse/dep5 b/.reuse/dep5 index c51c9ff..75aa168 100644 --- a/.reuse/dep5 +++ b/.reuse/dep5 @@ -3,10 +3,10 @@ Upstream-Name: SimWeights Upstream-Contact: the SimWeights contributors Source: https://github.com/icecube/simweights -Files: docs/*.svg +Files: docs/*.svg tests/flux_values.json Copyright: 2022 the SimWeights contributors License: BSD-2-Clause -Files: tests/flux_values.json -Copyright: 2022 the SimWeights contributors +Files: tests/test_fluxes/*.npz +Copyright: 2024 the SimWeights contributors License: BSD-2-Clause diff --git a/pyproject.toml b/pyproject.toml index e9d46af..ab70d5c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,7 +34,15 @@ requires-python = "~=3.7" [project.optional-dependencies] dev = ["pytest", "pre-commit", "reuse", "black", "ruff", "pylint", "mypy"] docs = ["sphinx", "sphinx-rtd-theme", "pandas"] -test = ["h5py", "tables < 3.8; python_version < '3.9'", "tables <= 3.9.2", "pandas", "uproot", "pytest-cov"] +test = [ + "h5py", + "tables < 3.8; python_version < '3.9'", + "tables", + "pandas", + "uproot", + "pytest-cov", + 'pytest-regressions' +] [project.scripts] simweights = "simweights.cmdline:main" @@ -89,7 +97,6 @@ target-version = "py38" fixable = ["I"] ignore = [ "ANN401", # any-type - "FBT", # flake8-boolean-trap "S101", # assert-used "COM812", # confilcts with ruff formatter "ISC001", # confilcts with ruff formatter @@ -118,3 +125,15 @@ select = ["ALL"] [tool.ruff.lint.pydocstyle] convention = "google" + +[tool.tox] +legacy_tox_ini = """ +[tox] +envlist = py3{8,9,10,11,12} +isolated_build = True + +[testenv] +passenv = SIMWEIGHTS_TESTDATA, HDF5_DIR +deps = .[test] +commands = pytest +""" diff --git a/src/simweights/__init__.py b/src/simweights/__init__.py index 66c631d..ed0cdcf 100644 --- a/src/simweights/__init__.py +++ b/src/simweights/__init__.py @@ -28,6 +28,7 @@ "GaisserHillas", "GlobalFitGST", "GlobalFitGST_IT", + "GlobalSplineFit", "GlobalSplineFit5Comp", "GlobalSplineFit_IT", "Hoerandel", @@ -54,6 +55,7 @@ GaisserHillas, GlobalFitGST, GlobalFitGST_IT, + GlobalSplineFit, GlobalSplineFit5Comp, GlobalSplineFit_IT, Hoerandel, diff --git a/src/simweights/_corsika_weighter.py b/src/simweights/_corsika_weighter.py index e90c969..430619b 100644 --- a/src/simweights/_corsika_weighter.py +++ b/src/simweights/_corsika_weighter.py @@ -13,11 +13,11 @@ from ._generation_surface import GenerationSurface, generation_surface from ._powerlaw import PowerLaw from ._spatial import NaturalRateCylinder -from ._utils import Column, constcol, get_column, get_table, has_table +from ._utils import Column, constcol, get_column, get_table, has_column, has_table from ._weighter import Weighter -def sframe_corsika_surface(table: Any, oversampling: bool) -> GenerationSurface: +def sframe_corsika_surface(table: Any) -> GenerationSurface: """Inspect the rows of a CORSIKA S-Frame table object to generate a surface object. This function works on files generated with either triggered CORSIKA or corsika-reader because @@ -40,7 +40,7 @@ def sframe_corsika_surface(table: Any, oversampling: bool) -> GenerationSurface: get_column(table, "max_energy")[i], "energy", ) - oversampling_val = get_column(table, "oversampling")[i] if oversampling else 1 + oversampling_val = get_column(table, "oversampling")[i] if has_column(table, "oversampling") else 1 pdgid = int(get_column(table, "primary_type")[i]) surfaces.append( get_column(table, "n_events")[i] @@ -108,7 +108,7 @@ def CorsikaWeighter(file_obj: Any, nfiles: float | None = None) -> Weighter: # ) raise RuntimeError(mesg) - surface = sframe_corsika_surface(get_table(file_obj, info_obj), oversampling=False) + surface = sframe_corsika_surface(get_table(file_obj, info_obj)) triggered = True elif nfiles is None: @@ -119,7 +119,7 @@ def CorsikaWeighter(file_obj: Any, nfiles: float | None = None) -> Weighter: # "nfiles and no I3CorsikaInfo table was found." ) raise RuntimeError(msg) - surface = sframe_corsika_surface(get_table(file_obj, info_obj), oversampling=True) + surface = sframe_corsika_surface(get_table(file_obj, info_obj)) triggered = False else: diff --git a/src/simweights/_fluxes.py b/src/simweights/_fluxes.py index 5207960..dc44b15 100644 --- a/src/simweights/_fluxes.py +++ b/src/simweights/_fluxes.py @@ -21,7 +21,7 @@ from pathlib import Path from typing import TYPE_CHECKING, Callable, Mapping, Sequence -from numpy import asfarray, bool_, broadcast_arrays, exp, float64, genfromtxt, int32, piecewise, sqrt +from numpy import asarray, bool_, broadcast_arrays, exp, float64, genfromtxt, int32, piecewise, sqrt from numpy import sum as nsum from scipy.interpolate import CubicSpline # pylint: disable=import-error @@ -398,7 +398,7 @@ def __init__( self: FixedFractionFlux, fractions: Mapping[PDGCode, float], basis: CosmicRayFlux | None = None, - normalized: bool = True, + normalized: bool = True, # noqa: FBT001, FBT002 ) -> None: """Flux that is a fixed fraction of another flux. @@ -422,7 +422,7 @@ def __call__(self: FixedFractionFlux, energy: ArrayLike, pdgid: ArrayLike) -> ND energy_arr, pdgid_arr = broadcast_arrays(energy, pdgid) fluxsum = sum(self.flux(energy_arr, p) for p in self.pdgids) cond = self._condition(energy_arr, pdgid_arr) - return asfarray(fluxsum * piecewise(energy, cond, self.fracs)) + return asarray(fluxsum * piecewise(energy, cond, self.fracs), dtype=float64) class _references: diff --git a/src/simweights/_generation_surface.py b/src/simweights/_generation_surface.py index 58dd694..24bbffb 100644 --- a/src/simweights/_generation_surface.py +++ b/src/simweights/_generation_surface.py @@ -101,13 +101,13 @@ def get_epdf( cols = {} shape = None for key, value in kwargs.items(): - cols[key] = np.asfarray(value) + cols[key] = np.asarray(value, dtype=np.float64) if shape is None: shape = cols[key].shape else: assert shape == cols[key].shape # type: ignore[unreachable] assert shape is not None - count = np.zeros(shape, dtype=np.float_) + count = np.zeros(shape, dtype=np.float64) # loop over particle type for ptype in np.unique(pdgid): mask = ptype == pdgid diff --git a/src/simweights/_powerlaw.py b/src/simweights/_powerlaw.py index 16e78f4..7786ef7 100644 --- a/src/simweights/_powerlaw.py +++ b/src/simweights/_powerlaw.py @@ -52,17 +52,17 @@ def __init__(self: PowerLaw, g: float, a: float, b: float, colname: str | None = self.columns = (colname,) def _pdf(self: PowerLaw, x: NDArray[np.float64]) -> NDArray[np.float64]: - return np.asfarray(x**self.g / self.integral) + return np.asarray(x**self.g / self.integral, dtype=np.float64) def _cdf(self: PowerLaw, x: NDArray[np.float64]) -> NDArray[np.float64]: if self.G == 0: - return np.asfarray(np.log(x / self.a) / self.integral) - return np.asfarray((x**self.G - self.a**self.G) / self.G / self.integral) + return np.asarray(np.log(x / self.a) / self.integral, dtype=np.float64) + return np.asarray((x**self.G - self.a**self.G) / self.G / self.integral, dtype=np.float64) def _ppf(self: PowerLaw, q: NDArray[np.float64]) -> NDArray[np.float64]: if self.G == 0: - return np.asfarray(self.a * np.exp(q * self.integral)) - return np.asfarray((q * self.G * self.integral + self.a**self.G) ** (1 / self.G)) + return np.asarray(self.a * np.exp(q * self.integral), dtype=np.float64) + return np.asarray((q * self.G * self.integral + self.a**self.G) ** (1 / self.G), np.float64) def pdf(self: PowerLaw, x: ArrayLike) -> NDArray[np.float64]: r"""Probability density function. @@ -74,7 +74,7 @@ def pdf(self: PowerLaw, x: ArrayLike) -> NDArray[np.float64]: Returns: array_like: Probability density function evaluated at `x` """ - xa = np.asfarray(x) + xa = np.asarray(x, dtype=np.float64) return np.piecewise(xa, [(xa >= self.a) & (xa <= self.b)], [self._pdf]) def cdf(self: PowerLaw, x: ArrayLike) -> NDArray[np.float64]: @@ -86,7 +86,7 @@ def cdf(self: PowerLaw, x: ArrayLike) -> NDArray[np.float64]: Returns: array_like: Cumulative distribution function evaluated at `x` """ - qa = np.asfarray(x) + qa = np.asarray(x, dtype=np.float64) return np.piecewise(qa, [qa < self.a, qa > self.b], [0, 1, self._cdf]) def ppf(self: PowerLaw, q: ArrayLike) -> NDArray[np.float64]: @@ -98,7 +98,7 @@ def ppf(self: PowerLaw, q: ArrayLike) -> NDArray[np.float64]: Returns: array_like: quantile corresponding to the lower tail probability `q`. """ - qa = np.asfarray(q) + qa = np.asarray(q, dtype=np.float64) return np.piecewise(qa, [(qa >= 0) & (qa <= 1)], [self._ppf, np.nan]) def rvs(self: PowerLaw, size: Any = None, random_state: SeedType = None) -> NDArray[np.float64]: @@ -116,7 +116,7 @@ def rvs(self: PowerLaw, size: Any = None, random_state: SeedType = None) -> NDAr Default is None. """ rand_state = check_random_state(random_state) - return self._ppf(np.asfarray(rand_state.uniform(0, 1, size))) + return self._ppf(np.asarray(rand_state.uniform(0, 1, size), dtype=np.float64)) def __repr__(self: PowerLaw) -> str: return f"{self.__class__.__name__}({self.g} ,{self.a}, {self.b})" diff --git a/src/simweights/_spatial.py b/src/simweights/_spatial.py index 3ecc6a1..e58df28 100644 --- a/src/simweights/_spatial.py +++ b/src/simweights/_spatial.py @@ -42,17 +42,18 @@ def projected_area(self: CylinderBase, cos_zen: ArrayLike) -> NDArray[np.float64 As seen from the angle described by cos_zen. """ - cosz = np.asfarray(cos_zen) + cosz = np.asarray(cos_zen, dtype=np.float64) assert np.all(cosz >= -1) assert np.all(cosz <= +1) - return np.asfarray(self._cap * np.abs(cosz) + self._side * np.sqrt(1 - cosz**2)) + return np.asarray(self._cap * np.abs(cosz) + self._side * np.sqrt(1 - cosz**2), dtype=np.float64) def _diff_etendue(self: CylinderBase, cos_zen: ArrayLike) -> NDArray[np.float64]: - cosz = np.asfarray(cos_zen) + cosz = np.asarray(cos_zen, dtype=np.float64) assert np.all(cosz >= -1) assert np.all(cosz <= +1) - return np.asfarray( + return np.asarray( np.pi * (self._cap * cosz * np.abs(cosz) + self._side * (cosz * np.sqrt(1 - cosz**2) - np.arccos(cosz))), + dtype=np.float64, ) def pdf(self: CylinderBase, cos_zen: ArrayLike) -> NDArray[np.float64]: @@ -91,7 +92,7 @@ def _pdf(self: UniformSolidAngleCylinder, cos_zen: NDArray[np.float64]) -> NDArr return 1 / (2 * np.pi * (self.cos_zen_max - self.cos_zen_min) * self.projected_area(cos_zen)) def pdf(self: UniformSolidAngleCylinder, cos_zen: ArrayLike) -> NDArray[np.float64]: - cosz = np.asfarray(cos_zen) + cosz = np.asarray(cos_zen, dtype=np.float64) return np.piecewise(cosz, [(cosz >= self.cos_zen_min) & (cosz <= self.cos_zen_max)], [self._pdf]) @@ -123,7 +124,7 @@ def __init__( self._normalization = 1 / self.etendue def pdf(self: NaturalRateCylinder, cos_zen: ArrayLike) -> NDArray[np.float64]: - cosz = np.asfarray(cos_zen) + cosz = np.asarray(cos_zen, dtype=np.float64) return np.piecewise( cosz, [(cosz >= self.cos_zen_min) & (cosz <= self.cos_zen_max)], @@ -157,7 +158,7 @@ def projected_area(self: CircleInjector, cos_zen: float) -> float: # noqa: ARG0 def pdf(self: CircleInjector, cos_zen: ArrayLike) -> NDArray[np.float64]: """The probability density function for the given zenith angle.""" - cosz = np.asfarray(cos_zen) + cosz = np.asarray(cos_zen, dtype=np.float64) return np.piecewise( cosz, [(cosz >= self.cos_zen_min) & (cosz <= self.cos_zen_max)], diff --git a/src/simweights/_utils.py b/src/simweights/_utils.py index 6759951..c63cbe4 100644 --- a/src/simweights/_utils.py +++ b/src/simweights/_utils.py @@ -29,7 +29,7 @@ def __init__(self: Column, colname: str | None = None) -> None: def pdf(self: Column, value: ArrayLike) -> NDArray[np.float64]: r"""Probability density function.""" - return 1 / np.asfarray(value) + return 1 / np.asarray(value, dtype=np.float64) def __eq__(self: Column, other: object) -> bool: return isinstance(other, Column) and self.columns == other.columns @@ -47,7 +47,7 @@ def __init__(self: Const, v: float) -> None: def pdf(self: Const) -> NDArray[np.float64]: r"""Probability density function.""" - return np.asfarray(self.v) + return np.asarray(self.v, dtype=np.float64) def __eq__(self: Const, other: object) -> bool: return isinstance(other, Const) and self.v == other.v @@ -84,11 +84,11 @@ def has_column(table: Any, name: str) -> bool: def get_column(table: Any, name: str) -> NDArray[np.float64]: """Helper function getting a column from a table, works with h5py, pytables, and pandas.""" if hasattr(table, "cols"): - return np.asfarray(getattr(table.cols, name)[:]) + return np.asarray(getattr(table.cols, name)[:], dtype=np.float64) column = table[name] if hasattr(column, "array") and callable(column.array): - return np.asfarray(column.array(library="np")) - return np.asfarray(column) + return np.asarray(column.array(library="np"), dtype=np.float64) + return np.asarray(column, dtype=np.float64) def constcol(table: Any, colname: str, mask: NDArray[np.bool_] | None = None) -> float: diff --git a/src/simweights/_weighter.py b/src/simweights/_weighter.py index 36c76db..2896109 100644 --- a/src/simweights/_weighter.py +++ b/src/simweights/_weighter.py @@ -199,7 +199,7 @@ def effective_area( assert np.array_equal(enbin, energy_bins) assert np.array_equal(czbin, cos_zenith_bins) e_width, z_width = np.meshgrid(np.ediff1d(enbin), np.ediff1d(czbin)) - return np.asfarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies)) + return np.asarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies), dtype=np.float64) def __add__(self: Weighter, other: Weighter | int) -> Weighter: if other == 0: diff --git a/tests/test_corsika_datasets.py b/tests/test_corsika_datasets.py index f949de2..77996ac 100755 --- a/tests/test_corsika_datasets.py +++ b/tests/test_corsika_datasets.py @@ -5,135 +5,119 @@ # SPDX-License-Identifier: BSD-2-Clause import os -import unittest +import sys +from pathlib import Path import h5py import numpy as np import pandas as pd +import pytest import tables import uproot from simweights import CorsikaWeighter, GaisserH4a from simweights._utils import constcol - -class TestCorsikaDatasets(unittest.TestCase): - @classmethod - def setUpClass(cls): - cls.flux = GaisserH4a() - datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None) - if not datadir: - cls.skipTest(None, "environment variable SIMWEIGHTS_TESTDATA not set") - cls.datadir = datadir + "/" - - def untriggered_weights(self, f): - cwm = f["CorsikaWeightMap"] - pflux = self.flux(cwm["PrimaryEnergy"], cwm["PrimaryType"]) - if (cwm["PrimarySpectralIndex"] == -1).any(): - assert (cwm["PrimarySpectralIndex"] == -1).all() - energy_integral = np.log(cwm["EnergyPrimaryMax"] / cwm["EnergyPrimaryMin"]) - else: - energy_integral = ( - cwm["EnergyPrimaryMax"] ** (cwm["PrimarySpectralIndex"] + 1) - - cwm["EnergyPrimaryMin"] ** (cwm["PrimarySpectralIndex"] + 1) - ) / (cwm["PrimarySpectralIndex"] + 1) - - energy_weight = cwm["PrimaryEnergy"] ** cwm["PrimarySpectralIndex"] - return 1e4 * pflux * energy_integral / energy_weight * cwm["AreaSum"] / (cwm["NEvents"] * cwm["OverSampling"]) - - def triggered_weights(self, f): - i3cw = f["I3CorsikaWeight"] - flux_val = self.flux(i3cw["energy"], i3cw["type"]) - info = f["I3PrimaryInjectorInfo"] - energy = i3cw["energy"] - epdf = np.zeros_like(energy, dtype=float) - - for ptype in np.unique(info["primary_type"]): - info_mask = info["primary_type"] == ptype - n_events = info["n_events"][info_mask].sum() - min_energy = constcol(info, "min_energy", info_mask) - max_energy = constcol(info, "max_energy", info_mask) - min_zenith = constcol(info, "min_zenith", info_mask) - max_zenith = constcol(info, "max_zenith", info_mask) - cylinder_height = constcol(info, "cylinder_height", info_mask) - cylinder_radius = constcol(info, "cylinder_radius", info_mask) - power_law_index = constcol(info, "power_law_index", info_mask) - - G = power_law_index + 1 - side = 2e4 * cylinder_radius * cylinder_height - cap = 1e4 * np.pi * cylinder_radius**2 - cos_minz = np.cos(min_zenith) - cos_maxz = np.cos(max_zenith) - ET1 = cap * cos_minz * np.abs(cos_minz) + side * (cos_minz * np.sqrt(1 - cos_minz**2) - min_zenith) - ET2 = cap * cos_maxz * np.abs(cos_maxz) + side * (cos_maxz * np.sqrt(1 - cos_maxz**2) - max_zenith) - etendue = np.pi * (ET1 - ET2) - - mask = ptype == i3cw["type"] - energy_term = energy[mask] ** power_law_index * G / (max_energy**G - min_energy**G) - epdf[mask] += n_events * energy_term / etendue - - return i3cw["weight"] * flux_val / epdf - - def cmp_dataset(self, triggered, fname, rate): - fname = self.datadir + fname - - if triggered: - nfiles = None - refweight = self.triggered_weights - else: - nfiles = 1 - refweight = self.untriggered_weights - - reffile = h5py.File(fname + ".hdf5", "r") - w0 = refweight(reffile) - - inputfiles = [ - ("h5py", reffile), - ("uproot", uproot.open(fname + ".root")), - ("tables", tables.open_file(fname + ".hdf5", "r")), - ("pandas", pd.HDFStore(fname + ".hdf5", "r")), - ] - - for name, infile in inputfiles: - with self.subTest(name): - wobj = CorsikaWeighter(infile, nfiles) - w = wobj.get_weights(self.flux) - self.assertAlmostEqual(w.sum(), rate) - np.testing.assert_allclose(w0, w, 1e-6) - - for _, infile in inputfiles: - infile.close() - - def test_012602(self): - self.cmp_dataset(False, "Level2_IC86.2015_corsika.012602.000000", 102.01712611701736) - - def test_020014(self): - self.cmp_dataset(False, "Level2_IC86.2015_corsika.020014.000000", 23.015500214424705) - - def test_020021(self): - self.cmp_dataset(False, "Level2_IC86.2015_corsika.020021.000000", 69.75465614509928) - - def test_020208(self): - self.cmp_dataset(False, "Level2_IC86.2016_corsika.020208.000001", 22.622983704306385) - - def test_020243(self): - self.cmp_dataset(False, "Level2_IC86.2016_corsika.020243.000001", 4.590586137762489) - - def test_020263(self): - self.cmp_dataset(False, "Level2_IC86.2016_corsika.020263.000000", 10.183937153798436) - - def test_020777(self): - self.cmp_dataset(False, "Level2_IC86.2016_corsika.020777.000000", 362.94284441826704) - - def test_020778(self): - self.cmp_dataset(False, "Level2_IC86.2016_corsika.020778.000000", 6.2654796956603) - - def test_020780(self): - self.cmp_dataset(False, "Level2_IC86.2016_corsika.020780.000000", 14.215947086098588) - - def test_21889(self): - self.cmp_dataset(True, "Level2_IC86.2016_corsika.021889.000000", 122.83809329321922) +flux = GaisserH4a() +datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None) +if datadir: + datadir = Path(datadir) + +datasets = [ + (False, "Level2_IC86.2015_corsika.012602.000000", 102.01712611701736), + (False, "Level2_IC86.2015_corsika.020014.000000", 23.015500214424705), + (False, "Level2_IC86.2015_corsika.020021.000000", 69.75465614509928), + (False, "Level2_IC86.2016_corsika.020208.000001", 22.622983704306385), + (False, "Level2_IC86.2016_corsika.020243.000001", 4.590586137762489), + (False, "Level2_IC86.2016_corsika.020263.000000", 10.183937153798436), + (False, "Level2_IC86.2016_corsika.020777.000000", 362.94284441826704), + (False, "Level2_IC86.2016_corsika.020778.000000", 6.2654796956603), + (False, "Level2_IC86.2016_corsika.020780.000000", 14.215947086098588), + (True, "Level2_IC86.2016_corsika.021889.000000", 122.83809329321922), +] + + +def untriggered_weights(f): + cwm = f["CorsikaWeightMap"] + pflux = flux(cwm["PrimaryEnergy"], cwm["PrimaryType"]) + if (cwm["PrimarySpectralIndex"] == -1).any(): + assert (cwm["PrimarySpectralIndex"] == -1).all() + energy_integral = np.log(cwm["EnergyPrimaryMax"] / cwm["EnergyPrimaryMin"]) + else: + energy_integral = ( + cwm["EnergyPrimaryMax"] ** (cwm["PrimarySpectralIndex"] + 1) + - cwm["EnergyPrimaryMin"] ** (cwm["PrimarySpectralIndex"] + 1) + ) / (cwm["PrimarySpectralIndex"] + 1) + + energy_weight = cwm["PrimaryEnergy"] ** cwm["PrimarySpectralIndex"] + return 1e4 * pflux * energy_integral / energy_weight * cwm["AreaSum"] / (cwm["NEvents"] * cwm["OverSampling"]) + + +def triggered_weights(f): + i3cw = f["I3CorsikaWeight"] + flux_val = flux(i3cw["energy"], i3cw["type"]) + info = f["I3PrimaryInjectorInfo"] + energy = i3cw["energy"] + epdf = np.zeros_like(energy, dtype=float) + + for ptype in np.unique(info["primary_type"]): + info_mask = info["primary_type"] == ptype + n_events = info["n_events"][info_mask].sum() + min_energy = constcol(info, "min_energy", info_mask) + max_energy = constcol(info, "max_energy", info_mask) + min_zenith = constcol(info, "min_zenith", info_mask) + max_zenith = constcol(info, "max_zenith", info_mask) + cylinder_height = constcol(info, "cylinder_height", info_mask) + cylinder_radius = constcol(info, "cylinder_radius", info_mask) + power_law_index = constcol(info, "power_law_index", info_mask) + + G = power_law_index + 1 + side = 2e4 * cylinder_radius * cylinder_height + cap = 1e4 * np.pi * cylinder_radius**2 + cos_minz = np.cos(min_zenith) + cos_maxz = np.cos(max_zenith) + ET1 = cap * cos_minz * np.abs(cos_minz) + side * (cos_minz * np.sqrt(1 - cos_minz**2) - min_zenith) + ET2 = cap * cos_maxz * np.abs(cos_maxz) + side * (cos_maxz * np.sqrt(1 - cos_maxz**2) - max_zenith) + etendue = np.pi * (ET1 - ET2) + + mask = ptype == i3cw["type"] + energy_term = energy[mask] ** power_law_index * G / (max_energy**G - min_energy**G) + epdf[mask] += n_events * energy_term / etendue + + return i3cw["weight"] * flux_val / epdf + + +@pytest.mark.parametrize(("triggered", "fname", "rate"), datasets) +@pytest.mark.skipif(not datadir, reason="environment variable SIMWEIGHTS_TESTDATA not set") +def test_dataset(triggered, fname, rate): + fname = datadir / fname + + if triggered: + nfiles = None + refweight = triggered_weights + else: + nfiles = 1 + refweight = untriggered_weights + + reffile = h5py.File(str(fname) + ".hdf5", "r") + w0 = refweight(reffile) + + inputfiles = [ + ("h5py", reffile), + ("uproot", uproot.open(str(fname) + ".root")), + ("tables", tables.open_file(str(fname) + ".hdf5", "r")), + ("pandas", pd.HDFStore(str(fname) + ".hdf5", "r")), + ] + + for _, infile in inputfiles: + wobj = CorsikaWeighter(infile, nfiles) + w = wobj.get_weights(flux) + assert w.sum() == pytest.approx(rate) + assert w0 == pytest.approx(w, 1e-6) + + for _, infile in inputfiles: + infile.close() if __name__ == "__main__": - unittest.main() + sys.exit(pytest.main(["-v", __file__, *sys.argv[1:]])) diff --git a/tests/test_fluxes.py b/tests/test_fluxes.py index 2b18fdc..2c161e3 100755 --- a/tests/test_fluxes.py +++ b/tests/test_fluxes.py @@ -5,123 +5,74 @@ # SPDX-License-Identifier: BSD-2-Clause import json -import unittest +import sys from pathlib import Path import numpy as np +import pytest -from simweights import _fluxes +import simweights +with (Path(__file__).parent / "flux_values.json").open() as f: + flux_values = json.load(f) E = np.logspace(2, 10, 9) - -class TestCosmicRayModels(unittest.TestCase): - @classmethod - def setUpClass(cls): - with (Path(__file__).parent / "flux_values.json").open() as f: - cls.flux_values = json.load(f) - - def flux_cmp(self, name, *args): - flux = getattr(_fluxes, name)(*args) - v1 = flux(*np.meshgrid(E, [int(i) for i in self.flux_values[name]])) - v2 = np.array(list(self.flux_values[name].values())) / 1e4 - np.testing.assert_allclose(v1, v2, 1e-13) - - # make sure you get zero for non CR primaries - np.testing.assert_allclose(flux(E, 22), 0) - - def test_Hoerandel(self): - self.flux_cmp("Hoerandel") - - def test_Hoerandel5(self): - self.flux_cmp("Hoerandel5") - - def test_Hoerandel_IT(self): - self.flux_cmp("Hoerandel_IT") - - def test_GaisserHillas(self): - self.flux_cmp("GaisserHillas") - - def test_GaisserH3a(self): - self.flux_cmp("GaisserH3a") - - def test_GaisserH4a(self): - self.flux_cmp("GaisserH4a") - - def test_GaisserH4a_IT(self): - self.flux_cmp("GaisserH4a_IT") - - def test_Honda2004(self): - self.flux_cmp("Honda2004") - - def test_TIG1996(self): - self.flux_cmp("TIG1996") - - def test_GlobalFitGST(self): - self.flux_cmp("GlobalFitGST") - - def test_GlobalFitGST_IT(self): - self.flux_cmp("GlobalFitGST_IT") - - def test_GlobalSplineFit(self): - self.flux_cmp("GlobalSplineFit") - - def test_GlobalSplineFit5Comp(self): - self.flux_cmp("GlobalSplineFit5Comp") - - def test_GlobalSplineFit_IT(self): - self.flux_cmp("GlobalSplineFit_IT") - - def test_FixedFractionFlux(self): - self.flux_cmp("FixedFractionFlux", {2212: 0.1, 1000020040: 0.2, 1000080160: 0.3, 1000260560: 0.4}) - self.flux_cmp( - "FixedFractionFlux", - {2212: 0.1, 1000020040: 0.2, 1000080160: 0.3, 1000260560: 0.4}, - _fluxes.GaisserH4a_IT(), - ) - - -def test_GlobalSplineFit_similar(): - """ - Test if GlobalSplineFit is similar to to other models within a factor of 500, - mainly to check if the units provided by the GST files match. - This can be not transparent in the file and web-interface. - If the units mismatch, we should expect at least a deviation of 1000 (one prefix) - or most likely a mismatch of 10 000 (cm^2 vs m^2). - """ - gsf = _fluxes.GlobalSplineFit() - - for name in ("GlobalFitGST", "GaisserH3a", "Hoerandel5"): - model = getattr(_fluxes, name)() - - same_pdgids = [pdgid for pdgid in gsf.pdgids if pdgid in model.pdgids] - - f_gsf = gsf(*np.meshgrid(E, [int(i) for i in same_pdgids])) - f = model(*np.meshgrid(E, [int(i) for i in same_pdgids])) - - assert 0.2 < np.mean(f / f_gsf) < 500 - - -def test_GlobalSplineFit5Comp_similar(): +flux_models = [ + simweights.Hoerandel(), + simweights.Hoerandel5(), + simweights.Hoerandel_IT(), + simweights.GaisserHillas(), + simweights.GaisserH3a(), + simweights.GaisserH4a(), + simweights.GaisserH4a_IT(), + simweights.Honda2004(), + simweights.TIG1996(), + simweights.GlobalFitGST(), + simweights.GlobalFitGST_IT(), + simweights.GlobalSplineFit(), + simweights.GlobalSplineFit5Comp(), + simweights.GlobalSplineFit_IT(), + simweights.FixedFractionFlux({2212: 0.1, 1000020040: 0.2, 1000080160: 0.3, 1000260560: 0.4}), + simweights.FixedFractionFlux({2212: 0.1, 1000020040: 0.2, 1000080160: 0.3, 1000260560: 0.4}, simweights.GaisserH4a_IT()), +] + + +@pytest.mark.parametrize("flux", flux_models, ids=[x.__class__.__name__ for x in flux_models]) +def test_flux_model(flux, ndarrays_regression): + # this is the old regression test it can stick around for a bit but will be deleted at a certain point + for pdgid in flux.pdgids: + v1 = flux(E, pdgid) + v2 = np.array(flux_values[flux.__class__.__name__][str(int(pdgid))]) / 1e4 + assert v1 == pytest.approx(v2, rel=1e-13) + + ndarrays_regression.check({pdgid.name: flux(E, pdgid) for pdgid in flux.pdgids}, default_tolerance={"rtol": 1e-13}) + # make sure you get zero for non CR primaries + assert flux(E, 22) == pytest.approx(0) + + +gsfmodels = [ + simweights.GlobalSplineFit(), + simweights.GlobalSplineFit5Comp(), + simweights.GlobalSplineFit_IT(), +] + + +@pytest.mark.parametrize("gsf", gsfmodels) +def test_GlobalSplineFit_similar(gsf): """ - Test if GlobalSplineFit is similar to to other models within a factor of 500, + Test if GlobalSplineFit is similar to to other models within a factor of 0.4, mainly to check if the units provided by the GST files match. + Sum all species before check because different models use different particles. This can be not transparent in the file and web-interface. If the units mismatch, we should expect at least a deviation of 1000 (one prefix) or most likely a mismatch of 10 000 (cm^2 vs m^2). """ - gsf = _fluxes.GlobalSplineFit5Comp() - for name in ("GlobalFitGST", "GaisserH3a", "Hoerandel5"): - model = getattr(_fluxes, name)() - - same_pdgids = [pdgid for pdgid in gsf.pdgids if pdgid in model.pdgids] - - f_gsf = gsf(*np.meshgrid(E, [int(i) for i in same_pdgids])) - f = model(*np.meshgrid(E, [int(i) for i in same_pdgids])) - - assert 0.2 < np.mean(f / f_gsf) < 500 + model = getattr(simweights, name)() + f_gsf = gsf(*np.meshgrid(E, gsf.pdgids)).sum(axis=0) + f = model(*np.meshgrid(E, model.pdgids)).sum(axis=0) + assert f == pytest.approx(f_gsf, rel=0.4) if __name__ == "__main__": - unittest.main() + sys.exit(pytest.main(["-v", __file__, *sys.argv[1:]])) diff --git a/tests/test_fluxes/test_flux_model_FixedFractionFlux0_.npz b/tests/test_fluxes/test_flux_model_FixedFractionFlux0_.npz new file mode 100644 index 0000000..e8be2b4 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_FixedFractionFlux0_.npz differ diff --git a/tests/test_fluxes/test_flux_model_FixedFractionFlux1_.npz b/tests/test_fluxes/test_flux_model_FixedFractionFlux1_.npz new file mode 100644 index 0000000..e8be2b4 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_FixedFractionFlux1_.npz differ diff --git a/tests/test_fluxes/test_flux_model_GaisserH3a_.npz b/tests/test_fluxes/test_flux_model_GaisserH3a_.npz new file mode 100644 index 0000000..495eed6 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_GaisserH3a_.npz differ diff --git a/tests/test_fluxes/test_flux_model_GaisserH4a_.npz b/tests/test_fluxes/test_flux_model_GaisserH4a_.npz new file mode 100644 index 0000000..3ae5147 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_GaisserH4a_.npz differ diff --git a/tests/test_fluxes/test_flux_model_GaisserH4a_IT_.npz b/tests/test_fluxes/test_flux_model_GaisserH4a_IT_.npz new file mode 100644 index 0000000..3639bc7 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_GaisserH4a_IT_.npz differ diff --git a/tests/test_fluxes/test_flux_model_GaisserHillas_.npz b/tests/test_fluxes/test_flux_model_GaisserHillas_.npz new file mode 100644 index 0000000..958a182 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_GaisserHillas_.npz differ diff --git a/tests/test_fluxes/test_flux_model_GlobalFitGST_.npz b/tests/test_fluxes/test_flux_model_GlobalFitGST_.npz new file mode 100644 index 0000000..7e94810 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_GlobalFitGST_.npz differ diff --git a/tests/test_fluxes/test_flux_model_GlobalFitGST_IT_.npz b/tests/test_fluxes/test_flux_model_GlobalFitGST_IT_.npz new file mode 100644 index 0000000..4618962 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_GlobalFitGST_IT_.npz differ diff --git a/tests/test_fluxes/test_flux_model_GlobalSplineFit5Comp_.npz b/tests/test_fluxes/test_flux_model_GlobalSplineFit5Comp_.npz new file mode 100644 index 0000000..ab596b4 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_GlobalSplineFit5Comp_.npz differ diff --git a/tests/test_fluxes/test_flux_model_GlobalSplineFit_.npz b/tests/test_fluxes/test_flux_model_GlobalSplineFit_.npz new file mode 100644 index 0000000..fbe3310 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_GlobalSplineFit_.npz differ diff --git a/tests/test_fluxes/test_flux_model_GlobalSplineFit_IT_.npz b/tests/test_fluxes/test_flux_model_GlobalSplineFit_IT_.npz new file mode 100644 index 0000000..ffe04c1 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_GlobalSplineFit_IT_.npz differ diff --git a/tests/test_fluxes/test_flux_model_Hoerandel5_.npz b/tests/test_fluxes/test_flux_model_Hoerandel5_.npz new file mode 100644 index 0000000..029d652 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_Hoerandel5_.npz differ diff --git a/tests/test_fluxes/test_flux_model_Hoerandel_.npz b/tests/test_fluxes/test_flux_model_Hoerandel_.npz new file mode 100644 index 0000000..ed13c42 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_Hoerandel_.npz differ diff --git a/tests/test_fluxes/test_flux_model_Hoerandel_IT_.npz b/tests/test_fluxes/test_flux_model_Hoerandel_IT_.npz new file mode 100644 index 0000000..ebb94c5 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_Hoerandel_IT_.npz differ diff --git a/tests/test_fluxes/test_flux_model_Honda2004_.npz b/tests/test_fluxes/test_flux_model_Honda2004_.npz new file mode 100644 index 0000000..f4973a9 Binary files /dev/null and b/tests/test_fluxes/test_flux_model_Honda2004_.npz differ diff --git a/tests/test_fluxes/test_flux_model_TIG1996_.npz b/tests/test_fluxes/test_flux_model_TIG1996_.npz new file mode 100644 index 0000000..c989a4b Binary files /dev/null and b/tests/test_fluxes/test_flux_model_TIG1996_.npz differ diff --git a/tests/test_genie_datasets.py b/tests/test_genie_datasets.py index 56f51b5..ba74c64 100755 --- a/tests/test_genie_datasets.py +++ b/tests/test_genie_datasets.py @@ -5,81 +5,74 @@ # SPDX-License-Identifier: BSD-2-Clause import os -import unittest +import sys from pathlib import Path import h5py import numpy as np import pandas as pd +import pytest import tables import uproot from simweights import GenieWeighter from simweights._utils import get_column, get_table +datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None) +datasets = [ + "upgrade_genie_step3_140021_000000", + "upgrade_genie_step3_141828_000000", + "GENIE_NuMu_IceCubeUpgrade_v58.22590.000000", +] +approx = pytest.approx -class TestNugenDatasets(unittest.TestCase): - @classmethod - def setUpClass(cls): - datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None) - if not datadir: - cls.skipTest(None, "environment variable SIMWEIGHTS_TESTDATA not set") - cls.datadir = datadir + "/" - def cmp_dataset(self, fname): - filename = Path(self.datadir) / fname - reffile = h5py.File(str(filename) + ".hdf5", "r") - wd = reffile["I3MCWeightDict"] +@pytest.mark.parametrize("fname", datasets) +@pytest.mark.skipif(not datadir, reason="environment variable SIMWEIGHTS_TESTDATA not set") +def test_dataset(fname): + filename = Path(datadir) / fname + reffile = h5py.File(str(filename) + ".hdf5", "r") + wd = reffile["I3MCWeightDict"] - solid_angle = 2 * np.pi * (np.cos(wd["MinZenith"]) - np.cos(wd["MaxZenith"])) - injection_area = np.pi * (wd["InjectionSurfaceR"] * 1e2) ** 2 - global_probability_scale = wd["GlobalProbabilityScale"] - genie_weight = wd["GENIEWeight"] + solid_angle = 2 * np.pi * (np.cos(wd["MinZenith"]) - np.cos(wd["MaxZenith"])) + injection_area = np.pi * (wd["InjectionSurfaceR"] * 1e2) ** 2 + global_probability_scale = wd["GlobalProbabilityScale"] + genie_weight = wd["GENIEWeight"] - pli = -wd["PowerLawIndex"][0] - energy_integral = ((10 ** wd["MaxEnergyLog"][0]) ** (pli + 1) - (10 ** wd["MinEnergyLog"][0]) ** (pli + 1)) / (pli + 1) - energy_factor = 1 / (wd["PrimaryNeutrinoEnergy"] ** pli / energy_integral) - one_weight = global_probability_scale * genie_weight * energy_factor * solid_angle * injection_area - np.testing.assert_allclose(one_weight, wd["OneWeight"]) - final_weight = wd["OneWeight"] / (get_column(get_table(reffile, "I3GenieInfo"), "n_flux_events")[0]) + pli = -wd["PowerLawIndex"][0] + energy_integral = ((10 ** wd["MaxEnergyLog"][0]) ** (pli + 1) - (10 ** wd["MinEnergyLog"][0]) ** (pli + 1)) / (pli + 1) + energy_factor = 1 / (wd["PrimaryNeutrinoEnergy"] ** pli / energy_integral) + one_weight = global_probability_scale * genie_weight * energy_factor * solid_angle * injection_area + np.testing.assert_allclose(one_weight, wd["OneWeight"]) + final_weight = wd["OneWeight"] / (get_column(get_table(reffile, "I3GenieInfo"), "n_flux_events")[0]) - fobjs = [ - reffile, - uproot.open(str(filename) + ".root"), - tables.open_file(str(filename) + ".hdf5", "r"), - pd.HDFStore(str(filename) + ".hdf5", "r"), - ] + fobjs = [ + reffile, + uproot.open(str(filename) + ".root"), + tables.open_file(str(filename) + ".hdf5", "r"), + pd.HDFStore(str(filename) + ".hdf5", "r"), + ] - for fobj in fobjs: - with self.subTest(lib=str(fobj)): - w = GenieWeighter(fobj) + for fobj in fobjs: + w = GenieWeighter(fobj) - pdf0 = next(iter(w.surface.spectra.values()))[0].dists[0] - np.testing.assert_allclose(1 / pdf0.v, global_probability_scale * solid_angle * injection_area, 1e-5) + pdf0 = next(iter(w.surface.spectra.values()))[0].dists[0] + assert 1 / pdf0.v == approx(global_probability_scale * solid_angle * injection_area, rel=1e-5) - np.testing.assert_allclose(w.get_weight_column("wght"), genie_weight) + assert w.get_weight_column("wght") == approx(genie_weight) - power_law = next(iter(w.surface.spectra.values()))[0].dists[2] - energy_term = 1 / power_law.pdf(w.get_weight_column("energy")) - np.testing.assert_allclose(energy_term, energy_factor) + power_law = next(iter(w.surface.spectra.values()))[0].dists[2] + energy_term = 1 / power_law.pdf(w.get_weight_column("energy")) + assert energy_term == approx(energy_factor) - one_weight = w.get_weight_column("wght") * energy_term * solid_angle * injection_area * global_probability_scale - np.testing.assert_allclose(one_weight, wd["OneWeight"], 1e-5) + one_weight = w.get_weight_column("wght") * energy_term * solid_angle * injection_area * global_probability_scale + assert one_weight == approx(wd["OneWeight"], rel=1e-5) - np.testing.assert_allclose(w.get_weights(1), final_weight, 1e-5) + assert w.get_weights(1) == approx(final_weight, rel=1e-5) - for fobj in fobjs: - fobj.close() - - def test_140021(self): - self.cmp_dataset("upgrade_genie_step3_140021_000000") - - def test_141828(self): - self.cmp_dataset("upgrade_genie_step3_141828_000000") - - def test_22590(self): - self.cmp_dataset("GENIE_NuMu_IceCubeUpgrade_v58.22590.000000") + for fobj in fobjs: + fobj.close() if __name__ == "__main__": - unittest.main() + sys.exit(pytest.main(["-v", __file__, *sys.argv[1:]])) diff --git a/tests/test_icetop_datasets.py b/tests/test_icetop_datasets.py index 3248e5a..aabfaaf 100755 --- a/tests/test_icetop_datasets.py +++ b/tests/test_icetop_datasets.py @@ -5,63 +5,58 @@ # SPDX-License-Identifier: BSD-2-Clause import os -import unittest +import sys from pathlib import Path import h5py import numpy as np import pandas as pd +import pytest import tables import uproot from simweights import IceTopWeighter - -class TestIceTopDatasets(unittest.TestCase): - @classmethod - def setUpClass(cls): - datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None) - if not datadir: - cls.skipTest(None, "environment variable SIMWEIGHTS_TESTDATA not set") - cls.datadir = datadir + "/" - - def cmp_dataset(self, fname): - filename = Path(self.datadir) / fname - reffile = h5py.File(str(filename) + ".hdf5", "r") - - assert len(reffile["I3TopInjectorInfo"]) == 1 - si = reffile["I3TopInjectorInfo"][0] - pri = reffile["MCPrimary"] - solid_angle = np.pi * (np.cos(si["min_zenith"]) ** 2 - np.cos(si["max_zenith"]) ** 2) - injection_area = np.pi * (si["sampling_radius"] * 1e2) ** 2 - energy_integral = np.log(si["max_energy"] / si["min_energy"]) # assuming E^-1 - energy_factor = energy_integral * pri["energy"] - final_weight = energy_factor * solid_angle * injection_area / si["n_events"] - - fobjs = [ - reffile, - uproot.open(str(filename) + ".root"), - tables.open_file(str(filename) + ".hdf5", "r"), - pd.HDFStore(str(filename) + ".hdf5", "r"), - ] - - for fobj in fobjs: - with self.subTest(lib=str(fobj)): - w = IceTopWeighter(fobj) - spatial = w.surface.spectra[2212][0].dists[1] - proj_area = spatial.projected_area(1) - np.testing.assert_allclose(proj_area, injection_area) - sw_etendue = 1 / spatial.pdf(1) - np.testing.assert_allclose(sw_etendue, solid_angle * injection_area, 1e-5) - np.testing.assert_allclose(energy_integral * w.get_weight_column("energy"), energy_factor) - np.testing.assert_allclose(w.get_weights(1), final_weight, 1e-5) - - for fobj in fobjs: - fobj.close() - - def test_12360(self): - self.cmp_dataset("Level3_IC86.2012_SIBYLL2.1_p_12360_E6.0_0") +datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None) +datasets = ["Level3_IC86.2012_SIBYLL2.1_p_12360_E6.0_0"] +approx = pytest.approx + + +@pytest.mark.parametrize("fname", datasets) +@pytest.mark.skipif(not datadir, reason="environment variable SIMWEIGHTS_TESTDATA not set") +def test_dataset(fname): + filename = Path(datadir) / fname + reffile = h5py.File(str(filename) + ".hdf5", "r") + + assert len(reffile["I3TopInjectorInfo"]) == 1 + si = reffile["I3TopInjectorInfo"][0] + pri = reffile["MCPrimary"] + solid_angle = np.pi * (np.cos(si["min_zenith"]) ** 2 - np.cos(si["max_zenith"]) ** 2) + injection_area = np.pi * (si["sampling_radius"] * 1e2) ** 2 + energy_integral = np.log(si["max_energy"] / si["min_energy"]) # assuming E^-1 + energy_factor = energy_integral * pri["energy"] + final_weight = energy_factor * solid_angle * injection_area / si["n_events"] + + fobjs = [ + reffile, + uproot.open(str(filename) + ".root"), + tables.open_file(str(filename) + ".hdf5", "r"), + pd.HDFStore(str(filename) + ".hdf5", "r"), + ] + + for fobj in fobjs: + w = IceTopWeighter(fobj) + spatial = w.surface.spectra[2212][0].dists[1] + proj_area = spatial.projected_area(1) + assert proj_area == approx(injection_area) + sw_etendue = 1 / spatial.pdf(1) + assert sw_etendue == approx(solid_angle * injection_area, 1e-5) + assert energy_factor == approx(energy_integral * w.get_weight_column("energy")) + assert final_weight == approx(w.get_weights(1), 1e-5) + + for fobj in fobjs: + fobj.close() if __name__ == "__main__": - unittest.main() + sys.exit(pytest.main(["-v", __file__, *sys.argv[1:]])) diff --git a/tests/test_nugen_datasets.py b/tests/test_nugen_datasets.py index f65737b..b1e9a50 100755 --- a/tests/test_nugen_datasets.py +++ b/tests/test_nugen_datasets.py @@ -5,123 +5,94 @@ # SPDX-License-Identifier: BSD-2-Clause import os -import unittest +import sys from pathlib import Path import h5py import numpy as np import pandas as pd +import pytest import tables import uproot from simweights import NuGenWeighter - -class TestNugenDatasets(unittest.TestCase): - @classmethod - def setUpClass(cls): - datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None) - if not datadir: - cls.skipTest(None, "environment variable SIMWEIGHTS_TESTDATA not set") - cls.datadir = datadir + "/" - - def cmp_dataset(self, fname): - filename = Path(self.datadir) / fname - reffile = h5py.File(str(filename) + ".hdf5", "r") - - wd = reffile["I3MCWeightDict"] - pdgid = wd["PrimaryNeutrinoType"][0] - - solid_angle = 2 * np.pi * (np.cos(wd["MinZenith"]) - np.cos(wd["MaxZenith"])) - if "SolidAngle" in wd.dtype.names: - np.testing.assert_allclose(solid_angle, wd["SolidAngle"]) - - if "InjectionAreaCGS" in wd.dtype.names: - injection_area = wd["InjectionAreaCGS"] - if "InjectionAreaNormCGS" in wd.dtype.names: - injection_area = wd["InjectionAreaNormCGS"] - - if "TotalWeight" in wd.dtype.names: - total_weight = wd["TotalWeight"] - elif "TotalInteractionProbabilityWeight" in wd.dtype.names: - total_weight = wd["TotalInteractionProbabilityWeight"] - - type_weight = wd["TypeWeight"] if "TypeWeight" in wd.dtype.names else 0.5 - w0 = wd["OneWeight"] / (wd["NEvents"] * type_weight) - - fobjs = [ - reffile, - uproot.open(str(filename) + ".root"), - tables.open_file(str(filename) + ".hdf5", "r"), - pd.HDFStore(str(filename) + ".hdf5", "r"), - ] - - for fobj in fobjs: - with self.subTest(lib=str(fobj)): - w = NuGenWeighter(fobj, nfiles=1) - - event_weight = w.get_weight_column("event_weight") - np.testing.assert_allclose(event_weight, total_weight) - - cylinder = w.surface.spectra[pdgid][0].dists[2] - proj_area = cylinder.projected_area(w.get_weight_column("cos_zen")) - np.testing.assert_allclose(proj_area, injection_area) - - sw_etendue = 1 / cylinder.pdf(w.get_weight_column("cos_zen")) - np.testing.assert_allclose(sw_etendue, solid_angle * injection_area, 1e-5) - - power_law = w.surface.spectra[pdgid][0].dists[1] - energy_factor = 1 / power_law.pdf(w.get_weight_column("energy")) - one_weight = w.get_weight_column("event_weight") * energy_factor * solid_angle * injection_area - np.testing.assert_allclose(one_weight, wd["OneWeight"]) - - np.testing.assert_allclose(w.get_weights(1), w0, 1e-5) - - for fobj in fobjs: - fobj.close() - - def test_20885(self): - self.cmp_dataset("Level2_IC86.2016_NuE.020885.000000") - - def test_20878(self): - self.cmp_dataset("Level2_IC86.2016_NuMu.020878.000000") - - def test_20895(self): - self.cmp_dataset("Level2_IC86.2016_NuTau.020895.000000") - - def test_12646(self): - self.cmp_dataset("Level2_IC86.2012_nugen_nue.012646.000000.clsim-base-4.0.5.0.99_eff") - - def test_11836(self): - self.cmp_dataset("Level2_IC86.2012_nugen_nutau.011836.000000.clsim-base-4.0.3.0.99_eff") - - def test_11477(self): - self.cmp_dataset("Level2_IC86.2012_nugen_nutau.011477.000000.clsim-base-4.0.3.0.99_eff") - - def test_11374(self): - self.cmp_dataset("Level2_IC86.2012_nugen_numu.011374.000050.clsim-base-4.0.3.0.99_eff") - - def test_11297(self): - self.cmp_dataset("Level2_nugen_nutau_IC86.2012.011297.000000") - - def test_11070(self): - self.cmp_dataset("Level2_nugen_numu_IC86.2012.011070.000000") - - def test_11069(self): - self.cmp_dataset("Level2_nugen_numu_IC86.2012.011069.000000") - - def test_11065(self): - self.cmp_dataset("Level2_IC86.2012_nugen_NuTau.011065.000001") - - def test_11029(self): - self.cmp_dataset("Level2_nugen_numu_IC86.2012.011029.000000") - - def test_20692(self): - self.cmp_dataset("Level2_IC86.2011_nugen_NuE.010692.000000") - - def test_10634(self): - self.cmp_dataset("Level2_IC86.2011_nugen_NuMu.010634.000000") +datasets = [ + "Level2_IC86.2016_NuE.020885.000000", + "Level2_IC86.2016_NuMu.020878.000000", + "Level2_IC86.2016_NuTau.020895.000000", + "Level2_IC86.2012_nugen_nue.012646.000000.clsim-base-4.0.5.0.99_eff", + "Level2_IC86.2012_nugen_nutau.011836.000000.clsim-base-4.0.3.0.99_eff", + "Level2_IC86.2012_nugen_nutau.011477.000000.clsim-base-4.0.3.0.99_eff", + "Level2_IC86.2012_nugen_numu.011374.000050.clsim-base-4.0.3.0.99_eff", + "Level2_nugen_nutau_IC86.2012.011297.000000", + "Level2_nugen_numu_IC86.2012.011070.000000", + "Level2_nugen_numu_IC86.2012.011069.000000", + "Level2_IC86.2012_nugen_NuTau.011065.000001", + "Level2_nugen_numu_IC86.2012.011029.000000", + "Level2_IC86.2011_nugen_NuE.010692.000000", + "Level2_IC86.2011_nugen_NuMu.010634.000000", +] +approx = pytest.approx +datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None) + + +@pytest.mark.parametrize("fname", datasets) +@pytest.mark.skipif(not datadir, reason="environment variable SIMWEIGHTS_TESTDATA not set") +def test_dataset(fname): + filename = Path(datadir) / fname + reffile = h5py.File(str(filename) + ".hdf5", "r") + + wd = reffile["I3MCWeightDict"] + pdgid = wd["PrimaryNeutrinoType"][0] + + solid_angle = 2 * np.pi * (np.cos(wd["MinZenith"]) - np.cos(wd["MaxZenith"])) + if "SolidAngle" in wd.dtype.names: + np.testing.assert_allclose(solid_angle, wd["SolidAngle"]) + + if "InjectionAreaCGS" in wd.dtype.names: + injection_area = wd["InjectionAreaCGS"] + if "InjectionAreaNormCGS" in wd.dtype.names: + injection_area = wd["InjectionAreaNormCGS"] + + if "TotalWeight" in wd.dtype.names: + total_weight = wd["TotalWeight"] + elif "TotalInteractionProbabilityWeight" in wd.dtype.names: + total_weight = wd["TotalInteractionProbabilityWeight"] + + type_weight = wd["TypeWeight"] if "TypeWeight" in wd.dtype.names else 0.5 + w0 = wd["OneWeight"] / (wd["NEvents"] * type_weight) + + fobjs = [ + reffile, + uproot.open(str(filename) + ".root"), + tables.open_file(str(filename) + ".hdf5", "r"), + pd.HDFStore(str(filename) + ".hdf5", "r"), + ] + + for fobj in fobjs: + w = NuGenWeighter(fobj, nfiles=1) + + event_weight = w.get_weight_column("event_weight") + assert event_weight == approx(total_weight) + + cylinder = w.surface.spectra[pdgid][0].dists[2] + proj_area = cylinder.projected_area(w.get_weight_column("cos_zen")) + assert proj_area == approx(injection_area) + + sw_etendue = 1 / cylinder.pdf(w.get_weight_column("cos_zen")) + assert sw_etendue == approx(solid_angle * injection_area, rel=1e-5) + + power_law = w.surface.spectra[pdgid][0].dists[1] + energy_factor = 1 / power_law.pdf(w.get_weight_column("energy")) + one_weight = w.get_weight_column("event_weight") * energy_factor * solid_angle * injection_area + assert one_weight == approx(wd["OneWeight"]) + + assert w0 == approx(w.get_weights(1), rel=1e-5) + + for fobj in fobjs: + fobj.close() if __name__ == "__main__": - unittest.main() + sys.exit(pytest.main(["-v", __file__, *sys.argv[1:]])) diff --git a/tox.ini b/tox.ini deleted file mode 100644 index 6942f1a..0000000 --- a/tox.ini +++ /dev/null @@ -1,12 +0,0 @@ -; SPDX-FileCopyrightText: © 2022 the SimWeights contributors -; -; SPDX-License-Identifier: BSD-2-Clause - -[tox] -envlist = py3{8,9,10,11,12} -isolated_build = True - -[testenv] -passenv = SIMWEIGHTS_TESTDATA, HDF5_DIR -deps = .[test] -commands = pytest