geo_vectors_utils.h

Go to the documentation of this file.

#ifndef LARCOREOBJ_SIMPLETYPESANDCONSTANTS_GEO_VECTORS_UTILS_H
#define LARCOREOBJ_SIMPLETYPESANDCONSTANTS_GEO_VECTORS_UTILS_H

// LArSoft libraries
#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"

// C/C++ standard library
#include <array>
#include <vector>
#include <iterator> // std::back_inserter()
#include <type_traits> // std::declval(), std::is_same<>, ...
#include <functional> // std::mem_fn()
#include <cassert>


namespace geo {
  
  namespace vect {
    
    //--------------------------------------------------------------------------
    namespace details {
      //------------------------------------------------------------------------
      template <typename Op, typename... T>
      struct AccumulateImpl;
      
      template
        <typename Op, typename First, typename Second, typename... Others>
      struct AccumulateImpl<Op, First, Second, Others...> {
        static auto compute(Op op, First&& a, Second&& b, Others&&... others)
          -> decltype(auto) 
          {
            return op(
              a,
              AccumulateImpl<Op, Second, Others...>::compute
                (op, std::forward<Second>(b), std::forward<Others>(others)...)
              );
          }
      }; // AccumulateImpl<>
      
      template <typename Op, typename T>
      struct AccumulateImpl<Op, T> {
        static auto compute(Op, T&& v) -> decltype(auto)
          { return std::forward<T>(v); }
      };
      
      template <typename Op, typename... T>
      auto extended_accumulate(Op op, T&&... args)
        { return AccumulateImpl<Op, T...>::compute(op, std::forward<T>(args)...); }
      
      
      template <typename... T>
      auto extended_and(T... args) -> decltype(auto)
        {
          auto and_op = [](auto&& a, auto&& b) { return a && b; };
          return extended_accumulate(and_op, std::forward<T>(args)...); 
        }
      
      //------------------------------------------------------------------------
      //
      // These two pages of metaprogramming madness are aimed to have objects
      // that can provide a uniform coordinate get/set for any type of vector,
      // with the minimal run-time overhead.
      // Complications arise from the fact that TVector3 setter returns void,
      // while GenVector vectors return a reference to the vector itself
      // (hence the need of the additional SetterResult template parameter,
      // and of quite some code to autodetect its value).
      // The makeXxxx() functions are aimed to enclose the additional
      // autodetection overlay, which should become unnecessary with C++17.
      //
      
      template <typename Vector>
      struct VectorScalar { using type = typename Vector::Scalar; };
      
      template <typename Vector>
      using VectorScalar_t = typename VectorScalar<Vector>::type;
      
      
      //------------------------------------------------------------------------
      template <typename Vector>
      struct HasGetter {
          private:
        
        template<typename Class>
        static constexpr bool TestX(decltype(std::declval<Class>().X())*) { return true; }
        template<typename Class>
        static constexpr bool TestX(...) { return false; }
        template<typename Class>
        static constexpr bool TestY(decltype(std::declval<Class>().Y())*) { return true; }
        template<typename Class>
        static constexpr bool TestY(...) { return false; }
        template<typename Class>
        static constexpr bool TestZ(decltype(std::declval<Class>().Z())*) { return true; }
        template<typename Class>
        static constexpr bool TestZ(...) { return false; }
        template<typename Class>
        static constexpr bool TestT(decltype(std::declval<Class>().T())*) { return true; }
        template<typename Class>
        static constexpr bool TestT(...) { return false; }
        
          public:
        static constexpr bool X = TestX<Vector>(nullptr);
        static constexpr bool Y = TestY<Vector>(nullptr);
        static constexpr bool Z = TestZ<Vector>(nullptr);
        static constexpr bool T = TestT<Vector>(nullptr);
      }; // struct HasGetter<>

      template <typename Vector>
      constexpr bool HasX() { return HasGetter<Vector>::X; }
      template <typename Vector>
      constexpr bool HasY() { return HasGetter<Vector>::Y; }
      template <typename Vector>
      constexpr bool HasZ() { return HasGetter<Vector>::Z; }
      template <typename Vector>
      constexpr bool HasT() { return HasGetter<Vector>::T; }

      template <typename Vector>
      constexpr unsigned int dimension() {
        return HasT<Vector>()? 4U
          : HasZ<Vector>()? 3U
          : HasY<Vector>()? 2U
          : HasX<Vector>()? 1U
          : 0U;
      } // dimension()
      
      
      template <typename Vector>
      using CoordinateArray_t
        = std::array<VectorScalar_t<Vector>, dimension<Vector>()>;
      
      
      template <typename T>
      struct MemberFuncReturnType {
        // with C++17 `result_type` is deprecated,
        // and a different way will be needed
        using type
          = typename decltype(std::mem_fn(std::declval<T>()))::result_type;
      }; // MemberFuncReturnType
      
      template <typename T>
      using MemberFuncReturn_t = typename MemberFuncReturnType<T>::type;
      
      template <typename T>
      struct MemberFuncClassType;
      
      template <typename Class, typename Func>
      struct MemberFuncClassType<Func Class::*> { using type = Class; };
      
      template <typename T>
      using MemberFuncClass_t = typename MemberFuncClassType<T>::type;
      
      
      template <typename Vector, typename SetterType = void>
      struct BaseCoordTypes;
      
      template <typename Vector>
      struct BaseCoordTypes<Vector, void> {
        using Vector_t = std::decay_t<Vector>;
        using Scalar_t = VectorScalar_t<Vector>;
        using Getter_t = Scalar_t (Vector_t::*)() const;
      }; // struct BaseCoordTypes<void>
      
      template <typename Vector, typename SetterType>
      struct BaseCoordTypes {
          private:
        using BaseTypes_t = BaseCoordTypes<Vector, void>;
        
          public:
        using Scalar_t = typename BaseTypes_t::Scalar_t;
        using Vector_t = typename BaseTypes_t::Vector_t;
        using Getter_t = typename BaseTypes_t::Getter_t;
        using Setter_t = SetterType;
        using SetterReturn_t = MemberFuncReturn_t<Setter_t>;
        static_assert(
          std::is_same<Setter_t, SetterReturn_t (Vector_t::*)(Scalar_t)>(),
          "Invalid setter type"
          );
      }; // struct BaseCoordTypes<>
      
      
      template <typename Vector>
      struct CoordGetterTraits {
          private:
        using BaseTypes_t = BaseCoordTypes<Vector>;
          public:
        using Vector_t = typename BaseTypes_t::Vector_t;
        using Scalar_t = typename BaseTypes_t::Scalar_t;
        using Getter_t = typename BaseTypes_t::Getter_t;
      }; // struct CoordGetterTraits
      
      
      template <typename Vector>
      class CoordGetter {
        using Traits_t = CoordGetterTraits<Vector>;
        
          public:
        using Vector_t = typename Traits_t::Vector_t;
        using Scalar_t = typename Traits_t::Scalar_t;
        using Getter_t = typename Traits_t::Getter_t;
        
        constexpr CoordGetter(Getter_t getter): fGetter(getter) {}
        
        Scalar_t operator() (Vector_t const& v) const { return get(v); }
        
        Scalar_t get(Vector_t const& v) const { return (v.*fGetter)(); }
        
          private:
        Getter_t fGetter; 
        
      }; // class CoordGetter<>
      
      template <typename Getter>
      constexpr auto makeCoordReader(Getter getter)
        {
          using Vector_t = std::remove_reference_t<MemberFuncClass_t<Getter>>;
          return CoordGetter<Vector_t>{ getter };
        }
      
      
      template <typename Vector, typename SetterType>
      struct CoordManagerTraits {
          private:
        using BaseTypes_t = BaseCoordTypes<Vector, SetterType>;
          public:
        using Vector_t = typename BaseTypes_t::Vector_t;
        using Scalar_t = typename BaseTypes_t::Scalar_t;
        using Getter_t = typename BaseTypes_t::Getter_t;
        using Setter_t = typename BaseTypes_t::Setter_t;
      }; // struct VectorCoordManagerTraits<>
      
      
      template <typename Vector, typename SetterType>
      class CoordManager: public CoordGetter<Vector> {
        using Base_t = CoordGetter<Vector>;
        using Traits_t = CoordManagerTraits<Vector, SetterType>;
        
          public:
        using Vector_t = typename Traits_t::Vector_t; // this is not constant
        using Scalar_t = typename Traits_t::Scalar_t;
        using Getter_t = typename Traits_t::Getter_t;
        using Setter_t = typename Traits_t::Setter_t;
        
        constexpr CoordManager(Getter_t getter, Setter_t setter)
          : Base_t(getter), fSetter(setter) {}
        
        void operator()(Vector_t& v, Scalar_t c) const { set(v, c); }
        
        void set(Vector_t& v, Scalar_t c) const { (v.*fSetter)(c); }
        
        void incr(Vector_t& v, Scalar_t c) const { set(v, Base_t::get(v) + c); }
        
        void decr(Vector_t& v, Scalar_t c) const { set(v, Base_t::get(v) - c); }
        
        void mult(Vector_t& v, Scalar_t f) const { set(v, Base_t::get(v) * f); }
        
        void div(Vector_t& v, Scalar_t f) const { set(v, Base_t::get(v) / f); }
        
          private:
        Setter_t fSetter; 
      }; // class CoordManager<>
      
      
      template <typename Getter, typename Setter>
      constexpr auto makeCoordManager(Getter getter, Setter setter)
        {
          using Vector_t = std::remove_reference_t<MemberFuncClass_t<Getter>>;
          return CoordManager<Vector_t, Setter>{getter, setter};
        }
      
      
      
      template <typename CoordHelper, typename StoredVector>
      class BoundCoordGetter {
        
          public:
        using Stored_t = StoredVector;
        
        using CoordHelper_t = CoordHelper;
        using Vector_t = typename CoordHelper_t::Vector_t;
        using Scalar_t = typename CoordHelper_t::Scalar_t;
        using Getter_t = typename CoordHelper_t::Getter_t;
        
        BoundCoordGetter(Stored_t& v, CoordHelper_t coordManager)
          : fCoord(coordManager), fVector(v) {}
          
        BoundCoordGetter(Stored_t& v, Getter_t getter)
          : fCoord(getter), fVector(v) {}
        
        Scalar_t get() const { return manager().get(vector()); }
        
        Scalar_t operator() () const { return get(); }
        
        operator Scalar_t() const { return manager().get(vector()); }
        
          protected:
        CoordHelper_t const& manager() const { return fCoord; }
        Stored_t& vector() const { return fVector; }
        
          private:
        CoordHelper_t fCoord; 
        Stored_t& fVector; 
      }; // class VectorCoordGetter<>
      
      
      template <typename CoordHelper, typename StoredVector>
      class BoundCoordManager
        : public BoundCoordGetter<CoordHelper, StoredVector>
      {
        using Base_t = BoundCoordGetter<CoordHelper, StoredVector>;
        
          public:
        using typename Base_t::Stored_t;
        
        using CoordHelper_t = CoordHelper;
        using Vector_t = typename CoordHelper_t::Vector_t;
        using Scalar_t = typename CoordHelper_t::Scalar_t;
        using Getter_t = typename CoordHelper_t::Getter_t;
        using Setter_t = typename CoordHelper_t::Setter_t;
        
        BoundCoordManager(Stored_t& v, CoordHelper_t coordManager)
          : Base_t(v, coordManager) {}
          
        BoundCoordManager(Stored_t& v, Getter_t getter, Setter_t setter)
          : Base_t(v, CoordHelper_t(getter, setter)) {}
        
        BoundCoordManager& operator= (Scalar_t c)
          { Base_t::manager().set(Base_t::vector(), c); return *this; }
        
        BoundCoordManager& operator+= (Scalar_t c)
          { Base_t::manager().incr(Base_t::vector(), c); return *this; }
        
        BoundCoordManager& operator-= (Scalar_t c)
          { Base_t::manager().decr(Base_t::vector(), c); return *this; }
        
        BoundCoordManager& operator*= (Scalar_t f)
          { Base_t::manager().mult(Base_t::vector(), f); return *this; }
        
        BoundCoordManager& operator/= (Scalar_t f)
          { Base_t::manager().div(Base_t::vector(), f); return *this; }
        
      }; // class BoundCoordManager
      
      //------------------------------------------------------------------------
      
    } // namespace details
    
    
    // BEGIN Geometry group ------------------------------------------------------
  
    namespace extra {
      
      template <typename T>
      T roundValue01(T value, T tol) {
        if (std::abs(value) < tol) return 0.;
        if (std::abs(std::abs(value) - 1.) < tol) return (value > 0.)? 1.: -1.;
        return value;
      } // roundValue01()
      
    } // namespace extra
    
    
    // --- BEGIN Vector coordinate access abstraction --------------------------

    template <typename Vector>
    constexpr unsigned int dimension() { return details::dimension<Vector>(); }
    template <typename Vector>
    constexpr unsigned int dimension(Vector&&) { return dimension<Vector>(); }
    
    template <typename Vector>
    constexpr std::array<std::size_t, geo::vect::dimension<Vector>()> indices();
    template <typename Vector>
    constexpr auto indices(Vector const&) -> decltype(indices<Vector>());
    
    template <typename Vector>
    using coordinate_t = details::VectorScalar_t<Vector>;
    
    template <typename Vector, typename Coords>
    constexpr Vector makeFromCoords(Coords&& coords);
    
    template <typename Vector>
    using CoordReader_t
      = decltype(details::makeCoordReader(&Vector::X));
    
    template <typename Vector>
    using CoordManager_t = decltype
      (details::makeCoordManager(&Vector::X, &Vector::SetX));
    
    template <typename Vector>
    static constexpr auto XcoordManager
      = details::makeCoordManager(&Vector::X, &Vector::SetX);
    
    template <typename Vector>
    static constexpr auto XcoordManager<Vector const>
      = details::makeCoordReader(&Vector::X);
    
    template <typename Vector>
    static constexpr auto const YcoordManager
      = details::makeCoordManager(&Vector::Y, &Vector::SetY);
    
    template <typename Vector>
    static constexpr auto YcoordManager<Vector const>
      = details::makeCoordReader(&Vector::Y);
    
    template <typename Vector>
    static constexpr auto ZcoordManager
      = details::makeCoordManager(&Vector::Z, &Vector::SetZ);
    
    template <typename Vector>
    static constexpr auto ZcoordManager<Vector const>
      = details::makeCoordReader(&Vector::Z);
    
    template <typename Vector>
    static constexpr auto TcoordManager
      = details::makeCoordManager(&Vector::T, &Vector::SetT);
    
    template <typename Vector>
    static constexpr auto TcoordManager<Vector const>
      = details::makeCoordReader(&Vector::T);
    
    
    template <typename Vector>
    constexpr auto coordManager(unsigned int n);
    
    template <typename Vector>
    constexpr auto coordManager(unsigned int n, Vector& v);
    
    
    template <typename Vector>
    constexpr auto coordManagers();
    template <typename Vector>
    constexpr auto coordManagers(Vector&&);
    
    template <typename Vector>
    constexpr auto coordReaders();
    template <typename Vector>
    constexpr auto coordReaders(Vector&&);
    
    
    template <typename Vector>
    constexpr auto bindCoord
      (Vector const& v, CoordReader_t<Vector> helper)
      {
        return details::BoundCoordGetter<CoordReader_t<Vector>, Vector const>
          (v, helper);
      }
    
    template <typename Vector>
    auto bindCoord(Vector& v, CoordManager_t<Vector> helper)
      -> details::BoundCoordManager<CoordManager_t<Vector>, Vector>
      { return { v, helper }; }
    
    
    template <typename Vector>
    auto Xcoord(Vector& v)
      { return bindCoord(v, XcoordManager<Vector>); }
    
    template <typename Vector>
    auto Ycoord(Vector& v)
      { return bindCoord<Vector>(v, YcoordManager<Vector>); }
    
    template <typename Vector>
    auto Zcoord(Vector& v)
      { return bindCoord<Vector>(v, ZcoordManager<Vector>); }
    
    template <typename Vector>
    auto Tcoord(Vector& v)
      { return bindCoord<Vector>(v, TcoordManager<Vector>); }
    
    template <typename Vector>
    auto coord(Vector& v, unsigned int n) noexcept;
    
    
    template <typename Vector>
    constexpr auto bindCoordManagers(Vector& v);
    
    template <typename Vector>
    constexpr auto bindCoordReaders(Vector const& v);
    
    
    template <typename Dest, typename Source>
    Dest convertTo(Source const& v);
    
    
    template <typename Dest, typename Source>
    std::vector<Dest> convertCollTo(std::vector<Source> const& coll);
    
    
    template <typename Vector, typename Pred>
    Vector transformCoords(Vector const& v, Pred&& pred);
    
    // --- END Vector coordinate access abstraction ----------------------------
    
    
    // --- BEGIN Functions for common vector operations ------------------------

    template <typename Vector, typename Scalar>
    Vector rounded01(Vector const& v, Scalar tol)
      {
        return transformCoords
          (v, [tol](auto c){ return extra::roundValue01(c, tol); });
      }
    
    template <typename Vector, typename Scalar>
    void round01(Vector& v, Scalar tol) { v = rounded01(v, tol); }
    
    
    template <typename Vector>
    bool isfinite(Vector const& v);
    
    template <typename Vector>
    Vector normalize(Vector const& v) { return v.Unit(); }
    
    template <typename Vector>
    Vector cross(Vector const& a, Vector const& b) { return a.Cross(b); }
    
    template <typename Vector>
    constexpr auto dot(Vector const& a, Vector const& b) { return a.Dot(b); }
    
    template <typename Vector>
    auto mag2(Vector const& v) { return v.Mag2(); }
    
    template <typename Vector>
    auto norm(Vector const& v) { return v.Mag(); }
    
    template <typename Vector>
    auto mixedProduct(Vector const& a, Vector const& b, Vector const& c)
      { return dot(cross(a, b), c); }
    
    
    // --- END Functions for common vector operations --------------------------
    
    
    template <unsigned int N = 3U>
    class MiddlePointAccumulatorDim {
      static constexpr unsigned int Dim = N; 
      std::array<Length_t, Dim> fSums; 
      double fW = 0.0; 
      
        public:
      MiddlePointAccumulatorDim() { fSums.fill(0.); }
      
      template <typename BeginIter, typename EndIter>
      MiddlePointAccumulatorDim(BeginIter begin, EndIter end)
        : MiddlePointAccumulatorDim()
        { add(begin, end); }
      
      
      // --- BEGIN Result query ------------------------------------------------
      
      bool empty() const { return fW == 0.0; }
      
      double weight() const { return fW; }
      
      template <typename Point>
      Point middlePointAs() const
      //  { return { fSums[0] / fW, fSums[1] / fW, fSums[2] / fW }; }
        { return makeWeightedPoint<Point>(1.0 / fW); }
      
      geo::Point_t middlePoint() const
        { return middlePointAs<geo::Point_t>(); }
      
      
      // --- END Result query --------------------------------------------------
      
      // --- BEGIN Addition of points ------------------------------------------
      
      template <typename Point>
      void add(Point const& p)
        {
          std::size_t ic = 0U;
          for (auto c: geo::vect::bindCoordManagers(p)) fSums[ic++] += c();
          fW += 1.0;
        }
      
      template <typename Point>
      void add(Point const& p, double weight)
        {
          std::size_t ic = 0U;
          for (auto c: geo::vect::bindCoordManagers(p))
            fSums[ic++] += weight * c();
          fW += weight;
        }
      
      template <typename BeginIter, typename EndIter>
      void add(BeginIter begin, EndIter end)
        { std::for_each(begin, end, [this](auto const& p){ this->add(p); }); }
      
      void clear() { fSums.fill(0.); fW = 0.0; }
      
      // --- END Addition of points --------------------------------------------
      
        private:
      using IndexSequence_t = std::make_index_sequence<Dim>;
      
      template <typename Point, std::size_t... I>
      Point makePointImpl(std::index_sequence<I...>) const
        { return { fSums.operator[](I)... }; }
      
      template <typename Point, std::size_t... I>
      Point makeWeightedPointImpl(double w, std::index_sequence<I...>) const
        { return { (fSums.operator[](I) * w)... }; }
      
      template <typename Point>
      Point makePoint() const
        { return geo::vect::makeFromCoords<Point>(fSums); }
      
      template <typename Point>
      Point makeWeightedPoint(double w) const
        { return makeWeightedPointImpl<Point>(w, IndexSequence_t{}); }
      
      
    }; // MiddlePointAccumulatorDim()
    
    
    using MiddlePointAccumulator = MiddlePointAccumulatorDim<3U>;
    
    
    // --- BEGIN Middle point functions ----------------------------------------
    
    template <typename Point, typename BeginIter, typename EndIter>
    Point middlePointAs(BeginIter begin, EndIter end)
      {
        constexpr auto Dim = geo::vect::dimension<Point>();
        return MiddlePointAccumulatorDim<Dim>(begin, end)
          .template middlePointAs<Point>();
      }
    
    template <typename BeginIter, typename EndIter>
    geo::Point_t middlePoint(BeginIter begin, EndIter end)
      { return middlePointAs<geo::Point_t>(begin, end); }
    
    template <typename Point>
    Point middlePoint(std::initializer_list<Point> points)
      {
        constexpr auto Dim = geo::vect::dimension<Point>();
        return MiddlePointAccumulatorDim<Dim>(points.begin(), points.end())
          .template middlePointAs<Point>();
      }
    
    // --- END Middle point functions ------------------------------------------
    
    
    // --- BEGIN Support for LArSoft geometry vectors --------------------------
    
    // import global definitions
    using ::geo::Vector_t;
    using ::geo::Point_t;
    
    template <typename Point>
    ::geo::Point_t toPoint(Point const& p) 
      { return geo::vect::convertTo<::geo::Point_t>(p); }
    
    template <typename Vector>
    ::geo::Vector_t toVector(Vector const& v) 
      { return geo::vect::convertTo<::geo::Vector_t>(v); }
    
    
    template <typename Point>
    std::vector<geo::Point_t> convertCollToPoint
      (std::vector<Point> const& coll)
      { return convertCollTo<geo::Point_t>(coll); }
    
    template <typename Vector>
    std::vector<geo::Vector_t> convertCollToVector
      (std::vector<Vector> const& coll)
      { return convertCollTo<geo::Vector_t>(coll); }
    
    
    template <typename Coords>
    GENVECTOR_CONSTEXPR ::geo::Point_t makePointFromCoords(Coords&& coords)
      { return makeFromCoords<::geo::Point_t>(std::forward<Coords>(coords)); }
    
    template <typename Coords>
    GENVECTOR_CONSTEXPR ::geo::Vector_t makeVectorFromCoords(Coords&& coords)
      { return makeFromCoords<::geo::Vector_t>(std::forward<Coords>(coords)); }
    
    // --- END Support for LArSoft geometry vectors ----------------------------
    
  
  } // namespace vect
  
  // END Geometry group --------------------------------------------------------
  
} // namespace geo


//------------------------------------------------------------------------------
//---  template specializations for standard geometry vectors
//---
namespace geo {
  namespace vect {
    
    //--------------------------------------------------------------------------
    template <>
    inline auto norm(geo::Vector_t const& v) { return v.R(); }
    
    //--------------------------------------------------------------------------
    
  } // namespace vect
} // namespace geo


//------------------------------------------------------------------------------
//---  template implementation

namespace geo {
  namespace vect {
    namespace details {
      
      //------------------------------------------------------------------------
      template <typename> struct AlwaysFalse: std::false_type {};
      
      
      // constexpr variant of forward (from StackExchange)
      template <typename T>
      constexpr T&& constexpr_forward(std::remove_reference_t<T>& t) 
        { return static_cast<T&&>(t); }
      
      template<typename T>
      constexpr T&& constexpr_forward(std::remove_reference_t<T>&& t) 
        {
          static_assert(!std::is_lvalue_reference<T>(),
            "template argument substituting T is an lvalue reference type");
          return static_cast<T&&>(t);
        }
      
      //------------------------------------------------------------------------
      template <typename Vector>
      using VectorIndices_t = std::make_index_sequence<dimension<Vector>()>;
      
      template <typename Vector>
      constexpr auto makeVectorIndices() { return VectorIndices_t<Vector>{}; }
      
      template <typename Vector>
      constexpr auto makeVectorIndices(Vector&&)
        { return makeVectorIndices<Vector>(); }
      
      
      template <typename T, T... Indices>
      constexpr auto makeIndexSeqImpl(std::integer_sequence<T, Indices...>)
        // BUG the double brace syntax is required to work around clang bug 21629
        // (https://bugs.llvm.org/show_bug.cgi?id=21629)
        { return std::array<T, sizeof...(Indices)>{{ Indices... }}; }
      
      // fill a sequence object with the first `N` values of type `T`.
      template <typename T, T N>
      constexpr auto makeIndexSeq()
        { return makeIndexSeqImpl<T>(std::make_integer_sequence<T, N>{}); }
      
      
      //------------------------------------------------------------------------
      template <std::size_t I, typename Data>
      constexpr auto accessElement(Data&& data) { return data[I]; }
      
      template <typename Vector, typename Coords, std::size_t... Indices>
      constexpr Vector makeFromCoordsImpl
        (Coords&& coords, std::index_sequence<Indices...>)
        {
          return
            { accessElement<Indices>(constexpr_forward<Coords>(coords))... };
        }
      
      
      //------------------------------------------------------------------------
      template <typename Vector>
      constexpr CoordManager_t<Vector> NoCoordManager{ nullptr, nullptr };
      
      template <typename Vector, unsigned int Dim = dimension<Vector>()>
      struct CoordManagerImpl;
      
      template <typename Vector>
      struct CoordManagerImpl<Vector, 1U> {
        static auto get(unsigned int n) noexcept
          { return (n == 0)? XcoordManager<Vector>: NoCoordManager<Vector>; }
      }; // CoordManagerImpl<1U>
      
      template <typename Vector>
      struct CoordManagerImpl<Vector, 2U> {
        static auto get(unsigned int n) noexcept
          {
            return (n == 1)
              ? YcoordManager<Vector>: CoordManagerImpl<Vector, 1U>::get(n); 
          }
      }; // CoordManagerImpl<2U>
      
      template <typename Vector>
      struct CoordManagerImpl<Vector, 3U> {
        static auto get(unsigned int n) noexcept
          {
            return (n == 2)
              ? ZcoordManager<Vector>: CoordManagerImpl<Vector, 2U>::get(n);
          }
      }; // CoordManagerImpl<3U>
      
      template <typename Vector>
      struct CoordManagerImpl<Vector, 4U> {
        static auto get(unsigned int n) noexcept
          {
            return (n == 3)
              ? TcoordManager<Vector>: CoordManagerImpl<Vector, 3U>::get(n);
          }
          
      }; // CoordManagerImpl<4U>
      
      
      //------------------------------------------------------------------------
      template <typename Vector, unsigned int N>
      struct CoordManagersImplBase {
        static constexpr unsigned int Dim = N;
        
        using Manager_t = decltype(XcoordManager<Vector>);
        using Return_t = std::array<Manager_t, Dim>;
        
        static_assert(dimension<Vector>() == Dim, "Inconsistent vector size.");
      }; // CoordManagersImplBase
      
      template <typename Vector, unsigned int N>
      struct CoordManagersImpl;
      
      template <typename Vector>
      struct CoordManagersImpl<Vector, 2U>
        : private CoordManagersImplBase<Vector, 2U>
      {
        using Base_t = CoordManagersImplBase<Vector, 2U>;
        using typename Base_t::Return_t;
        static constexpr Return_t get()
          {
            // BUG the double brace syntax is required to work around clang bug 21629
            // (https://bugs.llvm.org/show_bug.cgi?id=21629)
            return {{
                XcoordManager<Vector>
              , YcoordManager<Vector>
            }};
          }
      }; // CoordManagersImpl<2U>
      
      template <typename Vector>
      struct CoordManagersImpl<Vector, 3U>
        : private CoordManagersImplBase<Vector, 3U>
      {
        using Base_t = CoordManagersImplBase<Vector, 3U>;
        using typename Base_t::Return_t;
        static constexpr Return_t get()
          {
            // BUG the double brace syntax is required to work around clang bug 21629
            // (https://bugs.llvm.org/show_bug.cgi?id=21629)
            return {{
                XcoordManager<Vector>
              , YcoordManager<Vector>
              , ZcoordManager<Vector>
            }};
          }
      }; // CoordManagersImpl<3U>
      
      template <typename Vector>
      struct CoordManagersImpl<Vector, 4U>
        : private CoordManagersImplBase<Vector, 4U>
      {
        using Base_t = CoordManagersImplBase<Vector, 4U>;
        using typename Base_t::Return_t;
        static constexpr Return_t get()
          {
            // BUG the double brace syntax is required to work around clang bug 21629
            // (https://bugs.llvm.org/show_bug.cgi?id=21629)
            return {{
                XcoordManager<Vector>
              , YcoordManager<Vector>
              , ZcoordManager<Vector>
              , TcoordManager<Vector>
            }};
          }
      }; // CoordManagersImpl<4U>
      
      
      //------------------------------------------------------------------------
      template <typename Vector, unsigned int N>
      struct BindCoordManagersImplBase {
        static constexpr unsigned int Dim = N;
        
        using Manager_t = decltype(Xcoord(std::declval<Vector>()));
        using Return_t = std::array<Manager_t, Dim>;
        
        static_assert(dimension<Vector>() == Dim, "Inconsistent vector size.");
      }; // CoordManagersImplBase
      
      template <typename Vector, unsigned int N>
      struct BindCoordManagersImpl;
      
      template <typename Vector>
      struct BindCoordManagersImpl<Vector, 2U>
        : private BindCoordManagersImplBase<Vector, 2U>
      {
        using Base_t = CoordManagersImplBase<Vector, 2U>;
        using typename Base_t::Return_t;
        static Return_t bind(Vector& v)
          {
            // BUG the double brace syntax is required to work around clang bug 21629
            // (https://bugs.llvm.org/show_bug.cgi?id=21629)
            return {{
                Xcoord(v)
              , Ycoord(v)
            }};
          }
      }; // BindCoordManagersImpl<2U>
      
      template <typename Vector>
      struct BindCoordManagersImpl<Vector, 3U>
        : private BindCoordManagersImplBase<Vector, 3U>
      {
        using Base_t = BindCoordManagersImplBase<Vector, 3U>;
        using typename Base_t::Return_t;
        static Return_t bind(Vector& v)
          {
            // BUG the double brace syntax is required to work around clang bug 21629
            // (https://bugs.llvm.org/show_bug.cgi?id=21629)
            return {{
                Xcoord(v)
              , Ycoord(v)
              , Zcoord(v)
            }};
          }
      }; // BindCoordManagersImpl<3U>
      
      template <typename Vector>
      struct BindCoordManagersImpl<Vector, 4U>
        : private BindCoordManagersImplBase<Vector, 4U>
      {
        using Base_t = BindCoordManagersImplBase<Vector, 4U>;
        using typename Base_t::Return_t;
        static Return_t bind(Vector& v)
          {
            // BUG the double brace syntax is required to work around clang bug 21629
            // (https://bugs.llvm.org/show_bug.cgi?id=21629)
            return {{
                Xcoord(v)
              , Ycoord(v)
              , Zcoord(v)
              , Tcoord(v)
            }};
          }
      }; // BindCoordManagersImpl<4U>
      
      
      //------------------------------------------------------------------------
      template<typename Dest, typename Source>
      struct ConvertToImplBase {
        static_assert(dimension<Source>() == dimension<Dest>(),
          "Source and destination vectors must have the same dimension.");
      }; // struct ConvertToImplBase
      
      
      // special pass-through case
      template <typename Dest, typename Source, unsigned int Dim>
      struct ConvertToImpl {
        // trivial to do: open a feature request!
        static_assert(
          AlwaysFalse<Dest>(),
          "This vector dimensionality is not implemented yet."
          );
      }; // struct ConvertToImpl
      
      template <typename Dest, typename Source>
      struct ConvertToImpl<Dest, Source, 2U>
        : private ConvertToImplBase<Dest, Source>
      {
        static Dest convert(Source const& v)
          { return { Xcoord(v)(), Ycoord(v)() }; }
      }; // struct ConvertToImpl<2U>
      
      template <typename Dest, typename Source>
      struct ConvertToImpl<Dest, Source, 3U>
        : private ConvertToImplBase<Dest, Source>
      {
        static Dest convert(Source const& v)
          { return { Xcoord(v)(), Ycoord(v)(), Zcoord(v)() }; }
      }; // struct ConvertToImpl<3U>
      
      template <typename Dest, typename Source>
      struct ConvertToImpl<Dest, Source, 4U>
        : private ConvertToImplBase<Dest, Source>
      {
        static Dest convert(Source const& v)
          { return { Xcoord(v)(), Ycoord(v)(), Zcoord(v)(), Tcoord(v)() }; }
      }; // struct ConvertToImpl<4U>
      
      
      template
        <typename Dest, typename Source, unsigned int Dim = dimension<Source>()>
      struct ConvertToDispatcher: public ConvertToImpl<Dest, Source, Dim> {};
      
      // special pass-through case
      template <typename Vector, unsigned int Dim>
      struct ConvertToDispatcher<Vector, Vector, Dim> {
        static_assert
          (Dim == dimension<Vector>(), "Inconsistent vector dimension");
        static constexpr Vector convert(Vector const& v) { return v; }
      }; // struct ConvertToDispatcher<pass through>
      
      
      //------------------------------------------------------------------------
      template <typename Point, std::size_t... I>
      bool isfiniteImpl(Point const& point, std::index_sequence<I...>)
        { return extended_and(std::isfinite(coord(point, I).get())...); }
      
      //------------------------------------------------------------------------
      
    } // namespace details
  } // namespace vect
} // namespace geo


//------------------------------------------------------------------------------
template <typename Vector>
constexpr std::array<std::size_t, geo::vect::dimension<Vector>()>
  geo::vect::indices()
{
  return details::makeIndexSeq<std::size_t, dimension<Vector>()>();
}

template <typename Vector>
constexpr auto geo::vect::indices(Vector const&) -> decltype(indices<Vector>())
  { return indices<Vector>(); }


//------------------------------------------------------------------------
template <typename Vector, typename Coords>
constexpr Vector geo::vect::makeFromCoords(Coords&& coords) {
  using namespace geo::vect::details;
  return makeFromCoordsImpl<Vector>
    (constexpr_forward<Coords>(coords), makeVectorIndices<Vector>());
} // geo::vect::makeFromCoords()


//------------------------------------------------------------------------------
template <typename Vector>
constexpr auto geo::vect::coordManager(unsigned int n)
  { return details::CoordManagerImpl<Vector>::get(n); }


//------------------------------------------------------------------------------
template <typename Vector>
auto geo::vect::coord(Vector& v, unsigned int n) noexcept {
  return vect::bindCoord<Vector>(v, coordManager<Vector>(n));
}


//------------------------------------------------------------------------------
template <typename Vector>
constexpr auto geo::vect::coordManagers() {
  using PlainVector = std::remove_reference_t<Vector>;
  return
    details::CoordManagersImpl<PlainVector, dimension<PlainVector>()>::get();
} // geo::vect::coordManagers()

template <typename Vector>
constexpr auto geo::vect::coordManagers(Vector&&)
  { return coordManagers<Vector>(); }


template <typename Vector>
constexpr auto geo::vect::coordReaders() {
  using ConstVector = std::add_const_t<std::remove_reference_t<Vector>>;
  return 
    details::CoordManagersImpl<ConstVector, dimension<ConstVector>()>::get();
} // geo::vect::coordReaders()

template <typename Vector>
constexpr auto geo::vect::coordReaders(Vector&&)
  { return coordReaders<Vector>(); }


//------------------------------------------------------------------------------
template <typename Vector>
constexpr auto geo::vect::bindCoordManagers(Vector& v) {
  return details::BindCoordManagersImpl<Vector, dimension<Vector>()>::bind(v);
} // geo::vect::bindCoordManagers()


template <typename Vector>
constexpr auto geo::vect::bindCoordReaders(Vector const& v) {
  using ConstVector = std::add_const_t<std::remove_reference_t<Vector>>;
  return details::BindCoordManagersImpl<ConstVector, dimension<ConstVector>()>
    ::bind(v);
} // geo::vect::bindCoordReaders()


//------------------------------------------------------------------------------
template <typename Dest, typename Source>
Dest geo::vect::convertTo(Source const& v)
  { return details::ConvertToDispatcher<Dest, Source>::convert(v); }


//----------------------------------------------------------------------------
template <typename Dest, typename Source>
std::vector<Dest> geo::vect::convertCollTo(std::vector<Source> const& coll) {
  
  std::vector<Dest> dest;
  dest.reserve(coll.size());
  std::transform(coll.begin(), coll.end(), std::back_inserter(dest),
    geo::vect::convertTo<Dest, Source>);
  return dest;
  
} // geo::vect::convertCollTo()


//------------------------------------------------------------------------
template <typename Vector, typename Pred>
Vector geo::vect::transformCoords(Vector const& v, Pred&& pred) {
  details::CoordinateArray_t<Vector> values;
  unsigned int i = 0;
  for (auto c: bindCoordReaders(v)) values[i++] = pred(c());
  return makeFromCoords<Vector>(values);
} // geo::vect::transformCoords()


//------------------------------------------------------------------------------
template <typename Vector>
bool geo::vect::isfinite(Vector const& v)
  { return details::isfiniteImpl(v, details::makeVectorIndices<Vector>()); }

//------------------------------------------------------------------------------


#endif // LARCOREOBJ_SIMPLETYPESANDCONSTANTS_GEO_VECTORS_UTILS_H