Merge pull request #1442 from KrisThielemans/NoBoostIteratorAdaptor

Remove boost::iterator_adaptor from VectorWithOffset

src/include/stir/FullArrayIterator.h

@@ -22,6 +22,7 @@
 
 #include "stir/common.h"
 #include <iterator>
+#include "boost/iterator/iterator_adaptor.hpp"
 
 START_NAMESPACE_STIR
 
@@ -40,6 +41,7 @@ START_NAMESPACE_STIR
   class that uses FullArrayIterator implements end_all(). See the implementation of
   Array::end_all().
 
+  \todo use std::enable_if and std::is_convertible as opposed to boost::enable_if_convertible
   \internal
 */
 template <typename topleveliterT, typename restiterT, typename elemT, typename _Ref, typename _Ptr>

src/include/stir/VectorWithOffset.h

@@ -25,58 +25,10 @@
 
 #include "stir/shared_ptr.h"
 #include "stir/deprecated.h"
-#include "boost/iterator/iterator_adaptor.hpp"
-#include "boost/iterator/reverse_iterator.hpp"
+#include <iterator>
 
 START_NAMESPACE_STIR
 
-namespace detail
-{
-/*! \ingroup Array
-  \brief templated class for the iterators used by VectorWithOffset.
-
-  There should be no need to use this class yourself. Always use
-  VectorWithOffset::iterator or VectorWithOffset::const_iterator.
-*/
-template <class elemT>
-class VectorWithOffset_iter : public boost::iterator_adaptor<VectorWithOffset_iter<elemT> // Derived
-                                                             ,
-                                                             elemT* // Base
-                                                             ,
-                                                             boost::use_default // Value
-                                                             ,
-                                                             boost::random_access_traversal_tag // CategoryOrTraversal
-                                                             >
-{
-private:
-  // abbreviation of the type of this class
-  typedef VectorWithOffset_iter<elemT> self_t;
-
-public:
-  VectorWithOffset_iter()
-      : VectorWithOffset_iter::iterator_adaptor_(0)
-  {}
-
-  //! allow assignment from ordinary pointer
-  /*! really should be used within VectorWithOffset
-    It is explicit such that you can't do this by accident.
-  */
-  explicit VectorWithOffset_iter(elemT* p)
-      : VectorWithOffset_iter::iterator_adaptor_(p)
-  {}
-
-  //! some magic trickery to be able to assign iterators to const iterators, but not to incompatible types
-  /*! See the boost documentation for more info.
-   */
-  template <class OtherelemT>
-  VectorWithOffset_iter(VectorWithOffset_iter<OtherelemT> const& other,
-                        typename boost::enable_if_convertible<OtherelemT, elemT>::type* = 0)
-      : VectorWithOffset_iter::iterator_adaptor_(other.base())
-  {}
-};
-
-} // end of namespace detail
-
 /*!
   \ingroup Array
   \brief A templated class for vectors, but with indices starting not from 0
@@ -113,19 +65,16 @@ class VectorWithOffset
 {
 public:
   //! \name typedefs for iterator support
-  /*! Most of these should really not be needed because we use boost::iterator_adaptor now.
-      However, some are used directly in STIR code. (Maybe they shouldn't....)
-  */
   //@{
   typedef T value_type;
   typedef value_type& reference;
   typedef const value_type& const_reference;
   typedef ptrdiff_t difference_type;
-  typedef detail::VectorWithOffset_iter<T> iterator;
-  typedef detail::VectorWithOffset_iter<T const> const_iterator;
+  typedef T* iterator;
+  typedef T const* const_iterator;
 
-  typedef boost::reverse_iterator<iterator> reverse_iterator;
-  typedef boost::reverse_iterator<const_iterator> const_reverse_iterator;
+  typedef std::reverse_iterator<iterator> reverse_iterator;
+  typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
   //@}
   typedef size_t size_type;
 

src/include/stir/VectorWithOffset.inl

@@ -214,31 +214,31 @@ typename VectorWithOffset<T>::reverse_iterator
 VectorWithOffset<T>::rbegin()
 {
   this->check_state();
-  return boost::make_reverse_iterator(end());
+  return std::make_reverse_iterator(end());
 }
 
 template <class T>
 typename VectorWithOffset<T>::const_reverse_iterator
 VectorWithOffset<T>::rbegin() const
 {
   this->check_state();
-  return boost::make_reverse_iterator(end());
+  return std::make_reverse_iterator(end());
 }
 
 template <class T>
 typename VectorWithOffset<T>::reverse_iterator
 VectorWithOffset<T>::rend()
 {
   this->check_state();
-  return boost::make_reverse_iterator(begin());
+  return std::make_reverse_iterator(begin());
 }
 
 template <class T>
 typename VectorWithOffset<T>::const_reverse_iterator
 VectorWithOffset<T>::rend() const
 {
   this->check_state();
-  return boost::make_reverse_iterator(begin());
+  return std::make_reverse_iterator(begin());
 }
 
 template <class T>
@@ -567,10 +567,7 @@ void
 VectorWithOffset<T>::fill(const T& n)
 {
   this->check_state();
-  // TODO use std::fill() if we can use namespaces (to avoid name conflicts)
-  // std::fill(begin(), end(), n);
-  for (int i = this->get_min_index(); i <= this->get_max_index(); i++)
-    num[i] = n;
+  std::fill(this->begin(), this->end(), n);
   this->check_state();
 }
 

src/test/test_proj_data.cxx

@@ -111,23 +111,39 @@ ProjDataTests::run_tests_on_proj_data(ProjData& proj_data)
     // fill with values 0, 1, 2, ... to have something non-trivial
     std::iota(proj_data2.begin(), proj_data2.end(), 0.F);
     proj_data.fill(proj_data2);
+    const auto norm2 = proj_data2.norm();
     ProjDataInMemory proj_data3(proj_data2);
-    // check if equal by subtracting
-    proj_data3.sapyb(1.F, proj_data2, -1.F);
-    check(norm(proj_data3.begin(), proj_data3.end()) <= 0.0001F * norm(proj_data2.begin(), proj_data2.end()),
-          "fill(ProjDataInMemory&");
+    {
+      const auto norm3 = proj_data3.norm();
+      // check if equal by subtracting
+      proj_data3.sapyb(1.F, proj_data2, -1.F);
+      const auto normdiff = proj_data3.norm();
+      check(normdiff <= 0.0001F * norm2,
+            "fill(ProjDataInMemory&): norms org " + std::to_string(norm2) + " filled " + std::to_string(norm3) + " diff "
+                + std::to_string(normdiff));
+    }
     // test fill_from
-    fill_from(proj_data3, proj_data2.begin_all(), proj_data2.end_all());
-    // check if equal by subtracting
-    proj_data3.sapyb(1.F, proj_data2, -1.F);
-    check(norm(proj_data3.begin(), proj_data3.end()) <= 0.0001F * norm(proj_data2.begin(), proj_data2.end()),
-          "fill_from(ProjDataInMemory&");
+    {
+      fill_from(proj_data3, proj_data2.begin_all(), proj_data2.end_all());
+      const auto norm3 = proj_data3.norm();
+      // check if equal by subtracting
+      proj_data3.sapyb(1.F, proj_data2, -1.F);
+      const auto normdiff = proj_data3.norm();
+      check(normdiff <= 0.0001F * norm2,
+            "fill_from(ProjDataInMemory&): norms org " + std::to_string(norm2) + " filled " + std::to_string(norm3) + " diff "
+                + std::to_string(normdiff));
+    }
     // test copy_to
-    copy_to(proj_data2, proj_data3.begin_all());
-    // check if equal by subtracting
-    proj_data3.sapyb(1.F, proj_data2, -1.F);
-    check(norm(proj_data3.begin(), proj_data3.end()) <= 0.0001F * norm(proj_data2.begin(), proj_data2.end()),
-          "copy_to(ProjDataInMemory&");
+    {
+      copy_to(proj_data2, proj_data3.begin_all());
+      const auto norm3 = proj_data3.norm();
+      // check if equal by subtracting
+      proj_data3.sapyb(1.F, proj_data2, -1.F);
+      const auto normdiff = proj_data3.norm();
+      check(normdiff <= 0.0001F * norm2,
+            "copy_to(ProjDataInMemory&): norms org " + std::to_string(norm2) + " filled " + std::to_string(norm3) + " diff "
+                + std::to_string(normdiff));
+    }
   }
   timer.stop();
   std::cerr << "-- CPU Time " << timer.value() << '\n';

src/utilities/stir_timings.cxx

@@ -1,5 +1,5 @@
 /*
-    Copyright (C) 2023 University College London
+    Copyright (C) 2023, 2024 University College London
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0
@@ -50,9 +50,10 @@ static void
 print_usage_and_exit()
 {
   std::cerr << "\nUsage:\nstir_timings [--name some_string] [--threads num_threads] [--runs num_runs]\\\n"
-            << "\t[--skip-PP 1] [--skip-PMRT 1]\\\n"
+            << "\t[--skip-BB 1] [--skip-PP 1] [--skip-PMRT 1]\\\n"
             << "\t[--image image_filename]\\\n"
             << "\t--template-projdata template_proj_data_filename\n\n"
+            << "skip BB: basic building blocks; PP: Parallelproj; PMRT: ray-tracing matrix\n\n"
             << "Timings are reported to stdout as:\n"
             << "name\ttiming_name\tCPU_time_in_ms\twall-clock_time_in_ms\n";
   std::exit(EXIT_FAILURE);
@@ -68,13 +69,17 @@ class Timings : public TimedObject
   //! Use as prefix for all output
   std::string name;
   // variables that select timings
-  bool skip_PMRT;
-  bool skip_PP;
+  bool skip_BB;   //! skip basic building blocks
+  bool skip_PMRT; //! skip ProjMatrixByBinUsingRayTracing
+  bool skip_PP;   //! skip Parallelproj
 
   // variables used for running timings
-  shared_ptr<DiscretisedDensity<3, float>> image_sptr;
+  shared_ptr<VoxelsOnCartesianGrid<float>> image_sptr;
   shared_ptr<ProjData> output_proj_data_sptr;
   shared_ptr<ProjDataInMemory> mem_proj_data_sptr;
+  shared_ptr<ProjDataInMemory> mem_proj_data_sptr2;
+  std::vector<float> v1;
+  std::vector<float> v2;
   shared_ptr<ProjectorByBinPair> projectors_sptr;
   shared_ptr<ProjectorByBinPairUsingProjMatrixByBin> pmrt_projectors_sptr;
 #ifdef STIR_WITH_Parallelproj_PROJECTOR
@@ -88,7 +93,7 @@ class Timings : public TimedObject
   Timings(const std::string& image_filename, const std::string& template_proj_data_filename)
   {
     if (!image_filename.empty())
-      this->image_sptr = read_from_file<DiscretisedDensity<3, float>>(image_filename);
+      this->image_sptr = read_from_file<VoxelsOnCartesianGrid<float>>(image_filename);
 
     if (!template_proj_data_filename.empty())
       this->template_proj_data_sptr = ProjData::read_from_file(template_proj_data_filename);
@@ -107,12 +112,64 @@ class Timings : public TimedObject
     std::this_thread::sleep_for(std::chrono::milliseconds(1123));
   }
 
+  template <class T>
+  static void copy_add(T& t)
+  {
+    T c(t);
+    c += t;
+  }
+
+  template <class T>
+  static void copy_mult(T& t)
+  {
+    T c(t);
+    c *= t;
+  }
+
   void copy_image()
   {
     auto im = this->image_sptr->clone();
     delete im;
   }
 
+  void copy_add_image()
+  {
+    copy_add(*this->image_sptr);
+  }
+
+  void copy_mult_image()
+  {
+    copy_mult(*this->image_sptr);
+  }
+
+  void copy_std_vector()
+  {
+    std::copy(this->v1.begin(), this->v1.end(), this->v2.begin());
+  }
+  void create_std_vector()
+  {
+    std::vector<float> tmp(this->v1.size());
+    tmp[0] = 1; // assign something to avoid compiler warnings of unused variable
+  }
+  //! create proj_data in memory object
+  void create_proj_data_in_mem_no_init()
+  {
+    ProjDataInMemory tmp(this->template_proj_data_sptr->get_exam_info_sptr(),
+                         this->template_proj_data_sptr->get_proj_data_info_sptr(),
+                         /* initialise*/ false);
+  }
+  void create_proj_data_in_mem_init()
+  {
+    ProjDataInMemory tmp(this->template_proj_data_sptr->get_exam_info_sptr(),
+                         this->template_proj_data_sptr->get_proj_data_info_sptr(),
+                         /* initialise*/ true);
+  }
+  //! call ProjDataInMemory::fill(ProjDataInMemory&)
+  void copy_only_proj_data_mem_to_mem()
+  {
+    this->mem_proj_data_sptr2->fill(*this->mem_proj_data_sptr);
+  }
+
   //! copy from output_proj_data_sptr to new Interfile file
   void copy_proj_data_file_to_file()
   {
@@ -149,6 +206,16 @@ class Timings : public TimedObject
     tmp.fill(*this->mem_proj_data_sptr);
   }
 
+  void copy_add_proj_data_mem()
+  {
+    copy_add(*this->mem_proj_data_sptr);
+  }
+
+  void copy_mult_proj_data_mem()
+  {
+    copy_mult(*this->mem_proj_data_sptr);
+  }
+
   void projector_setup()
   {
     this->projectors_sptr->set_up(this->template_proj_data_sptr->get_proj_data_info_sptr(), this->image_sptr);
@@ -201,12 +268,32 @@ Timings::run_all(const unsigned runs)
 {
   this->init();
   // this->run_it(&Timings::sleep, "sleep", runs*1);
-  this->run_it(&Timings::copy_image, "copy_image", runs * 20);
   this->output_proj_data_sptr->fill(1.F);
-  this->run_it(&Timings::copy_proj_data_mem_to_mem, "copy_proj_data_mem_to_mem", runs * 2);
-  this->run_it(&Timings::copy_proj_data_mem_to_file, "copy_proj_data_mem_to_file", runs * 2);
-  this->run_it(&Timings::copy_proj_data_file_to_mem, "copy_proj_data_file_to_mem", runs * 2);
-  this->run_it(&Timings::copy_proj_data_file_to_file, "copy_proj_data_file_to_file", runs * 2);
+  if (!this->skip_BB)
+    {
+      this->mem_proj_data_sptr2
+          = std::make_shared<ProjDataInMemory>(this->exam_info_sptr, this->template_proj_data_sptr->get_proj_data_info_sptr());
+      this->v1.resize(this->template_proj_data_sptr->size_all());
+      this->v2.resize(this->template_proj_data_sptr->size_all());
+      this->run_it(&Timings::copy_image, "copy_image", runs * 20);
+      this->run_it(&Timings::copy_add_image, "copy_add_image", runs * 20);
+      this->run_it(&Timings::copy_mult_image, "copy_mult_image", runs * 20);
+      // reference timings: std::vector should be fast
+      this->run_it(&Timings::create_std_vector, "create_vector_of_size_projdata", runs * 2);
+      this->run_it(&Timings::copy_std_vector, "copy_std_vector_of_size_projdata", runs * 2);
+      v1.clear();
+      v2.clear();
+      this->run_it(&Timings::create_proj_data_in_mem_no_init, "create_proj_data_in_mem_no_init", runs * 2);
+      this->run_it(&Timings::create_proj_data_in_mem_init, "create_proj_data_in_mem_init", runs * 2);
+      this->run_it(&Timings::copy_only_proj_data_mem_to_mem, "copy_proj_data_mem_to_mem", runs * 2);
+      this->run_it(&Timings::copy_proj_data_mem_to_mem, "create_copy_proj_data_mem_to_mem", runs * 2);
+      this->mem_proj_data_sptr2.reset(); // no longer used
+      this->run_it(&Timings::copy_proj_data_mem_to_file, "create_copy_proj_data_mem_to_file", runs * 2);
+      this->run_it(&Timings::copy_proj_data_file_to_mem, "create_copy_proj_data_file_to_mem", runs * 2);
+      this->run_it(&Timings::copy_proj_data_file_to_file, "create_copy_proj_data_file_to_file", runs * 2);
+      this->run_it(&Timings::copy_add_proj_data_mem, "copy_add_proj_data_mem", runs * 2);
+      this->run_it(&Timings::copy_mult_proj_data_mem, "copy_mult_proj_data_mem", runs * 2);
+    }
   this->objective_function_sptr.reset(new PoissonLogLikelihoodWithLinearModelForMeanAndProjData<DiscretisedDensity<3, float>>);
   this->objective_function_sptr->set_proj_data_sptr(this->mem_proj_data_sptr);
   // this->objective_function.set_num_subsets(proj_data_sptr->get_num_views()/2);
@@ -321,6 +408,7 @@ main(int argc, char** argv)
   std::string prog_name = argv[0];
   unsigned num_runs = 3;
   int num_threads = get_default_num_threads();
+  bool skip_BB = false;
   bool skip_PMRT = false;
   bool skip_PP = false;
   // prefix output with this string
@@ -340,6 +428,8 @@ main(int argc, char** argv)
         num_runs = std::atoi(argv[1]);
       else if (!strcmp(argv[0], "--threads"))
         num_threads = std::atoi(argv[1]);
+      else if (!strcmp(argv[0], "--skip-BB"))
+        skip_BB = std::atoi(argv[1]) != 0;
       else if (!strcmp(argv[0], "--skip-PMRT"))
         skip_PMRT = std::atoi(argv[1]) != 0;
       else if (!strcmp(argv[0], "--skip-PP"))
@@ -358,6 +448,7 @@ main(int argc, char** argv)
 
   Timings timings(image_filename, template_proj_data_filename);
   timings.name = name;
+  timings.skip_BB = skip_BB;
   timings.skip_PMRT = skip_PMRT;
   timings.skip_PP = skip_PP;