Merge branch 'feature/max_files_option' of github.com:astrorama/SourceXtractorPlusPlus into feature/max_files_option

Author: marcschefer

.gitignore

@@ -11,3 +11,4 @@ __pycache__
 /packages
 /doc/build/
 /build/
+.vscode

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

SEBenchmarks/src/program/BenchVariablePsfStack.cpp

@@ -0,0 +1,136 @@
+/** 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 <map>
+#include <string>
+#include <random>
+
+#include <boost/program_options.hpp>
+#include <boost/timer/timer.hpp>
+#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<double> random_dist{0.0, 1000.0};
+
+public:
+
+  po::options_description defineSpecificProgramOptions() override {
+    po::options_description options{};
+    options.add_options()
+      ("iterations", po::value<int>()->default_value(100000), "Number of getPsf calls to benchmark")
+      ("measures", po::value<int>()->default_value(3), "Number of timing measurements to take")
+      ("fits-file", po::value<std::string>()->default_value(""), "FITS file containing PSF stack");
+    return options;
+  }
+
+  Elements::ExitCode mainMethod(std::map<std::string, po::variable_value> &args) override {
+
+    auto iterations = args["iterations"].as<int>();
+    auto measures = args["measures"].as<int>();
+    auto fits_file = args["fits-file"].as<std::string>();
+
+    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<CCfits::FITS> fitsPtr = nullptr;
+    if (!fits_file.empty()) {
+      try {
+        fitsPtr = std::make_shared<CCfits::FITS>(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();
+        return Elements::ExitCode::DATAERR;
+      }
+    } else {
+      logger.error() << "No FITS file provided";
+      return Elements::ExitCode::USAGE;
+    }
+    
+    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<std::vector<double>> 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();
+      return Elements::ExitCode::DATAERR;
+    }
+
+    return Elements::ExitCode::OK;
+  }
+};
+
+MAIN_FOR(BenchVariablePsfStack)

SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h

@@ -45,6 +45,30 @@ struct ImageCoordinate {
 
   ImageCoordinate() : m_x(0), m_y(0) {}
   ImageCoordinate(double x, double y) : m_x(x), m_y(y) {}
+
+  ImageCoordinate operator*(double scalar) const {
+    return ImageCoordinate(m_x * scalar, m_y * scalar);
+  }
+
+  ImageCoordinate operator+(const ImageCoordinate& other) const {
+    return ImageCoordinate(m_x + other.m_x, m_y + other.m_y);
+  }
+
+  ImageCoordinate& operator+=(const ImageCoordinate& other) {
+    m_x += other.m_x;
+    m_y += other.m_y;
+    return *this;
+  }
+
+  ImageCoordinate operator-(const ImageCoordinate& other) const {
+    return ImageCoordinate(m_x - other.m_x, m_y - other.m_y);
+  }
+
+  ImageCoordinate& operator-=(const ImageCoordinate& other) {
+    m_x -= other.m_x;
+    m_y -= other.m_y;
+    return *this;
+  }
 };
 
 class CoordinateSystem {

SEFramework/SEFramework/Output/OutputRegistry.h

@@ -169,7 +169,7 @@ class OutputRegistry {
         return converter(source.getProperty<PropertyType>(i));
       };
     }
-    Euclid::Table::Row::cell_type operator()(const SourceInterface& source) {
+    Euclid::Table::Row::cell_type operator()(const SourceInterface& source) const {
       return m_convert_func(source, index);
     }
     std::size_t index = 0;

SEFramework/SEFramework/Psf/VariablePsfStack.h

@@ -25,8 +25,10 @@
 #define _SEIMPLEMENTATION_PSF_VARIABLEPSFSTACK_H_
 
 #include <CCfits/CCfits>
+#include <memory>
 #include <SEFramework/Image/VectorImage.h>
 #include <SEFramework/Psf/Psf.h>
+#include <SEUtils/KdTree.h>
 
 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<SeFloat> m_ra_values;
-  std::vector<SeFloat> m_dec_values;
-  std::vector<SeFloat> m_x_values;
-  std::vector<SeFloat> m_y_values;
-  std::vector<double>  m_gridx_values;
-  std::vector<double>  m_gridy_values;
+  std::vector<PsfPosition> m_positions;
+  std::unique_ptr<KdTree<PsfPosition>> m_kdtree;
 
   std::vector<std::string> 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<VariablePsfStack::PsfPosition> {
+  static double getCoord(const VariablePsfStack::PsfPosition& pos, size_t index) {
+    return (index == 0) ? static_cast<double>(pos.x) : static_cast<double>(pos.y);
+  }
+};
+
 }  // namespace SourceXtractor
 
 #endif  //_SEIMPLEMENTATION_PSF_VARIABLEPSFSTACK_H_

SEFramework/src/lib/Output/OutputRegistry.cpp

@@ -49,21 +49,26 @@ auto OutputRegistry::getSourceToRowConverter(const std::vector<std::string>& ena
   return [this, out_prop_list](const SourceInterface& source) {
     std::vector<ColumnInfo::info_type> info_list {};
     std::vector<Row::cell_type> cell_values {};
+
+    const auto& property_to_names_map = m_property_to_names_map;
+    const auto& name_to_col_info_map = m_name_to_col_info_map;
+    const auto& name_to_converter_map = m_name_to_converter_map;
+
     for (const auto& property : out_prop_list) {
-      if (m_property_to_names_map.count(property) == 0) {
+      if (property_to_names_map.count(property) == 0) {
         throw Elements::Exception() << "Missing column generator for " << property.name();
       }
-      for (const auto& name : m_property_to_names_map.at(property)) {
-        auto& col_info = m_name_to_col_info_map.at(name);
-        info_list.emplace_back(name, m_name_to_converter_map.at(name).first,
-                               col_info.unit, col_info.description);
-        cell_values.emplace_back(m_name_to_converter_map.at(name).second(source));
+      for (const auto& name : property_to_names_map.at(property)) {
+        auto& col_info = name_to_col_info_map.at(name);
+        info_list.push_back(ColumnInfo::info_type {name, name_to_converter_map.at(name).first,
+                               col_info.unit, col_info.description});
+        cell_values.push_back(name_to_converter_map.at(name).second(source));
       }
     }
     if (info_list.empty()) {
       throw Elements::Exception() << "The given configuration would not generate any output";
     }
-    return Row {std::move(cell_values), std::make_shared<ColumnInfo>(move(info_list))};
+    return Row {std::move(cell_values), std::make_shared<ColumnInfo>(std::move(info_list))};
   };
 }
 

SEFramework/src/lib/Psf/VariablePsf.cpp

@@ -21,10 +21,14 @@
  *      Author: Alejandro Álvarez Ayllón
  */
 
-#include <ElementsKernel/Exception.h>
 #include <algorithm>
+#include <ElementsKernel/Logging.h>
+#include <ElementsKernel/Exception.h>
+#include "SEUtils/IsNan.h"
+
 #include "SEFramework/Psf/VariablePsf.h"
 
+static auto stack_logger = Elements::Logging::getLogger("PSFExPsf");
 
 namespace SourceXtractor {
 
@@ -130,6 +134,16 @@ void VariablePsf::selfTest() {
     if (coeff->getWidth() != psf_width || coeff->getHeight() != psf_height) {
       throw Elements::Exception() << "Malformed variable PSF, coefficient matrices do not have the same dimensions";
     }
+    for (auto x = 0; x < psf_width; ++x){
+      for (auto y = 0; y < psf_height; ++y) {
+        if (fastmath_isnan(coeff->at(x, y))) {
+          throw Elements::Exception() << "Malformed variable PSF, coefficient matrices contains NANs";
+        }
+        else if (fastmath_isinf(coeff->at(x, y))){
+          throw Elements::Exception() << "Malformed variable PSF, coefficient matrices contains INFs";
+        }
+      }
+    }
   }
 }