From 2d15aea71e04b41021a7ea168e5ce19e82caa74e Mon Sep 17 00:00:00 2001 From: Pete Peterson Date: Wed, 7 Aug 2024 10:55:01 -0400 Subject: [PATCH] Multithread finding min/max and add timers --- .../src/AlignAndFocusPowderSlim.cpp | 130 ++++++++++++++---- 1 file changed, 102 insertions(+), 28 deletions(-) diff --git a/Framework/DataHandling/src/AlignAndFocusPowderSlim.cpp b/Framework/DataHandling/src/AlignAndFocusPowderSlim.cpp index 59e7020096b8..8bde87849352 100644 --- a/Framework/DataHandling/src/AlignAndFocusPowderSlim.cpp +++ b/Framework/DataHandling/src/AlignAndFocusPowderSlim.cpp @@ -18,6 +18,7 @@ #include "MantidKernel/VectorHelper.h" #include "MantidNexus/NexusIOHelper.h" #include "tbb/parallel_for.h" +#include "tbb/parallel_reduce.h" #include #include @@ -68,38 +69,107 @@ const std::vector AlignAndFocusPowderSlim::seeAlso() const { return //---------------------------------------------------------------------------------------------- namespace { // anonymous -class ProcessEventsTask { +class Histogrammer { public: - ProcessEventsTask(const std::vector *binedges, const double bin_divisor, const double bin_offset, - const std::vector *detids, const std::vector *tofs, - const AlignAndFocusPowderSlim::BankCalibration *calibration, - std::vector *y_temp) - : m_bin_divisor(bin_divisor), m_bin_offset(bin_offset), m_binedges(binedges), m_detids(detids), m_tofs(tofs), - m_calibration(calibration), y_temp(y_temp) {} + Histogrammer(const std::vector *binedges, const double width, const bool linear_bins) : m_binedges(binedges) { + m_xmin = binedges->front(); + m_xmax = binedges->back(); + + if (linear_bins) { + m_findBin = DataObjects::EventList::findLinearBin; + m_bin_divisor = 1. / width; + m_bin_offset = m_xmin * m_bin_divisor; + } else { + m_findBin = DataObjects::EventList::findLogBin; + m_bin_divisor = 1. / log1p(abs(width)); // use this to do change of base + m_bin_offset = log(m_xmin) * m_bin_divisor; + } + } + + boost::optional findBin(const double tof) const { + // return boost::none; + if (tof < m_xmin || tof >= m_xmax) { + return boost::none; + } else { + return m_findBin(*m_binedges, tof, m_bin_divisor, m_bin_offset, true); + } + } + +private: + double m_bin_divisor; + double m_bin_offset; + double m_xmin; + double m_xmax; + const std::vector *m_binedges; + boost::optional (*m_findBin)(const MantidVec &, const double, const double, const double, const bool); +}; + +template class ProcessEventsTask { +public: + ProcessEventsTask(const Histogrammer *histogrammer, const std::vector *detids, + const std::vector *tofs, const AlignAndFocusPowderSlim::BankCalibration *calibration, + std::vector *y_temp) + : m_histogrammer(histogrammer), m_detids(detids), m_tofs(tofs), m_calibration(calibration), y_temp(y_temp) {} void operator()(const tbb::blocked_range &range) const { for (size_t i = range.begin(); i < range.end(); ++i) { const auto detid = static_cast(m_detids->at(i)); const auto tof = static_cast(m_tofs->at(i)) * m_calibration->value(detid); - if (!(tof < m_binedges->front() || tof >= m_binedges->back())) { - const auto binnum = DataObjects::EventList::findLinearBin(*m_binedges, static_cast(tof), m_bin_divisor, - m_bin_offset, true); - if (binnum) - y_temp->at(binnum.get())++; - } + const auto binnum = m_histogrammer->findBin(tof); + if (binnum) + y_temp->at(binnum.get())++; } } private: - const double m_bin_divisor; - const double m_bin_offset; - const std::vector *m_binedges; + const Histogrammer *m_histogrammer; const std::vector *m_detids; const std::vector *m_tofs; const AlignAndFocusPowderSlim::BankCalibration *m_calibration; - std::vector *y_temp; + std::vector *y_temp; +}; + +template class MinMax { + const std::vector *vec; + +public: + Type minval; + Type maxval; + void operator()(const tbb::blocked_range &range) { + const auto [minele, maxele] = std::minmax_element(vec->cbegin() + range.begin(), vec->cbegin() + range.end()); + if (*minele < minval) + minval = *minele; + if (*maxele > maxval) + maxval = *maxele; + } + + MinMax(MinMax &other, tbb::split) + : vec(other.vec), minval(std::numeric_limits::max()), maxval(std::numeric_limits::min()) {} + + MinMax(const std::vector *vec) + : vec(vec), minval(std::numeric_limits::max()), maxval(std::numeric_limits::min()) {} + + void join(const MinMax &other) { + if (other.minval < minval) + minval = other.minval; + if (other.maxval > maxval) + maxval = other.maxval; + } }; + +template std::pair parallel_minmax(const std::vector *vec) { + constexpr size_t grainsize{2000}; + + if (vec->size() < grainsize) { + const auto [minval, maxval] = std::minmax_element(vec->cbegin(), vec->cend()); + return std::make_pair(*minval, *maxval); + } else { + MinMax finder(vec); + tbb::parallel_reduce(tbb::blocked_range(0, vec->size(), grainsize), finder); + return std::make_pair(finder.minval, finder.maxval); + } +} } // anonymous namespace //---------------------------------------------------------------------------------------------- @@ -130,7 +200,9 @@ void AlignAndFocusPowderSlim::exec() { constexpr double xmax{39950}; HistogramData::BinEdges XValues_new(0); - UNUSED_ARG(Kernel::VectorHelper::createAxisFromRebinParams({xmin, 10, xmax}, XValues_new.mutableRawData(), true, + const double binWidth{10.}; + const bool linearBins = bool(binWidth > 0.); + UNUSED_ARG(Kernel::VectorHelper::createAxisFromRebinParams({xmin, binWidth, xmax}, XValues_new.mutableRawData(), true, false, xmin, xmax)); const size_t numBins = XValues_new.size() - 1; MatrixWorkspace_sptr wksp = WorkspaceFactory::Instance().create("Workspace2D", numHist, numBins + 1, numBins); @@ -198,6 +270,7 @@ void AlignAndFocusPowderSlim::exec() { for (const std::string &classEntry : classEntries) { if (std::regex_match(classEntry, groups, classRegex)) { const std::string entry_name(groups[2].str()); + const auto startTimeBank = std::chrono::high_resolution_clock::now(); // skip entries with junk data if (entry_name == "bank_error_events" || entry_name == "bank_unmapped_events") @@ -217,24 +290,25 @@ void AlignAndFocusPowderSlim::exec() { addTimer("readDetID" + entry_name, startTime, std::chrono::high_resolution_clock::now()); } - const auto startTimeProcess = std::chrono::high_resolution_clock::now(); - const auto [minval, maxval] = std::minmax_element(event_detid->cbegin(), event_detid->cend()); - BankCalibration calibration(*minval, *maxval, m_calibration); - - auto binFinder = DataObjects::EventList::findLinearBin; - const auto divisor{.1}; - const auto offset{xmin * divisor}; + const auto startTimeSetup = std::chrono::high_resolution_clock::now(); + const auto [minval, maxval] = parallel_minmax(event_detid.get()); + BankCalibration calibration(static_cast(minval), static_cast(maxval), m_calibration); auto &spectrum = wksp->getSpectrum(specnum); - const auto &x_values = spectrum.readX(); + Histogrammer histogrammer(&spectrum.readX(), binWidth, linearBins); const auto numEvent = event_time_of_flight->size(); + // std::atomic allows for multi-threaded accumulation and who cares about floats when you are just + // counting things std::vector y_temp(spectrum.dataY().size()); - ProcessEventsTask task(&x_values, divisor, offset, event_detid.get(), event_time_of_flight.get(), &calibration, - &y_temp); + addTimer("setup" + entry_name, startTimeSetup, std::chrono::high_resolution_clock::now()); + + const auto startTimeProcess = std::chrono::high_resolution_clock::now(); + ProcessEventsTask task(&histogrammer, event_detid.get(), event_time_of_flight.get(), &calibration, &y_temp); tbb::parallel_for(tbb::blocked_range(0, numEvent), task); auto &y_values = spectrum.dataY(); std::copy(y_temp.cbegin(), y_temp.cend(), y_values.begin()); addTimer("proc" + entry_name, startTimeProcess, std::chrono::high_resolution_clock::now()); + addTimer(entry_name, startTimeBank, std::chrono::high_resolution_clock::now()); h5file.closeGroup(); specnum++;