Add Radial and S2 Generators

include/integratorxx/generators/radial_factory.hpp

@@ -0,0 +1,78 @@
+#pragma once
+#include <integratorxx/quadratures/radial.hpp>
+
+namespace IntegratorXX {
+
+/// High-level specification of radial quadratures
+enum class RadialQuad : uint32_t {
+  Becke             = 0x0010, 
+  MurrayHandyLaming = 0x0020,
+  MuraKnowles       = 0x0030, 
+  TreutlerAhlrichs  = 0x0040
+};
+
+template <typename RadQuadType>
+RadialQuad radial_from_type() {
+  if constexpr (detail::is_becke_v<RadQuadType>) return RadialQuad::Becke;
+  if constexpr (detail::is_mk_v<RadQuadType>   ) return RadialQuad::MuraKnowles;
+  if constexpr (detail::is_mhl_v<RadQuadType>)   return RadialQuad::MurrayHandyLaming;
+  if constexpr (detail::is_ta_v<RadQuadType>)    return RadialQuad::TreutlerAhlrichs;
+
+  throw std::runtime_error("Unrecognized Radial Quadrature");
+};
+
+RadialQuad radial_from_string(std::string name);
+
+namespace detail {
+
+template <typename RadialTraitsType, typename... Args>
+std::unique_ptr<RadialTraits> make_radial_traits(Args&&... args) {
+  using traits_type = RadialTraitsType;
+  if constexpr (std::is_constructible_v<traits_type,Args...>)
+    return std::make_unique<traits_type>(std::forward<Args>(args)...);
+  else return nullptr;
+}
+
+}
+
+template <typename... Args>
+std::unique_ptr<RadialTraits> make_radial_traits(RadialQuad rq, Args&&... args) {
+  std::unique_ptr<RadialTraits> ptr;
+  switch(rq) {
+    case RadialQuad::Becke:
+      ptr = 
+        detail::make_radial_traits<BeckeRadialTraits>(std::forward<Args>(args)...);
+      break;
+    case RadialQuad::MurrayHandyLaming:
+      ptr = 
+        detail::make_radial_traits<MurrayHandyLamingRadialTraits<2>>(std::forward<Args>(args)...);
+      break;
+    case RadialQuad::MuraKnowles:
+      ptr = 
+        detail::make_radial_traits<MuraKnowlesRadialTraits>(std::forward<Args>(args)...);
+      break;
+    case RadialQuad::TreutlerAhlrichs:
+      ptr = 
+        detail::make_radial_traits<TreutlerAhlrichsRadialTraits>(std::forward<Args>(args)...);
+      break;
+  }
+
+  if(!ptr) throw std::runtime_error("RadialTraits Construction Failed");
+  return ptr;
+}
+
+
+struct RadialFactory {
+
+  using radial_grid_ptr = std::shared_ptr<
+    QuadratureBase<
+      std::vector<double>,
+      std::vector<double>
+    >
+  >;
+
+  static radial_grid_ptr generate(RadialQuad rq, const RadialTraits& traits);
+
+};
+
+}

include/integratorxx/generators/s2_factory.hpp

@@ -0,0 +1,41 @@
+#pragma once
+#include <integratorxx/quadratures/s2.hpp>
+
+#include <memory>
+
+namespace IntegratorXX {
+
+/// High-level specification of angular quadratures
+enum class AngularQuad : uint32_t {
+  AhrensBeylkin     = 0x0100,
+  Delley            = 0x0200,
+  LebedevLaikov     = 0x0300,
+  Womersley         = 0x0400
+};
+
+template <typename AngQuadType>
+AngularQuad angular_from_type() {
+  if constexpr (detail::is_ahrens_beyklin_v<AngQuadType>) return AngularQuad::AhrensBeylkin;
+  if constexpr (detail::is_delley_v<AngQuadType>   ) return AngularQuad::Delley;
+  if constexpr (detail::is_lebedev_laikov_v<AngQuadType>)   return AngularQuad::LebedevLaikov;
+  if constexpr (detail::is_womersley_v<AngQuadType>)    return AngularQuad::Womersley;
+
+  throw std::runtime_error("Unrecognized Angular Quadrature");
+};
+
+AngularQuad angular_from_string(std::string name);
+
+struct S2Factory {
+
+  using s2_grid_ptr = std::shared_ptr<
+    QuadratureBase<
+      std::vector<std::array<double,3>>,
+      std::vector<double>
+    >
+  >;
+
+  static s2_grid_ptr generate(AngularQuad aq, size_t npts);
+
+};
+
+}

include/integratorxx/generators/spherical_factory.hpp

@@ -2,53 +2,13 @@
 #include <memory>
 #include <vector>
 #include <array>
-#include <integratorxx/quadratures/radial.hpp>
-#include <integratorxx/quadratures/s2.hpp>
+#include <integratorxx/generators/radial_factory.hpp>
+#include <integratorxx/generators/s2_factory.hpp>
 #include <integratorxx/composite_quadratures/spherical_quadrature.hpp>
 #include <integratorxx/composite_quadratures/pruned_spherical_quadrature.hpp>
 
 namespace IntegratorXX {
 
-/// High-level specification of radial quadratures
-enum class RadialQuad : uint32_t {
-  Becke             = 0x0010, 
-  MurrayHandyLaming = 0x0020,
-  MuraKnowles       = 0x0030, 
-  TreutlerAhlrichs  = 0x0040
-};
-
-template <typename RadQuadType>
-RadialQuad radial_from_type() {
-  if constexpr (detail::is_becke_v<RadQuadType>) return RadialQuad::Becke;
-  if constexpr (detail::is_mk_v<RadQuadType>   ) return RadialQuad::MuraKnowles;
-  if constexpr (detail::is_mhl_v<RadQuadType>)   return RadialQuad::MurrayHandyLaming;
-  if constexpr (detail::is_ta_v<RadQuadType>)    return RadialQuad::TreutlerAhlrichs;
-
-  throw std::runtime_error("Unrecognized Radial Quadrature");
-};
-
-RadialQuad radial_from_string(std::string name);
-
-/// High-level specification of angular quadratures
-enum class AngularQuad : uint32_t {
-  AhrensBeylkin     = 0x0100,
-  Delley            = 0x0200,
-  LebedevLaikov     = 0x0300,
-  Womersley         = 0x0400
-};
-
-template <typename AngQuadType>
-AngularQuad angular_from_type() {
-  if constexpr (detail::is_ahrens_beyklin_v<AngQuadType>) return AngularQuad::AhrensBeylkin;
-  if constexpr (detail::is_delley_v<AngQuadType>   ) return AngularQuad::Delley;
-  if constexpr (detail::is_lebedev_laikov_v<AngQuadType>)   return AngularQuad::LebedevLaikov;
-  if constexpr (detail::is_womersley_v<AngQuadType>)    return AngularQuad::Womersley;
-
-  throw std::runtime_error("Unrecognized Angular Quadrature");
-};
-
-AngularQuad angular_from_string(std::string name);
-
 /// High-level specification of pruning schemes for spherical quadratures
 enum class PruningScheme {
   Unpruned, /// Unpruned quadrature
@@ -57,22 +17,29 @@ enum class PruningScheme {
 };
 
 // TODO: Make these strong (non-convertible) types
-using RadialScale = double;
-using RadialSize  = size_t;
+//using RadialScale = double;
+//using RadialSize  = size_t;
 using AngularSize = size_t;
-
+using radial_traits_ptr = std::unique_ptr<RadialTraits>;
 
 /// Generic specification of an unpruned spherical quadrature
 struct UnprunedSphericalGridSpecification {
-  RadialQuad  radial_quad;  ///< Radial quadrature specification
-  RadialSize  radial_size;  ///< Number of radial quadrature points
-  RadialScale radial_scale; ///< Radial scaling factor
+  RadialQuad   radial_quad;   ///< Radial quadrature specification
+  radial_traits_ptr radial_traits; ///< Radial traits (order, scaling factors, etc)
 
   AngularQuad angular_quad; /// Angular quadrature specification
   AngularSize angular_size; /// Number of angular quadrature points
+
+  UnprunedSphericalGridSpecification(RadialQuad, const RadialTraits&, AngularQuad,
+    AngularSize);
+
+  UnprunedSphericalGridSpecification(const UnprunedSphericalGridSpecification& other) :
+    radial_quad(other.radial_quad), radial_traits(other.radial_traits ? other.radial_traits->clone() : nullptr),
+    angular_quad(other.angular_quad), angular_size(other.angular_size) {};
 };
 
 
+
 /// Specification of a pruned region of an spherical quadrature
 struct PruningRegion {
   size_t idx_st; ///< Starting radial index for pruned region
@@ -91,15 +58,34 @@ struct PruningRegion {
 
 struct PrunedSphericalGridSpecification {
   RadialQuad  radial_quad;  ///< Radial quadrature specification
-  RadialSize  radial_size;  ///< Number of radial quadrature points
-  RadialScale radial_scale; ///< Radial scaling factor
+  radial_traits_ptr radial_traits; ///< Radial traits (order, scaling factors, etc)
 
   std::vector<PruningRegion> pruning_regions; ///< List of pruning regions over the radial quadrature
+
+  PrunedSphericalGridSpecification() = default;
+
+  template <typename... Arg>
+  PrunedSphericalGridSpecification(RadialQuad rq, radial_traits_ptr&& traits, Arg&&... arg) :
+    radial_quad(rq), radial_traits(std::move(traits)), pruning_regions(std::forward<Arg>(arg)...) { }
+  template <typename... Arg>
+  PrunedSphericalGridSpecification(RadialQuad rq, const RadialTraits& traits, Arg&&... arg) :
+    PrunedSphericalGridSpecification(rq, traits.clone(), std::forward<Arg>(arg)...) { }
+
+  PrunedSphericalGridSpecification(const PrunedSphericalGridSpecification& other) :
+    PrunedSphericalGridSpecification(other.radial_quad,
+      other.radial_traits ? other.radial_traits->clone() : nullptr,
+      other.pruning_regions) { }
+
+  PrunedSphericalGridSpecification& operator=(const PrunedSphericalGridSpecification& other) {
+    radial_quad = other.radial_quad;
+    radial_traits = other.radial_traits ? other.radial_traits->clone() : nullptr;
+    pruning_regions = other.pruning_regions;
+    return *this;
+  }
 
   inline bool operator==(const PrunedSphericalGridSpecification& other) const noexcept {
     return radial_quad == other.radial_quad and
-           radial_size == other.radial_size and
-           radial_scale == other.radial_scale and 
+           (radial_traits ? (other.radial_traits and radial_traits->compare(*other.radial_traits)) : !other.radial_traits) and 
            pruning_regions == other.pruning_regions;
   }
 };
@@ -176,10 +162,10 @@ struct SphericalGridFactory {
   }
 
 
-  static spherical_grid_ptr generate_unpruned_grid( RadialQuad, RadialSize, 
-    RadialScale, AngularQuad, AngularSize );
-  static spherical_grid_ptr generate_pruned_grid( RadialQuad, RadialSize, 
-    RadialScale, const std::vector<PruningRegion>&);
+  static spherical_grid_ptr generate_unpruned_grid( RadialQuad, const RadialTraits&, 
+    AngularQuad, AngularSize );
+  static spherical_grid_ptr generate_pruned_grid( RadialQuad, const RadialTraits&,
+    const std::vector<PruningRegion>&);
 
 
   static spherical_grid_ptr generate_grid(UnprunedSphericalGridSpecification gs);

include/integratorxx/quadratures/radial/becke.hpp

@@ -13,8 +13,10 @@ namespace IntegratorXX {
  *  J. Chem. Phys. 88, 2547 (1988)
  *  DOI: https://doi.org/10.1063/1.454033
  */
-class BeckeRadialTraits {
-  double R_;
+class BeckeRadialTraits : public RadialTraits {
+
+  size_t npts_; ///< Number of grid points
+  double R_; ///< Radial scaling factor
 
  public:
   /**
@@ -24,7 +26,22 @@ class BeckeRadialTraits {
    *
    *  @param[in] R     Radial scaling factor
    */
-  BeckeRadialTraits(double R = 1.0) : R_(R) {}
+  BeckeRadialTraits(size_t npts, double R = 1.0) : npts_(npts), R_(R) {}
+
+  size_t npts() const noexcept { return npts_; }
+
+  std::unique_ptr<RadialTraits> clone() const {
+    return std::make_unique<BeckeRadialTraits>(*this);
+  }
+
+  bool compare(const RadialTraits& other) const noexcept {
+    auto ptr = dynamic_cast<const BeckeRadialTraits*>(&other);
+    return ptr ? *this == *ptr : false;
+  }
+
+  bool operator==(const BeckeRadialTraits& other) const noexcept {
+    return npts_ == other.npts_ and R_ == other.R_;
+  }
 
   /**
    *  @brief Transformation rule for the Becke radial quadratures

include/integratorxx/quadratures/radial/mhl.hpp

@@ -17,13 +17,30 @@ namespace IntegratorXX {
  *            Typically taken to be 2.
  */
 template <size_t M>
-class MurrayHandyLamingRadialTraits {
+class MurrayHandyLamingRadialTraits : public RadialTraits {
 
+  size_t npts_; ///< Number of grid points
   double R_; ///< Radial scaling factor
+
 
 public:
 
-  MurrayHandyLamingRadialTraits(double R = 1.0) : R_(R) {}
+  MurrayHandyLamingRadialTraits(size_t npts, double R = 1.0) : npts_(npts), R_(R) {}
+
+  size_t npts() const noexcept { return npts_; }
+
+  std::unique_ptr<RadialTraits> clone() const {
+    return std::make_unique<MurrayHandyLamingRadialTraits>(*this);
+  }
+
+  bool compare(const RadialTraits& other) const noexcept {
+    auto ptr = dynamic_cast<const MurrayHandyLamingRadialTraits*>(&other);
+    return ptr ? *this == *ptr : false;
+  }
+
+  bool operator==(const MurrayHandyLamingRadialTraits& other) const noexcept {
+    return npts_ == other.npts_ and R_ == other.R_;
+  }
 
   /**
    *  @brief Transformation rule for the MHL radial quadrature

include/integratorxx/quadratures/radial/muraknowles.hpp

@@ -136,13 +136,29 @@ struct quadrature_traits<
 
 #else
 
-class MuraKnowlesRadialTraits {
+class MuraKnowlesRadialTraits : public RadialTraits {
 
+  size_t npts_; ///< Number of grid points
   double R_; ///< Radial scaling factor
 
 public:
 
-  MuraKnowlesRadialTraits(double R = 1.0) : R_(R) { }
+  MuraKnowlesRadialTraits(size_t npts, double R = 1.0) : npts_(npts), R_(R) { }
+
+  size_t npts() const noexcept { return npts_; }
+
+  std::unique_ptr<RadialTraits> clone() const {
+    return std::make_unique<MuraKnowlesRadialTraits>(*this);
+  }
+
+  bool compare(const RadialTraits& other) const noexcept {
+    auto ptr = dynamic_cast<const MuraKnowlesRadialTraits*>(&other);
+    return ptr ? *this == *ptr : false;
+  }
+
+  bool operator==(const MuraKnowlesRadialTraits& other) const noexcept {
+    return npts_ == other.npts_ and R_ == other.R_;
+  }
 
   template <typename PointType>
   inline auto radial_transform(PointType x) const noexcept {