MultiHistogram.hpp

Go to the documentation of this file.

// Copyright 2024 Sebastian Eibl
// Copyright 2026 Julian Friedrich Hille
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     https://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
 
#pragma once
 
#include <string>
 
#include "assert/assert.hpp"
#include "datatypes.hpp"
#include "functional.hpp"
 
namespace mrmd
{
namespace data
{
 
struct MultiHistogram
{
    MultiHistogram(const std::string& label,
                   const real_t minArg,
                   const real_t maxArg,
                   idx_t numBinsArg,
                   idx_t numHistogramsArg)
        : min(minArg),
          max(maxArg),
          numBins(numBinsArg),
          numHistograms(numHistogramsArg),
          binSize((max - min) / real_c(numBins)),
          inverseBinSize(1_r / binSize),
          data(label, numBins, numHistograms)
    {
        MRMD_HOST_CHECK_GREATER(maxArg, minArg);
        MRMD_HOST_CHECK_GREATEREQUAL(numBinsArg, 0);
        MRMD_HOST_CHECK_GREATEREQUAL(numHistogramsArg, 0);
    }
 
    MultiHistogram(const std::string& label, const MultiHistogram& histogram)
        : MultiHistogram(
              label, histogram.min, histogram.max, histogram.numBins, histogram.numHistograms)
    {
        Kokkos::deep_copy(data, histogram.data);
    }
 
    KOKKOS_INLINE_FUNCTION idx_t getBin(const real_t& val) const
    {
        auto bin = idx_c(std::floor((val - min) * inverseBinSize));
        if (bin < 0) bin = -1;
        if (bin >= numBins) bin = -1;
        return bin;
    }
 
    KOKKOS_INLINE_FUNCTION real_t getBinPosition(idx_t binIdx) const
    {
        assert(binIdx >= 0);
        assert(binIdx < numBins);
        auto binPosition = min + (real_c(binIdx) + 0.5_r) * binSize;
        return binPosition;
    }
 
    const real_t min;
    const real_t max;
    const idx_t numBins;
    const idx_t numHistograms;
    const real_t binSize;
    const real_t inverseBinSize;
    MultiView data;
 
    MultiHistogram& operator+=(const MultiHistogram& rhs);
    MultiHistogram& operator-=(const MultiHistogram& rhs);
    MultiHistogram& operator*=(const MultiHistogram& rhs);
    MultiHistogram& operator/=(const MultiHistogram& rhs);
 
    void scale(const real_t& scalingFactor);
    void scale(const ScalarView& scalingFactor);
    void makeSymmetric();
};
 
void cumulativeMovingAverage(data::MultiHistogram& average,
                             const data::MultiHistogram& current,
                             const real_t movingAverageFactor = 10_r);
 
data::MultiHistogram gradient(const data::MultiHistogram& input, const bool periodic = false);
 
data::MultiHistogram smoothen(data::MultiHistogram& input,
                              const real_t& sigma,
                              const real_t& range,
                              const bool periodic = false);
 
ScalarView createGrid(const MultiHistogram& input);
 
template <class BinaryOp>
void transform(const MultiHistogram& input1,
               const MultiHistogram& input2,
               MultiHistogram& output,
               const BinaryOp& binary_op)
{
    MRMD_HOST_CHECK_EQUAL(input1.numHistograms, input2.numHistograms);
    MRMD_HOST_CHECK_EQUAL(input1.numHistograms, output.numHistograms);
    MRMD_HOST_CHECK_EQUAL(input1.numBins, input2.numBins);
    MRMD_HOST_CHECK_EQUAL(input1.numBins, output.numBins);
 
    auto input1Data = input1.data;
    auto input2Data = input2.data;
    auto outputData = output.data;
 
    auto policy =
        Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {input1.numBins, input1.numHistograms});
    auto kernel = KOKKOS_LAMBDA(const idx_t idx, const idx_t jdx)
    {
        outputData(idx, jdx) = binary_op(input1Data(idx, jdx), input2Data(idx, jdx));
    };
    Kokkos::parallel_for("MultiHistogram::transform", policy, kernel);
    Kokkos::fence();
}
 
template <OneCoordinatePredicate Pred>
void replace_if_bin_position(MultiHistogram& hist, const Pred& pred, real_t newValue)
{
    auto policy =
        Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {hist.numBins, hist.numHistograms});
    auto kernel = KOKKOS_LAMBDA(const idx_t binIdx, const idx_t histIdx)
    {
        if (pred(hist.getBinPosition(binIdx)))
        {
            hist.data(binIdx, histIdx) = newValue;
        }
    };
    Kokkos::parallel_for("replace_if_bin_position", policy, kernel);
    Kokkos::fence();
}
 
}  // namespace data
}  // namespace mrmd