BinHisto.h

Go to the documentation of this file.

#ifndef BIN_HISTO_H
#define BIN_HISTO_H
 
#include "BinningTools.h"  // For postSum computation
#include "Ippl.h"
 
#include <iomanip>  // for std::setw, std::setprecision, etc. (debug output)
 
namespace ParticleBinning {
 
    template <bool UseDualView, typename ViewType>
    struct DeviceViewTraits;
 
    template <typename ViewType>
    struct DeviceViewTraits<true, ViewType> {
        using h_type = typename ViewType::t_host;
        using d_type = typename ViewType::t_dev;
    };
 
    template <typename ViewType>
    struct DeviceViewTraits<false, ViewType> {
        using h_type = ViewType;
        using d_type = ViewType;
    };
 
    template <
            typename size_type, typename bin_index_type, typename value_type,
            bool UseDualView = false, class... Properties>
    class Histogram {
    public:
        // Histogram counts view type(s)
        using view_type = std::conditional_t<
                UseDualView, Kokkos::DualView<size_type*, Properties...>,
                Kokkos::View<size_type*, Properties...>>;
        using dview_type = typename DeviceViewTraits<UseDualView, view_type>::d_type;
        using hview_type = typename DeviceViewTraits<UseDualView, view_type>::h_type;
 
        // Histogram widths view type(s)
        using width_view_type = std::conditional_t<
                UseDualView, Kokkos::DualView<value_type*, Properties...>,
                Kokkos::View<value_type*, Properties...>>;
        using hwidth_view_type = typename DeviceViewTraits<UseDualView, width_view_type>::h_type;
        using dwidth_view_type = typename DeviceViewTraits<UseDualView, width_view_type>::d_type;
 
        // View types to map bin indices
        template <class... Args>
        using index_transform_type  = Kokkos::View<bin_index_type*, Args...>;
        using dindex_transform_type = index_transform_type<
                Kokkos::DefaultExecutionSpace>;  // Do not remove: needed inside
                                                 // AdaptBins::genAdaptiveHistogram()
        using hindex_transform_type = index_transform_type<Kokkos::HostSpace>;
        // using hash_type             =
        // ippl::detail::hash_type<Kokkos::DefaultExecutionSpace::memory_space>;
 
        Histogram() = default;
 
        Histogram(
                std::string debug_name, bin_index_type numBins, value_type totalBinWidth,
                value_type binningAlpha, value_type binningBeta, value_type desiredWidth)
            : debug_name_m(debug_name),
              numBins_m(numBins),
              totalBinWidth_m(totalBinWidth),
              binningAlpha_m(binningAlpha),
              binningBeta_m(binningBeta),
              desiredWidth_m(desiredWidth) {
            // Initialize the histogram, bin widths and post sum view
            instantiateHistograms();
            initTimers();
        }
 
        /*~Histogram() {
            //std::cout << "Histogram " << debug_name_m << " destroyed." << std::endl;
        } */
 
        Histogram(const Histogram& other) { copyFields(other); }
 
        Histogram& operator=(const Histogram& other) {
            if (this == &other) return *this;
            copyFields(other);
            return *this;
        }
 
        size_type getNPartInBin(bin_index_type binIndex) {
            if constexpr (UseDualView) {
                return histogram_m.view_host()(binIndex);
            } else if (std::is_same<typename hview_type::memory_space, Kokkos::HostSpace>::value) {
                return histogram_m(binIndex);
            } else {
                std::cerr << "Warning: Accessing BinHisto.getNPartInBin without DualView might be "
                             "inefficient!"
                          << std::endl;
                Kokkos::View<size_type, Kokkos::HostSpace> host_scalar("host_scalar");
                Kokkos::deep_copy(host_scalar, Kokkos::subview(histogram_m, binIndex));
                return host_scalar();
            }
        }
 
        size_type getCurrentBinCount() const { return numBins_m; }
 
        view_type getHistogram() { return histogram_m; }
        view_type getHistogram() const { return histogram_m; }
 
        view_type getPostSum() { return postSum_m; }
 
        width_view_type getBinWidths() const { return binWidths_m; }
 
        template <typename Histogram_t>
        void copyBinWidths(const Histogram_t& other) {
            using other_dwidth_view_type = typename Histogram_t::dwidth_view_type;
            Kokkos::deep_copy(
                    getDeviceView<dwidth_view_type>(binWidths_m),
                    other.template getDeviceView<other_dwidth_view_type>(other.getBinWidths()));
            if constexpr (UseDualView) {
                binWidths_m.modify_device();
 
                IpplTimings::startTimer(bDeviceSyncronizationT);
                binWidths_m.sync_host();
                IpplTimings::stopTimer(bDeviceSyncronizationT);
            }
        }
 
        void init() {  // const value_type constBinWidth
            // Assumes you have initialized histogram_m from the outside!
            sync();
            initConstBinWidths(totalBinWidth_m);
            initPostSum();
        }
 
        void initConstBinWidths(const value_type constBinWidth) {
            dwidth_view_type dWidthView = getDeviceView<dwidth_view_type>(binWidths_m);
            const value_type binWidth   = constBinWidth / numBins_m;
            using execution_space       = typename dwidth_view_type::execution_space;
 
            IpplTimings::startTimer(bHistogramInitT);
            // Note: "Kokkos::deep_copy(getDeviceView<dwidth_view_type>(binWidths_m), constBinWidth
            // / numBins_m);" would work too!
            Kokkos::parallel_for(
                    "InitConstBinWidths", Kokkos::RangePolicy<execution_space>(0, numBins_m),
                    KOKKOS_LAMBDA(const size_t i) { dWidthView(i) = binWidth; });
            IpplTimings::stopTimer(bHistogramInitT);
 
            if constexpr (UseDualView) {
                binWidths_m.modify_device();
                IpplTimings::startTimer(bDeviceSyncronizationT);
                binWidths_m.sync_host();
                IpplTimings::stopTimer(bDeviceSyncronizationT);
            }
        }
 
        void initPostSum() {
            IpplTimings::startTimer(bHistogramInitT);
            computeFixSum<dview_type>(
                    getDeviceView<dview_type>(histogram_m), getDeviceView<dview_type>(postSum_m));
            IpplTimings::stopTimer(bHistogramInitT);
 
            if constexpr (UseDualView) {
                postSum_m.modify_device();
                IpplTimings::startTimer(bDeviceSyncronizationT);
                postSum_m.sync_host();
                IpplTimings::stopTimer(bDeviceSyncronizationT);
            }
        }
 
        Kokkos::RangePolicy<> getBinIterationPolicy(
                const bin_index_type& binIndex1, const bin_index_type numBins = 1) {
            if constexpr (UseDualView) {
                return Kokkos::RangePolicy<>(
                        postSum_m.view_host()(binIndex1),
                        postSum_m.view_host()(binIndex1 + numBins));
            } else {
                std::cerr << "Warning: Accessing BinHisto.getBinIterationPolicy without DualView "
                             "might be inefficient!"
                          << std::endl;
                Kokkos::View<bin_index_type[2], Kokkos::HostSpace> host_ranges("host_scalar");
                Kokkos::deep_copy(
                        host_ranges,
                        Kokkos::subview(postSum_m, std::make_pair(binIndex1, binIndex1 + numBins)));
                return Kokkos::RangePolicy<>(host_ranges(0), host_ranges(1));
            }
        }
 
        /*
        Below are methods used for syncing the histogram view between host and device.
        If a normal View is used, they have no effect.
        Only necessary for the histogram view, since the binWidths and postSum views
        are only modified inside this class.
        */
 
        void sync() {
            IpplTimings::startTimer(bDeviceSyncronizationT);
            if constexpr (UseDualView) {
                if (histogram_m.need_sync_host() && histogram_m.need_sync_device()) {
                    std::cerr << "Warning: Histogram was modified on host AND device -- "
                                 "overwriting changes on host."
                              << std::endl;
                }
                if (histogram_m.need_sync_host()) {
                    histogram_m.sync_host();
                } else if (histogram_m.need_sync_device()) {
                    histogram_m.sync_device();
                }  // else do nothing
            }
            IpplTimings::stopTimer(bDeviceSyncronizationT);
        }
 
        void modify_device() {
            if constexpr (UseDualView) histogram_m.modify_device();
        }
 
        void modify_host() {
            if constexpr (UseDualView) histogram_m.modify_host();
        }
 
        template <typename return_type, typename HistogramType>
        static constexpr return_type getDeviceView(HistogramType histo) {
            if constexpr (UseDualView) {
                return histo.view_device();
            } else {
                return histo;
            }
        }
 
        template <typename return_type, typename HistogramType>
        static constexpr return_type getHostView(HistogramType histo) {
            if constexpr (UseDualView) {
                return histo.view_host();
            } else {
                return histo;
            }
        }
 
        hindex_transform_type mergeBins();
 
    private:
        value_type adaptiveBinningCostFunction(
                const size_type& sumCount, const value_type& sumWidth,
                const size_type& totalNumParticles);
 
        void initTimers() {
            bDeviceSyncronizationT = IpplTimings::getTimer("bDeviceSyncronization");
            bHistogramInitT        = IpplTimings::getTimer("bHistogramInit");
            bMergeBinsT            = IpplTimings::getTimer("bMergeBins");
        }
 
        void instantiateHistograms() {
            histogram_m = view_type("histogram", numBins_m);
            binWidths_m = width_view_type("binWidths", numBins_m);
            postSum_m   = view_type("postSum", numBins_m + 1);
        }
 
        void copyFields(const Histogram& other);
 
    private:
        std::string debug_name_m;    
        bin_index_type numBins_m;    
        value_type totalBinWidth_m;  
 
        value_type binningAlpha_m;  
        value_type binningBeta_m;  
        value_type desiredWidth_m;  
 
        view_type histogram_m;        
        width_view_type binWidths_m;  
        view_type postSum_m;  
 
        IpplTimings::TimerRef bDeviceSyncronizationT;
        IpplTimings::TimerRef bHistogramInitT;
        IpplTimings::TimerRef bMergeBinsT;
 
        /*
        Here are just some debug functions, like a nice histogram output formatted as python numpy
        arrays.
        */
    public:
        void printHistogram(std::ostream& os = std::cout) {
            if (ippl::Comm->rank() != 0) return;
            hview_type countsHost       = getHostView<hview_type>(histogram_m);
            hwidth_view_type widthsHost = getHostView<hwidth_view_type>(binWidths_m);
 
            // 3) Print header
            os << "Histogram \"" << debug_name_m << "\" with " << numBins_m
               << " bins. BinWidth = " << totalBinWidth_m << ".\n\n";
 
            // Format columns: BinIndex, Count, Width
            // Adjust widths as needed
            os << std::left << std::setw(10) << "Bin" << std::right << std::setw(12) << "Count"
               << std::setw(16) << "Width\n";
 
            os << std::string(38, '-') << "\n";
            // (38 dashes or however many you prefer to underline)
 
            // 4) Print each bin
            for (bin_index_type i = 0; i < numBins_m; ++i) {
                os << std::left << std::setw(10) << i  // bin index left-aligned
                   << std::right << std::setw(12) << countsHost(i) << std::fixed << std::setw(16)
                   << std::setprecision(6)
                   << static_cast<double>(widthsHost(i))  // in case 'value_type' is double/float
                   << "\n";
            }
 
            // os << "-----------------------------------------" << endl;
            os << std::endl;  // extra newline at the end
        }
 
        void printPythonArrays() const {
            if (ippl::Comm->rank() != 0) return;
            hview_type hostCounts       = getHostView<hview_type>(histogram_m);
            hwidth_view_type hostWidths = getHostView<hwidth_view_type>(binWidths_m);
            // TODO: if I leave this here, it may need a deep_copy to make it save for every
            // execution space
 
            Inform msg("AdaptBins::printPythonArrays");
 
            // Output counts as a Python NumPy array
            msg << level5 << "bin_counts = np.array([";
            for (bin_index_type i = 0; i < numBins_m; ++i) {
                msg << level5 << hostCounts(i);
                if (i < numBins_m - 1) msg << level5 << ", ";
            }
            msg << level5 << "])" << endl;
 
            // Output widths as a Python NumPy array
            msg << level5 << "bin_widths = np.array([";
            for (bin_index_type i = 0; i < numBins_m; ++i) {
                msg << level5 << std::fixed << std::setprecision(6) << hostWidths(i);
                if (i < numBins_m - 1) msg << level5 << ", ";
            }
            msg << level5 << "])" << endl;
        }
    };
 
}  // namespace ParticleBinning
 
#include "BinHisto.tpp"
 
#endif  // BIN_HISTO_H