BinHisto.h

Go to the documentation of this file.

#ifndef BIN_HISTO_H
#define BIN_HISTO_H
 
#include "Ippl.h"
#include "BinningTools.h" // For postSum computation 
 
#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.h_view(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 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.h_view(binIndex1), postSum_m.h_view(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
 
            // Output counts as a Python NumPy array
            std::cout << "bin_counts = np.array([";
            for (bin_index_type i = 0; i < numBins_m; ++i) {
                std::cout << hostCounts(i);
                if (i < numBins_m - 1) std::cout << ", ";
            }
            std::cout << "])" << std::endl;
 
            // Output widths as a Python NumPy array
            std::cout << "bin_widths = np.array([";
            for (bin_index_type i = 0; i < numBins_m; ++i) {
                std::cout << std::fixed << std::setprecision(6) << hostWidths(i);
                if (i < numBins_m - 1) std::cout << ", ";
            }
            std::cout << "])" << std::endl;
        }
 
    };
    
}
 
#include "BinHisto.tpp"
 
#endif // BIN_HISTO_H