diff --git a/include/picongpu/plugins/binning/axis/Symlog.hpp b/include/picongpu/plugins/binning/axis/Symlog.hpp new file mode 100644 index 0000000000..fb96f0dcf3 --- /dev/null +++ b/include/picongpu/plugins/binning/axis/Symlog.hpp @@ -0,0 +1,263 @@ +/* Copyright 2025-2025 Tapish Narwal + * + * This file is part of PIConGPU. + * + * PIConGPU is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PIConGPU is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with PIConGPU. + * If not, see . + */ + +#pragma once + +#include "picongpu/plugins/binning/UnitConversion.hpp" +#include "picongpu/plugins/binning/axis/Axis.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace picongpu::plugins::binning::axis +{ + /** + * Symlog (symmetric log) axis with linear behavior near zero and logarithmic behavior for large values. + * Similar to matplotlib's symlog scale. + * The axis has a linear region within [-linearThreshold, linearThreshold] and logarithmic regions outside. + * Bins are distributed to maintain smooth transitions at the threshold boundaries. + * If overflow bins are enabled, allocates 2 extra bins for under and overflow. + */ + template + class SymlogAxis + { + private: + /** + * Transform value to symlog space + * For |val| <= threshold: returns val/threshold (normalized linear) + * For |val| > threshold: returns sign(val) * (1 + log(|val|/threshold)) + */ + static double symlogTransform(double val, double threshold) + { + if(std::abs(val) <= threshold) + { + return val / threshold; + } + else + { + return std::copysign(1.0 + std::log(std::abs(val) / threshold), val); + } + } + + /** + * Inverse transform from symlog space to real space + */ + static double symlogInverseTransform(double y, double threshold) + { + if(std::abs(y) <= 1.0) + { + return y * threshold; + } + else + { + return std::copysign(threshold * std::exp(std::abs(y) - 1.0), y); + } + } + + void initBinEdges() + { + binEdges.reserve(axisSplit.nBins + 1); + + // Transform min and max to symlog space + double yMin = symlogTransform(axisSplit.m_range.min, linearThresholdSI); + double yMax = symlogTransform(axisSplit.m_range.max, linearThresholdSI); + double yRange = yMax - yMin; + + // Create evenly spaced bins in transformed space + for(size_t i = 0; i <= axisSplit.nBins; i++) + { + double y = yMin + (i * yRange) / axisSplit.nBins; + double edge = symlogInverseTransform(y, linearThresholdSI); + binEdges.emplace_back(edge); + } + } + + public: + using Type = T_Attribute; + using ScalingType = std::conditional_t< + std::is_integral_v, + std::conditional_t, + T_Attribute>; + + AxisSplitting axisSplit; + /** Axis name, written out to OpenPMD */ + std::string label; + /** Units(Dimensionality) of the axis */ + std::array units; + /** Linear threshold in SI units */ + double linearThresholdSI; + /** axisSplit.nBins + 1 bin edges in SI units */ + std::vector binEdges; + + struct SymlogAxisKernel + { + /** Function to place particle on axis */ + T_AttrFunctor getAttributeValue; + /** Linear threshold in PIC units */ + ScalingType linearThresholdPIC; + /** 1/linearThresholdPIC */ + ScalingType invLinearThreshold; + /** Range in PIC units in symlog space */ + Range symlogPicRange; + // 1/yRange + ScalingType invYRange; + /** Enable or disable allocation of extra bins for out of range particles */ + bool overflowEnabled; + /** Number of bins excluding overflow */ + uint32_t nBins; + + constexpr SymlogAxisKernel( + T_AttrFunctor const& attrFunc, + AxisSplitting const& axisSplit, + double linearThreshold, + std::array const& unitsArr) + : getAttributeValue{attrFunc} + , linearThresholdPIC{static_cast(toPICUnits(linearThreshold, unitsArr))} + , invLinearThreshold{ScalingType(1) / linearThresholdPIC} + , symlogPicRange{symlogTransformKernel(toPICUnits(axisSplit.m_range.min, unitsArr)), symlogTransformKernel(toPICUnits(axisSplit.m_range.max, unitsArr))} + , invYRange{ ScalingType(1) / (symlogPicRange.max - symlogPicRange.min)} + , overflowEnabled{axisSplit.enableOverflowBins} + , nBins{axisSplit.nBins} + { + } + + /** + * Kernel-side symlog transform + */ + ALPAKA_FN_HOST_ACC ScalingType symlogTransformKernel(ScalingType val) const + { + ScalingType absVal = math::abs(val); + + // Compute both branches + ScalingType linear = val * invLinearThreshold; + ScalingType logarithmic = math::copysign(ScalingType(1) + math::log(absVal * invLinearThreshold), val); + + // Check if in linear region + ScalingType isLinear = ScalingType(absVal <= linearThresholdPIC); + + return isLinear ? linear : logarithmic; + } + + template + ALPAKA_FN_HOST_ACC std::pair getBinIdx(Args const&... args) const + { + auto val = getAttributeValue(args...); + auto sVal = static_cast(val); + ScalingType y = symlogTransformKernel(sVal); + + uint32_t inRange = (y >= symlogPicRange.min) & (y < symlogPicRange.max); + // uint32_t belowMin = (y < symlogPicRange.min); + uint32_t aboveMax = (y >= symlogPicRange.max); + + ScalingType normalized = (y - symlogPicRange.min) * invYRange; + + uint32_t rangeIdx = math::min(static_cast(normalized * nBins), nBins - 1); + + if(overflowEnabled) + { + // underflow*0 + inRange*(rangeIdx+1) + aboveMax*(nBins-1) + // I dont compute underflow/belowMin since it is always zero + uint32_t binIdx = inRange * (rangeIdx + 1) + aboveMax * (nBins + 1); + return {true, binIdx}; + } + else + { + return {inRange, rangeIdx}; + } + } + }; + + SymlogAxisKernel sAK; + + SymlogAxis( + AxisSplitting const& axSplit, + T_AttrFunctor const& attrFunctor, + double linearThreshold, + std::string const& label, + std::array const& units) + : axisSplit{axSplit} + , label{label} + , units{units} + , linearThresholdSI{linearThreshold} + , sAK{attrFunctor, axisSplit, linearThreshold, units} + { + if(linearThreshold <= 0) + { + throw std::domain_error("Linear threshold must be positive"); + } + initBinEdges(); + } + + constexpr uint32_t getNBins() const + { + return sAK.nBins; + } + + double getUnitConversion() const + { + return getConversionFactor(units); + } + + SymlogAxisKernel getAxisKernel() const + { + return sAK; + } + + /** + * @return bin edges in SI units + */ + std::vector getBinEdgesSI() const + { + return binEdges; + } + }; + + /** + * @details Creates a symlog (symmetric log) axis + * @tparam T_Attribute Type of the attribute + * @param axisSplitting Range and number of bins + * @param linearThreshold Threshold for linear region (in SI units) + * @param functorDescription Functor and metadata + */ + template + HINLINE auto createSymlog( + AxisSplitting const& axSplit, + double linearThreshold, + T_FunctorDescription const& functorDesc) + { + static_assert( + std::is_same_v, + "Access functor return type and range type should be the same"); + + return SymlogAxis( + axSplit, + functorDesc.functor, + linearThreshold, + functorDesc.name, + functorDesc.units); + } + +} // namespace picongpu::plugins::binning::axis diff --git a/include/pmacc/math/functions/Copysign.hpp b/include/pmacc/math/functions/Copysign.hpp new file mode 100644 index 0000000000..beefb3c119 --- /dev/null +++ b/include/pmacc/math/functions/Copysign.hpp @@ -0,0 +1,33 @@ +/* Copyright 2025-2025 Tapish Narwal + * + * This file is part of PMacc. + * + * PMacc is free software: you can redistribute it and/or modify + * it under the terms of either the GNU General Public License or + * the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PMacc is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License and the GNU Lesser General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License + * and the GNU Lesser General Public License along with PMacc. + * If not, see . + */ + + +#pragma once + +#include "pmacc/math/functions/Common.hpp" + +#include + +namespace pmacc::math +{ + //! Create value with the magnitude of x and the sign of y. + ALPAKA_BINARY_MATH_FN(copysign, alpaka::math::ConceptMathCopysign, Copysign) +} // namespace pmacc::math diff --git a/include/pmacc/math/math.hpp b/include/pmacc/math/math.hpp index 8252653605..9fc92ab4e0 100644 --- a/include/pmacc/math/math.hpp +++ b/include/pmacc/math/math.hpp @@ -24,6 +24,7 @@ #include "pmacc/math/functions/Abs.hpp" #include "pmacc/math/functions/Comparison.hpp" +#include "pmacc/math/functions/Copysign.hpp" #include "pmacc/math/functions/Erf.hpp" #include "pmacc/math/functions/Exp.hpp" #include "pmacc/math/functions/Log.hpp" diff --git a/share/picongpu/unit/non_dimensional_tests/binningAxis.cpp b/share/picongpu/unit/non_dimensional_tests/binningAxis.cpp index 483598ca73..4cd9d14c38 100644 --- a/share/picongpu/unit/non_dimensional_tests/binningAxis.cpp +++ b/share/picongpu/unit/non_dimensional_tests/binningAxis.cpp @@ -22,6 +22,7 @@ #include "picongpu/plugins/binning/FunctorDescription.hpp" #include "picongpu/plugins/binning/axis/LinearAxis.hpp" #include "picongpu/plugins/binning/axis/LogAxis.hpp" +#include "picongpu/plugins/binning/axis/Symlog.hpp" #include #include @@ -249,6 +250,73 @@ TEST_CASE("LogAxis Bin Edges with Units", "[axis][LogAxis][Edges][Units]") runAxisEdgesTestCases("Double", logVelocityTestCases, createLog, velocityFunctorDesc); } +TEST_CASE("SymlogAxis Bin Edges", "[axis][SymlogAxis][Edges]") +{ + auto createSymlog = [](auto const& split, auto threshold, auto const& funcDesc) + { return axis::createSymlog(split, threshold, funcDesc); }; + + SECTION("Double Type") + { + std::vector> const symlogDoubleTestCases + = {{"Symmetric range with threshold", + {{-100.0, 100.0}, 5, false}, + {-100.0, -18.47, -10.0, 10.0, 18.47, 100.0}, + false}, + // {"Asymmetric range positive heavy", + // {{-10.0, 1000.0}, 6, false}, + // {/* calculated edges based on transform */}, + // false}, + { + "Small range within threshold", + {{-5.0, 5.0}, 4, false}, + {-5.0, -2.5, 0.0, 2.5, 5.0}, // Should be linear + // false}, + // {"Zero crossing range", {{-50.0, 200.0}, 8, true}, {/* calculated edges */}, false + }}; + + auto functorDesc = createFunctorDescription(identityFunctor, "test"); + // Test cases with hardcoded 10 as linear region + for(auto const& testCase : symlogDoubleTestCases) + { + if(!testCase.expectThrow) + { + SECTION("Double: " + testCase.description) + { + auto axis = createSymlog(testCase.split, 10.0, functorDesc); + auto binEdges = axis.getBinEdgesSI(); + CHECK(binEdges.size() == testCase.split.nBins + 1); + CHECK(binEdges.front() == Catch::Approx(testCase.split.m_range.min)); + CHECK(binEdges.back() == Catch::Approx(testCase.split.m_range.max)); + + // Check monotonic increasing + for(size_t i = 1; i < binEdges.size(); ++i) + { + CHECK(binEdges[i] > binEdges[i - 1]); + } + } + } + } + } + + SECTION("Integer Type") + { + std::vector> const symlogIntTestCases + = {{"Symmetric range large", {{-1000, 1000}, 10, false}, {/* calculated edges */}, false}, + {"Positive only range", {{1, 1000}, 8, false}, {/* calculated edges */}, false}}; + + auto functorDesc = createFunctorDescription(identityFunctor, "test"); + for(auto const& testCase : symlogIntTestCases) + { + SECTION("Integer: " + testCase.description) + { + auto axis = createSymlog(testCase.split, 10, functorDesc); + auto binEdges = axis.getBinEdgesSI(); + CHECK(binEdges.size() == testCase.split.nBins + 1); + } + } + } +} + // --- Test Cases for Axis::getBinIdx --- template @@ -510,6 +578,107 @@ TEST_CASE("LogAxisKernel::getBinIdx", "[axis][LogAxis][Kernel]") } } +TEST_CASE("SymlogAxisKernel::getBinIdx", "[axis][SymlogAxis][Kernel]") +{ + SECTION("Double Type") + { + SECTION("Overflow Enabled, Symmetric Range [-100.0, 100.0], 10 bins, threshold 10.0") + { + axis::AxisSplitting split{{-100.0, 100.0}, 10, true}; + auto axis = axis::createSymlog(split, 10.0, createFunctorDescription(identityFunctor, "test")); + auto kernel = axis.getAxisKernel(); + + // Total bins = 10 + 2 = 12. Bins: 0 (under), 1-10 (in range), 11 (over) + // Transform range: yMin ~ -3.303, yMax ~ 3.303, yRange ~ 6.606 + + std::vector> tests = { + {"Value at zero", 0.0, {true, 6}}, // y=0 + {"Small positive in linear region", 5.0, {true, 6}}, // y=0.5 + {"Small negative in linear region", -5.0, {true, 5}}, // y=-0.5 + {"At positive threshold", 10.0, {true, 7}}, // y=1.0 + {"At negative threshold", -10.0, {true, 4}}, // y=-1.0 + {"Large positive value", 50.0, {true, 9}}, // y~2.609 + {"Large negative value", -50.0, {true, 2}}, // y~-2.609 + {"Value exactly max", 100.0, {true, 11}}, // At max -> overflow bin + {"Value exactly min", -100.0, {true, 1}}, // At min -> first regular bin + {"Value above max", 200.0, {true, 11}}, // Above max -> overflow bin + {"Value below min", -200.0, {true, 0}}, // Below min -> underflow bin + }; + runBinIdxTests(kernel, tests); + } + + SECTION("Overflow Disabled, Asymmetric Range [-50.0, 200.0], 8 bins, threshold 5.0") + { + axis::AxisSplitting split{{-50.0, 200.0}, 8, false}; + auto axis = axis::createSymlog(split, 5.0, createFunctorDescription(identityFunctor, "test")); + auto kernel = axis.getAxisKernel(); + + // Total bins = 8. Bins: 0-7 (in range) + // Transform range: yMin ~ -(1+ln(10)) ~ -3.303, yMax ~ 1+ln(40) ~ 4.689 + + std::vector> tests = { + {"Value at zero", 0.0, {true, 3}}, // y=0, in middle-left of range + {"Value in linear region", 2.5, {true, 3}}, // y=0.5, slightly right of zero + {"Value outside linear region positive", 20.0, {true, 5}}, // y~2.386 + {"Value outside linear region negative", -10.0, {true, 1}}, // y~-1.693 + {"Value exactly max", 200.0, {false, 0}}, // At max -> don't bin + {"Value exactly min", -50.0, {true, 0}}, // At min -> bin 0 + {"Value above max", 300.0, {false, 0}}, // Above max -> don't bin + {"Value below min", -100.0, {false, 0}}, // Below min -> don't bin + }; + runBinIdxTests(kernel, tests); + } + } + + SECTION("Integer Type") + { + SECTION("Overflow Enabled, Symmetric Range [-1000, 1000], 12 bins, threshold 10") + { + axis::AxisSplitting split{{-1000, 1000}, 12, true}; + auto axis = axis::createSymlog(split, 10, createFunctorDescription(identityFunctor, "test")); + auto kernel = axis.getAxisKernel(); + + // Total bins = 12 + 2 = 14. Bins: 0 (under), 1-12 (in range), 13 (over) + // Transform range: yMin ~ -(1+ln(100)) ~ -5.605, yMax ~ 5.605 + + std::vector> tests = { + {"Value at zero", 0, {true, 7}}, // y=0 -> bin 7 + {"Small positive in linear region", 5, {true, 7}}, // y=0.5, slightly right + {"Small negative in linear region", -5, {true, 6}}, // y=-0.5, slightly left + {"At positive threshold", 10, {true, 8}}, // y=1.0 + {"At negative threshold", -10, {true, 5}}, // y=-1.0 + {"Large positive value", 500, {true, 12}}, // y~4.912 + {"Large negative value", -300, {true, 2}}, // y~-4.689 + {"Value exactly max", 1000, {true, 13}}, // At max -> overflow bin + {"Value exactly min", -1000, {true, 1}}, // At min -> first regular bin + }; + runBinIdxTests(kernel, tests); + } + + SECTION("Overflow Disabled, Range [0, 1000], 10 bins, threshold 10") + { + axis::AxisSplitting split{{0, 1000}, 10, false}; + auto axis = axis::createSymlog(split, 10, createFunctorDescription(identityFunctor, "test")); + auto kernel = axis.getAxisKernel(); + + // Total bins = 10. Bins: 0-9 (in range) + // Transform range: yMin = 0, yMax ~ 5.605 + + std::vector> tests = { + {"Value at zero", 0, {true, 0}}, // y=0, at minimum + {"Small value in linear region", 5, {true, 0}}, // y=0.5 + {"At threshold", 10, {true, 1}}, // y=1.0 + {"Medium value", 100, {true, 5}}, // y~3.303 + {"Large value", 500, {true, 8}}, // y~4.912 + {"Value exactly max", 1000, {false, 0}}, // At max -> don't bin + {"Value below min", -10, {false, 0}}, // Below min -> don't bin + {"Value above max", 2000, {false, 0}}, // Above max -> don't bin + }; + runBinIdxTests(kernel, tests); + } + } +} + TEST_CASE("LinearAxisKernel::getBinIdx", "[axis][LinearAxis][Kernel]") { SECTION("Double Type")