From 9bf1525ab767bba6ea465c889eddf16b3fd70cf9 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 27 Jun 2025 17:24:43 +0200 Subject: [PATCH 1/3] Optimize VariablePsfStack using KdTree --- SEBenchmarks/CMakeLists.txt | 2 + .../src/program/BenchVariablePsfStack.cpp | 156 +++++++++++ .../SEFramework/Psf/VariablePsfStack.h | 39 ++- SEFramework/src/lib/Psf/VariablePsfStack.cpp | 88 +++--- SEUtils/CMakeLists.txt | 3 + SEUtils/SEUtils/KdTree.h | 2 + SEUtils/SEUtils/_impl/KdTree.icpp | 48 ++++ SEUtils/tests/src/KdTree_test.cpp | 258 ++++++++++++++++++ 8 files changed, 553 insertions(+), 43 deletions(-) create mode 100644 SEBenchmarks/src/program/BenchVariablePsfStack.cpp create mode 100644 SEUtils/tests/src/KdTree_test.cpp diff --git a/SEBenchmarks/CMakeLists.txt b/SEBenchmarks/CMakeLists.txt index 8c45f204c..0b5b7dd31 100644 --- a/SEBenchmarks/CMakeLists.txt +++ b/SEBenchmarks/CMakeLists.txt @@ -56,6 +56,8 @@ elements_add_executable(BenchRendering src/program/BenchRendering.cpp LINK_LIBRARIES SEFramework SEImplementation ${Boost_LIBRARIES}) elements_add_executable(BenchBackgroundModel src/program/BenchBackgroundModel.cpp LINK_LIBRARIES SEFramework SEImplementation ${Boost_LIBRARIES}) +elements_add_executable(BenchVariablePsfStack src/program/BenchVariablePsfStack.cpp + LINK_LIBRARIES SEFramework ${Boost_LIBRARIES}) #=============================================================================== # Declare the Boost tests here diff --git a/SEBenchmarks/src/program/BenchVariablePsfStack.cpp b/SEBenchmarks/src/program/BenchVariablePsfStack.cpp new file mode 100644 index 000000000..2e0c823e0 --- /dev/null +++ b/SEBenchmarks/src/program/BenchVariablePsfStack.cpp @@ -0,0 +1,156 @@ +/** Copyright © 2019-2025 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License as published by the Free + * Software Foundation; either version 3.0 of the License, or (at your option) + * any later version. + * + * This library 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 Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ + +/** + * @file src/program/BenchVariablePsfStack.cpp + * @date 06/27/25 + * @author marc schefer + */ + +#include +#include +#include + +#include +#include +#include "ElementsKernel/ProgramHeaders.h" +#include "ElementsKernel/Real.h" +#include "SEFramework/Psf/VariablePsfStack.h" + +namespace po = boost::program_options; +namespace timer = boost::timer; +using namespace SourceXtractor; + +static Elements::Logging logger = Elements::Logging::getLogger("BenchVariablePsfStack"); + +class BenchVariablePsfStack : public Elements::Program { +private: + std::default_random_engine random_generator; + std::uniform_real_distribution random_dist{0.0, 1000.0}; + +public: + + po::options_description defineSpecificProgramOptions() override { + po::options_description options{}; + options.add_options() + ("iterations", po::value()->default_value(100000), "Number of getPsf calls to benchmark") + ("measures", po::value()->default_value(10), "Number of timing measurements to take") + ("fits-file", po::value()->default_value(""), "FITS file containing PSF stack (optional, will use nullptr if not provided)"); + return options; + } + + Elements::ExitCode mainMethod(std::map &args) override { + + auto iterations = args["iterations"].as(); + auto measures = args["measures"].as(); + auto fits_file = args["fits-file"].as(); + + logger.info() << "Benchmarking VariablePsfStack::getPsf() with " << iterations << " iterations"; + logger.info() << "Taking " << measures << " timing measurements"; + + // Initialize VariablePsfStack with FITS file if provided, otherwise nullptr + std::shared_ptr fitsPtr = nullptr; + if (!fits_file.empty()) { + try { + fitsPtr = std::make_shared(fits_file); + logger.info() << "Using FITS file: " << fits_file; + } catch (const std::exception& e) { + logger.error() << "Failed to load FITS file '" << fits_file << "': " << e.what(); + logger.info() << "Continuing with nullptr - this will likely cause exceptions during getPsf() calls"; + } + } else { + logger.info() << "No FITS file provided - using nullptr (will likely cause exceptions)"; + } + + try { + VariablePsfStack psfStack(fitsPtr); + + logger.info() << "VariablePsfStack loaded successfully with " << psfStack.getNumberOfPsfs() << " PSFs"; + logger.info() << "PSF size: " << psfStack.getWidth() << "x" << psfStack.getHeight(); + logger.info() << "Pixel sampling: " << psfStack.getPixelSampling(); + + std::cout << "Iterations,Measurement,Time_nanoseconds" << std::endl; + + for (int m = 0; m < measures; ++m) { + logger.info() << "Measurement " << (m + 1) << "/" << measures; + + timer::cpu_timer timer; + timer.stop(); + + // Prepare test values for getPsf calls + std::vector> testValues; + testValues.reserve(iterations); + + for (int i = 0; i < iterations; ++i) { + testValues.push_back({random_dist(random_generator), random_dist(random_generator)}); + } + + // Start timing + timer.start(); + + for (int i = 0; i < iterations; ++i) { + try { + auto psf = psfStack.getPsf(testValues[i]); + // Prevent compiler optimization by using the result + volatile auto width = psf->getWidth(); + (void)width; // Suppress unused variable warning + } catch (const std::exception& e) { + // Expected to fail with nullptr, but we still measure the timing + // until the exception is thrown + } + } + + timer.stop(); + + auto elapsed_ns = timer.elapsed().wall; + std::cout << iterations << "," << (m + 1) << "," << elapsed_ns << std::endl; + + logger.info() << "Time for " << iterations << " calls: " << (elapsed_ns / 1e9) << " seconds"; + logger.info() << "Average time per call: " << (elapsed_ns / iterations) << " nanoseconds"; + } + + } catch (const std::exception& e) { + logger.error() << "Error initializing VariablePsfStack: " << e.what(); + logger.info() << "This is expected with nullptr parameter - fix the FITS file parameter later"; + + // Still run a basic timing test to measure overhead + std::cout << "Running basic timing test without actual PSF operations..." << std::endl; + std::cout << "Iterations,Measurement,Time_nanoseconds" << std::endl; + + for (int m = 0; m < measures; ++m) { + timer::cpu_timer timer; + timer.stop(); + + timer.start(); + for (int i = 0; i < iterations; ++i) { + // Just measure the overhead of generating random values and vector creation + std::vector values = {random_dist(random_generator), random_dist(random_generator)}; + volatile auto size = values.size(); + (void)size; + } + timer.stop(); + + auto elapsed_ns = timer.elapsed().wall; + std::cout << iterations << "," << (m + 1) << "," << elapsed_ns << std::endl; + } + } + + return Elements::ExitCode::OK; + } +}; + +MAIN_FOR(BenchVariablePsfStack) diff --git a/SEFramework/SEFramework/Psf/VariablePsfStack.h b/SEFramework/SEFramework/Psf/VariablePsfStack.h index 836b9cf11..c71fa0456 100644 --- a/SEFramework/SEFramework/Psf/VariablePsfStack.h +++ b/SEFramework/SEFramework/Psf/VariablePsfStack.h @@ -25,8 +25,10 @@ #define _SEIMPLEMENTATION_PSF_VARIABLEPSFSTACK_H_ #include +#include #include #include +#include namespace SourceXtractor { @@ -42,6 +44,18 @@ namespace SourceXtractor { */ class VariablePsfStack final : public Psf { public: + /** + * @brief Structure to hold PSF position data + */ + struct PsfPosition { + SeFloat ra; + SeFloat dec; + SeFloat x; + SeFloat y; + double gridx; + double gridy; + }; + /** * Constructor */ @@ -83,6 +97,13 @@ class VariablePsfStack final : public Psf { return m_components; }; + /** + * @return The number of PSFs loaded in the stack + */ + long getNumberOfPsfs() const { + return m_nrows; + }; + /** * */ @@ -99,12 +120,8 @@ class VariablePsfStack final : public Psf { long m_nrows; - std::vector m_ra_values; - std::vector m_dec_values; - std::vector m_x_values; - std::vector m_y_values; - std::vector m_gridx_values; - std::vector m_gridy_values; + std::vector m_positions; + std::unique_ptr> m_kdtree; std::vector m_components = {"X_IMAGE", "Y_IMAGE"}; @@ -119,6 +136,16 @@ class VariablePsfStack final : public Psf { void selfTest(); }; +/** + * @brief KdTree traits specialization for PsfPosition + */ +template <> +struct KdTreeTraits { + static double getCoord(const VariablePsfStack::PsfPosition& pos, size_t index) { + return (index == 0) ? static_cast(pos.x) : static_cast(pos.y); + } +}; + } // namespace SourceXtractor #endif //_SEIMPLEMENTATION_PSF_VARIABLEPSFSTACK_H_ diff --git a/SEFramework/src/lib/Psf/VariablePsfStack.cpp b/SEFramework/src/lib/Psf/VariablePsfStack.cpp index 2276e55d7..9b65da752 100644 --- a/SEFramework/src/lib/Psf/VariablePsfStack.cpp +++ b/SEFramework/src/lib/Psf/VariablePsfStack.cpp @@ -21,6 +21,7 @@ * Author: Martin Kuemmel */ #include +#include #include #include #include "SEFramework/Psf/VariablePsfStack.h" @@ -73,18 +74,33 @@ void VariablePsfStack::setup(std::shared_ptr pFits) { // read the nrows value m_nrows = position_data.rows(); + // Temporary vectors for reading data + std::vector ra_values, dec_values, x_values, y_values; + std::vector gridx_values, gridy_values; + try { // read in all the EXT specific columns - position_data.column("GRIDX", false).read(m_gridx_values, 0, m_nrows); - position_data.column("GRIDY", false).read(m_gridy_values, 0, m_nrows); + position_data.column("GRIDX", false).read(gridx_values, 0, m_nrows); + position_data.column("GRIDY", false).read(gridy_values, 0, m_nrows); } catch (CCfits::Table::NoSuchColumn) { - position_data.column("X_CENTER", false).read(m_gridx_values, 0, m_nrows); - position_data.column("Y_CENTER", false).read(m_gridy_values, 0, m_nrows); + position_data.column("X_CENTER", false).read(gridx_values, 0, m_nrows); + position_data.column("Y_CENTER", false).read(gridy_values, 0, m_nrows); + } + position_data.column("RA", false).read(ra_values, 0, m_nrows); + position_data.column("DEC", false).read(dec_values, 0, m_nrows); + position_data.column("X", false).read(x_values, 0, m_nrows); + position_data.column("Y", false).read(y_values, 0, m_nrows); + + // Populate the positions vector + m_positions.reserve(m_nrows); + for (long i = 0; i < m_nrows; ++i) { + m_positions.push_back({ra_values[i], dec_values[i], x_values[i], y_values[i], gridx_values[i], gridy_values[i]}); + } + + // Build KdTree for fast nearest neighbor searches + if (!m_positions.empty()) { + m_kdtree = std::make_unique>(m_positions); } - position_data.column("RA", false).read(m_ra_values, 0, m_nrows); - position_data.column("DEC", false).read(m_dec_values, 0, m_nrows); - position_data.column("X", false).read(m_x_values, 0, m_nrows); - position_data.column("Y", false).read(m_y_values, 0, m_nrows); } catch (CCfits::FitsException& e) { throw Elements::Exception() << "Error loading stacked PSF file: " << e.message(); @@ -95,54 +111,54 @@ void VariablePsfStack::selfTest() { int naxis1, naxis2; // read in the min/max grid values in x/y - const auto x_grid_minmax = std::minmax_element(begin(m_gridx_values), end(m_gridx_values)); - const auto y_grid_minmax = std::minmax_element(begin(m_gridy_values), end(m_gridy_values)); + auto x_grid_minmax = std::minmax_element(m_positions.begin(), m_positions.end(), + [](const PsfPosition& a, const PsfPosition& b) { return a.gridx < b.gridx; }); + auto y_grid_minmax = std::minmax_element(m_positions.begin(), m_positions.end(), + [](const PsfPosition& a, const PsfPosition& b) { return a.gridy < b.gridy; }); // read the image size m_pFits->extension(1).readKey("NAXIS1", naxis1); m_pFits->extension(1).readKey("NAXIS2", naxis2); // make sure all PSF in the grid are there - if (*x_grid_minmax.first - m_grid_offset < 1) - throw Elements::Exception() << "The PSF at the smallest x-grid starts at: " << *x_grid_minmax.first - m_grid_offset; - if (*y_grid_minmax.first - m_grid_offset < 1) - throw Elements::Exception() << "The PSF at the smallest y-grid starts at: " << *y_grid_minmax.first - m_grid_offset; - if (*x_grid_minmax.second + m_grid_offset > naxis1) - throw Elements::Exception() << "The PSF at the largest x-grid is too large: " << *x_grid_minmax.second + m_grid_offset + if (x_grid_minmax.first->gridx - m_grid_offset < 1) + throw Elements::Exception() << "The PSF at the smallest x-grid starts at: " << x_grid_minmax.first->gridx - m_grid_offset; + if (y_grid_minmax.first->gridy - m_grid_offset < 1) + throw Elements::Exception() << "The PSF at the smallest y-grid starts at: " << y_grid_minmax.first->gridy - m_grid_offset; + if (x_grid_minmax.second->gridx + m_grid_offset > naxis1) + throw Elements::Exception() << "The PSF at the largest x-grid is too large: " << x_grid_minmax.second->gridx + m_grid_offset << " NAXIS1: " << naxis1; - if (*y_grid_minmax.second + m_grid_offset > naxis2) - throw Elements::Exception() << "The PSF at the largest y-grid is too large: " << *y_grid_minmax.second + m_grid_offset - << " NAXIS2: " << naxis1; + if (y_grid_minmax.second->gridy + m_grid_offset > naxis2) + throw Elements::Exception() << "The PSF at the largest y-grid is too large: " << y_grid_minmax.second->gridy + m_grid_offset + << " NAXIS2: " << naxis2; } std::shared_ptr> VariablePsfStack::getPsf(const std::vector& values) const { - long index_min_distance = 0; - double min_distance = 1.0e+32; - // make sure there are only two positions if (values.size() != 2) throw Elements::Exception() << "There can be only two positional value for the stacked PSF!"; - // find the position of minimal distance - for (int act_index = 0; act_index < m_nrows; act_index++) { - double act_distance = (values[0] - m_x_values[act_index]) * (values[0] - m_x_values[act_index]) + - (values[1] - m_y_values[act_index]) * (values[1] - m_y_values[act_index]); - if (act_distance < min_distance) { - index_min_distance = act_index; - min_distance = act_distance; - } - } + // Use KdTree to find the nearest PSF position + KdTree::Coord coord; + coord.coord[0] = values[0]; // x coordinate + coord.coord[1] = values[1]; // y coordinate + + PsfPosition nearest_position = m_kdtree->findNearest(coord); + + // Calculate distance for logging + double min_distance = (values[0] - nearest_position.x) * (values[0] - nearest_position.x) + + (values[1] - nearest_position.y) * (values[1] - nearest_position.y); + // give some feedback stack_logger.debug() << "Distance: " << sqrt(min_distance) << " (" << values[0] << "," << values[1] << ")<-->(" - << m_x_values[index_min_distance] << "," << m_y_values[index_min_distance] - << ") index: " << index_min_distance; + << nearest_position.x << "," << nearest_position.y << ")"; // get the first and last pixels for the PSF to be extracted // NOTE: CCfits has 1-based indices, also the last index is *included* in the reading // NOTE: the +0.5 forces a correct cast/ceiling - std::vector first_vertex{long(m_gridx_values[index_min_distance]+.5) - long(m_grid_offset), long(m_gridy_values[index_min_distance]+.5) - long(m_grid_offset)}; + std::vector first_vertex{long(nearest_position.gridx+.5) - long(m_grid_offset), long(nearest_position.gridy+.5) - long(m_grid_offset)}; stack_logger.debug() << "First vertex: ( " << first_vertex[0] << ", " << first_vertex[1] << ") First vertex alternative: " << - m_gridx_values[index_min_distance]-m_grid_offset << " " << m_gridy_values[index_min_distance]-m_grid_offset << + nearest_position.gridx-m_grid_offset << " " << nearest_position.gridy-m_grid_offset << " grid offset:" << m_grid_offset; std::vector last_vertex{first_vertex[0] + long(m_psf_size) - 1, first_vertex[1] +long( m_psf_size) - 1}; @@ -155,8 +171,6 @@ std::shared_ptr> VariablePsfStack::getPsf(const std::vector m_pFits->extension(1).read(stamp_data, first_vertex, last_vertex, stride); } - //stack_logger.info() << "DDD ( " << first_vertex[0] << ", " << first_vertex[1] << ") --> ( " << last_vertex[0] << ", " << last_vertex[1] << "): " << stamp_data.size(); - // create and return the psf image return VectorImage::create(m_psf_size, m_psf_size, std::begin(stamp_data), std::end(stamp_data)); } diff --git a/SEUtils/CMakeLists.txt b/SEUtils/CMakeLists.txt index 9f3c453aa..64bcc6852 100644 --- a/SEUtils/CMakeLists.txt +++ b/SEUtils/CMakeLists.txt @@ -66,6 +66,9 @@ elements_add_unit_test(Misc_test tests/src/Misc_test.cpp elements_add_unit_test(QuadTree_test tests/src/QuadTree_test.cpp LINK_LIBRARIES SEUtils TYPE Boost) +elements_add_unit_test(KdTree_test tests/src/KdTree_test.cpp + LINK_LIBRARIES SEUtils + TYPE Boost) if(GMOCK_FOUND) elements_add_unit_test(Observable_test tests/src/Observable_test.cpp diff --git a/SEUtils/SEUtils/KdTree.h b/SEUtils/SEUtils/KdTree.h index 5c4fc3546..3e739c2eb 100644 --- a/SEUtils/SEUtils/KdTree.h +++ b/SEUtils/SEUtils/KdTree.h @@ -21,6 +21,7 @@ #include #include #include +#include namespace SourceXtractor { @@ -49,6 +50,7 @@ class KdTree { explicit KdTree(const std::vector& data); std::vector findPointsWithinRadius(Coord coord, double radius) const; + T findNearest(Coord coord) const; private: class Node; diff --git a/SEUtils/SEUtils/_impl/KdTree.icpp b/SEUtils/SEUtils/_impl/KdTree.icpp index 606b939cf..ce45e4d22 100644 --- a/SEUtils/SEUtils/_impl/KdTree.icpp +++ b/SEUtils/SEUtils/_impl/KdTree.icpp @@ -21,6 +21,7 @@ template class KdTree::Node { public: virtual std::vector findPointsWithinRadius(Coord coord, double radius) const = 0; + virtual T findNearest(Coord coord, double& best_dist_sq) const = 0; virtual ~Node() = default; }; @@ -45,6 +46,24 @@ public: return selection; } + virtual T findNearest(Coord coord, double& best_dist_sq) const { + T best_point = m_data[0]; + best_dist_sq = std::numeric_limits::max(); + + for (const auto& entry : m_data) { + double square_dist = 0.0; + for (size_t i = 0; i < N; i++) { + double delta = Traits::getCoord(entry, i) - coord.coord[i]; + square_dist += delta * delta; + } + if (square_dist < best_dist_sq) { + best_dist_sq = square_dist; + best_point = entry; + } + } + return best_point; + } + private: const std::vector m_data; }; @@ -101,6 +120,29 @@ public: } } + virtual T findNearest(Coord coord, double& best_dist_sq) const { + // Determine which side to search first + bool search_left_first = coord.coord[m_axis] < m_split_value; + auto first_child = search_left_first ? m_left_child : m_right_child; + auto second_child = search_left_first ? m_right_child : m_left_child; + + // Search the first side + T best_point = first_child->findNearest(coord, best_dist_sq); + + // Check if we need to search the other side + double axis_dist = coord.coord[m_axis] - m_split_value; + if (axis_dist * axis_dist < best_dist_sq) { + double second_best_dist_sq = best_dist_sq; + T second_best = second_child->findNearest(coord, second_best_dist_sq); + if (second_best_dist_sq < best_dist_sq) { + best_dist_sq = second_best_dist_sq; + best_point = second_best; + } + } + + return best_point; + } + private: size_t m_axis; double m_split_value; @@ -124,4 +166,10 @@ std::vector KdTree::findPointsWithinRadius(Coord coord, double radiu return m_root->findPointsWithinRadius(coord, radius); } +template +T KdTree::findNearest(Coord coord) const { + double best_dist_sq; + return m_root->findNearest(coord, best_dist_sq); +} + } diff --git a/SEUtils/tests/src/KdTree_test.cpp b/SEUtils/tests/src/KdTree_test.cpp new file mode 100644 index 000000000..5c5cf9e1b --- /dev/null +++ b/SEUtils/tests/src/KdTree_test.cpp @@ -0,0 +1,258 @@ +/** Copyright © 2021 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License as published by the Free + * Software Foundation; either version 3.0 of the License, or (at your option) + * any later version. + * + * This library 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 Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ + +#include +#include + +#include "SEUtils/KdTree.h" + +namespace SourceXtractor { + +struct TestPoint { + double m_x, m_y; + int m_id; // For identification in tests + + TestPoint(double x, double y, int id = 0) : m_x(x), m_y(y), m_id(id) {} +}; + +bool operator==(const TestPoint& a, const TestPoint& b) { + return a.m_x == b.m_x && a.m_y == b.m_y && a.m_id == b.m_id; +} + +template <> +struct KdTreeTraits { + static double getCoord(const TestPoint& t, size_t index) { + if (index == 0) { + return t.m_x; + } else { + return t.m_y; + } + } +}; + +// Test point for 3D tests +struct TestPoint3D { + double m_x, m_y, m_z; + int m_id; + + TestPoint3D(double x, double y, double z, int id = 0) : m_x(x), m_y(y), m_z(z), m_id(id) {} +}; + +bool operator==(const TestPoint3D& a, const TestPoint3D& b) { + return a.m_x == b.m_x && a.m_y == b.m_y && a.m_z == b.m_z && a.m_id == b.m_id; +} + +template <> +struct KdTreeTraits { + static double getCoord(const TestPoint3D& t, size_t index) { + if (index == 0) { + return t.m_x; + } else if (index == 1) { + return t.m_y; + } else { + return t.m_z; + } + } +}; + +BOOST_AUTO_TEST_SUITE (KdTree_test) + +BOOST_AUTO_TEST_CASE( smoke_test ) { + std::vector data = {{1.0, 2.0, 1}, {3.0, 4.0, 2}, {5.0, 6.0, 3}}; + KdTree tree(data); +} + +BOOST_AUTO_TEST_CASE( empty_tree_test ) { + std::vector data = {}; + // This should not crash - testing edge case + BOOST_CHECK_NO_THROW(KdTree tree(data)); +} + +BOOST_AUTO_TEST_CASE( single_point_test ) { + std::vector data = {{1.0, 2.0, 1}}; + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 1.0; + coord.coord[1] = 2.0; + + auto result = tree.findPointsWithinRadius(coord, 0.1); + BOOST_CHECK_EQUAL(result.size(), 1); + BOOST_CHECK(result[0] == TestPoint(1.0, 2.0, 1)); + + auto nearest = tree.findNearest(coord); + BOOST_CHECK(nearest == TestPoint(1.0, 2.0, 1)); +} + +BOOST_AUTO_TEST_CASE( simple_radius_search_test ) { + std::vector data = { + {1.0, 1.0, 1}, {2.0, 2.0, 2}, {3.0, 3.0, 3}, {10.0, 10.0, 4} + }; + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 2.0; + coord.coord[1] = 2.0; + + // Search with radius 2.0 - should find points at (1,1), (2,2), (3,3) + auto result = tree.findPointsWithinRadius(coord, 2.0); + BOOST_CHECK_EQUAL(result.size(), 3); + + // Search with radius 1.0 - should find only point at (2,2) + result = tree.findPointsWithinRadius(coord, 1.0); + BOOST_CHECK_EQUAL(result.size(), 1); + BOOST_CHECK(result[0] == TestPoint(2.0, 2.0, 2)); + + // Search with large radius - should find all points + result = tree.findPointsWithinRadius(coord, 15.0); + BOOST_CHECK_EQUAL(result.size(), 4); +} + +BOOST_AUTO_TEST_CASE( nearest_neighbor_test ) { + std::vector data = { + {1.0, 1.0, 1}, {2.0, 2.0, 2}, {3.0, 3.0, 3}, {10.0, 10.0, 4} + }; + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 2.1; + coord.coord[1] = 1.9; + + auto nearest = tree.findNearest(coord); + BOOST_CHECK(nearest == TestPoint(2.0, 2.0, 2)); + + // Test corner case - equidistant points + coord.coord[0] = 1.5; + coord.coord[1] = 1.5; + nearest = tree.findNearest(coord); + // Should return one of the two closest points (1,1) or (2,2) + BOOST_CHECK(nearest == TestPoint(1.0, 1.0, 1) || nearest == TestPoint(2.0, 2.0, 2)); +} + +BOOST_AUTO_TEST_CASE( large_dataset_test ) { + std::vector data; + + // Create a grid of 1000 points + for (int i = 0; i < 1000; i++) { + double x = 5.0 + (i % 100) / 10.0; + double y = 2.0 + (i / 100) / 5.0; + data.emplace_back(x, y, i); + } + + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 7.0; + coord.coord[1] = 3.0; + + // Test radius search on large dataset + auto result = tree.findPointsWithinRadius(coord, 1.0); + BOOST_CHECK_GT(result.size(), 0); + + // Test nearest neighbor on large dataset + auto nearest = tree.findNearest(coord); + BOOST_CHECK_GT(nearest.m_id, -1); // Should find a valid point +} + +BOOST_AUTO_TEST_CASE( three_dimensional_test ) { + std::vector data = { + {1.0, 1.0, 1.0, 1}, {2.0, 2.0, 2.0, 2}, {3.0, 3.0, 3.0, 3}, {10.0, 10.0, 10.0, 4} + }; + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 2.0; + coord.coord[1] = 2.0; + coord.coord[2] = 2.0; + + // Test radius search in 3D + auto result = tree.findPointsWithinRadius(coord, 2.0); + BOOST_CHECK_EQUAL(result.size(), 3); // Should find first 3 points + + // Test nearest neighbor in 3D + coord.coord[0] = 2.1; + coord.coord[1] = 1.9; + coord.coord[2] = 2.1; + + auto nearest = tree.findNearest(coord); + BOOST_CHECK(nearest == TestPoint3D(2.0, 2.0, 2.0, 2)); +} + +BOOST_AUTO_TEST_CASE( small_leaf_size_test ) { + std::vector data = { + {1.0, 1.0, 1}, {2.0, 2.0, 2}, {3.0, 3.0, 3}, {4.0, 4.0, 4}, {5.0, 5.0, 5} + }; + + // Use small leaf size to force tree structure + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 3.0; + coord.coord[1] = 3.0; + + auto result = tree.findPointsWithinRadius(coord, 1.5); + BOOST_CHECK_GE(result.size(), 1); + + auto nearest = tree.findNearest(coord); + BOOST_CHECK(nearest == TestPoint(3.0, 3.0, 3)); +} + +BOOST_AUTO_TEST_CASE( edge_case_identical_coordinates ) { + std::vector data = { + {1.0, 1.0, 1}, {1.0, 1.0, 2}, {1.0, 1.0, 3}, {2.0, 2.0, 4} + }; + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 1.0; + coord.coord[1] = 1.0; + + auto result = tree.findPointsWithinRadius(coord, 0.1); + BOOST_CHECK_EQUAL(result.size(), 3); // Should find all three identical points + + auto nearest = tree.findNearest(coord); + // Should return one of the three identical points + BOOST_CHECK(nearest.m_x == 1.0 && nearest.m_y == 1.0); +} + +BOOST_AUTO_TEST_CASE( stress_test_performance ) { + std::vector data; + + // Create 10000 random-ish points + for (int i = 0; i < 10000; i++) { + double x = (i * 17) % 1000 / 10.0; // Pseudo-random distribution + double y = (i * 23) % 1000 / 10.0; + data.emplace_back(x, y, i); + } + + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 50.0; + coord.coord[1] = 50.0; + + // This should complete reasonably quickly due to tree structure + auto result = tree.findPointsWithinRadius(coord, 10.0); + BOOST_CHECK_GE(result.size(), 0); + + auto nearest = tree.findNearest(coord); + BOOST_CHECK_GE(nearest.m_id, 0); +} + +BOOST_AUTO_TEST_SUITE_END () + +} From da6f69c6b7ddc99f5bfc15ad7038a74640f6a247 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Mon, 30 Jun 2025 11:42:47 +0200 Subject: [PATCH 2/3] benchmark code cleanup --- .../src/program/BenchVariablePsfStack.cpp | 32 ++++--------------- 1 file changed, 6 insertions(+), 26 deletions(-) diff --git a/SEBenchmarks/src/program/BenchVariablePsfStack.cpp b/SEBenchmarks/src/program/BenchVariablePsfStack.cpp index 2e0c823e0..e3de2eb86 100644 --- a/SEBenchmarks/src/program/BenchVariablePsfStack.cpp +++ b/SEBenchmarks/src/program/BenchVariablePsfStack.cpp @@ -48,8 +48,8 @@ class BenchVariablePsfStack : public Elements::Program { po::options_description options{}; options.add_options() ("iterations", po::value()->default_value(100000), "Number of getPsf calls to benchmark") - ("measures", po::value()->default_value(10), "Number of timing measurements to take") - ("fits-file", po::value()->default_value(""), "FITS file containing PSF stack (optional, will use nullptr if not provided)"); + ("measures", po::value()->default_value(3), "Number of timing measurements to take") + ("fits-file", po::value()->default_value(""), "FITS file containing PSF stack"); return options; } @@ -70,10 +70,11 @@ class BenchVariablePsfStack : public Elements::Program { logger.info() << "Using FITS file: " << fits_file; } catch (const std::exception& e) { logger.error() << "Failed to load FITS file '" << fits_file << "': " << e.what(); - logger.info() << "Continuing with nullptr - this will likely cause exceptions during getPsf() calls"; + return Elements::ExitCode::DATAERR; } } else { - logger.info() << "No FITS file provided - using nullptr (will likely cause exceptions)"; + logger.error() << "No FITS file provided"; + return Elements::ExitCode::USAGE; } try { @@ -125,28 +126,7 @@ class BenchVariablePsfStack : public Elements::Program { } catch (const std::exception& e) { logger.error() << "Error initializing VariablePsfStack: " << e.what(); - logger.info() << "This is expected with nullptr parameter - fix the FITS file parameter later"; - - // Still run a basic timing test to measure overhead - std::cout << "Running basic timing test without actual PSF operations..." << std::endl; - std::cout << "Iterations,Measurement,Time_nanoseconds" << std::endl; - - for (int m = 0; m < measures; ++m) { - timer::cpu_timer timer; - timer.stop(); - - timer.start(); - for (int i = 0; i < iterations; ++i) { - // Just measure the overhead of generating random values and vector creation - std::vector values = {random_dist(random_generator), random_dist(random_generator)}; - volatile auto size = values.size(); - (void)size; - } - timer.stop(); - - auto elapsed_ns = timer.elapsed().wall; - std::cout << iterations << "," << (m + 1) << "," << elapsed_ns << std::endl; - } + return Elements::ExitCode::DATAERR; } return Elements::ExitCode::OK; From 8f1e8a4b1859d2388de86427cb82934a2f4b1783 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Mon, 30 Jun 2025 11:44:36 +0200 Subject: [PATCH 3/3] update .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 8b9266d5a..7b4781cce 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,4 @@ __pycache__ /packages /doc/build/ /build/ +.vscode