coordinates.hh

Go to the documentation of this file.

#pragma once
 
// standard library
#include <cmath>
#include <cstddef>
#include <cstdio>
#include <cstdlib>
#include <limits>
#include <stdexcept>
#include <string>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>
 
// DiFfRG
#include <DiFfRG/common/kokkos.hh>
#include <DiFfRG/common/utils.hh>
 
namespace DiFfRG
{
  template <typename T>
  concept is_coordinates = requires(T t) {
    typename T::ctype;
    T::dim;
    t.size();
    t.forward(device::array<size_t, T::dim>{});
    { t.from_linear_index(size_t{}) } -> std::same_as<device::array<size_t, T::dim>>;
  };
 
  namespace internal
  {
    template <typename T, typename = void> struct coord_periodic : std::false_type {
    };
    template <typename T>
    struct coord_periodic<T, std::void_t<decltype(T::periodic)>> : std::bool_constant<T::periodic> {
    };
  } // namespace internal
 
  template <typename T> inline constexpr bool is_periodic_coordinate_v = internal::coord_periodic<T>::value;
 
  namespace internal
  {
    template <typename C, size_t i, typename = void> struct axis_periodic : std::false_type {
    };
    template <typename C, size_t i>
    struct axis_periodic<C, i, std::void_t<typename C::template coordinate_type<i>>>
        : std::bool_constant<is_periodic_coordinate_v<typename C::template coordinate_type<i>>> {
    };
  } // namespace internal
 
  template <typename C, size_t i> inline constexpr bool is_periodic_axis_v = internal::axis_periodic<C, i>::value;
 
  namespace internal
  {
    template <typename NT> std::string round_trip(const NT value)
    {
      char buffer[40];
      for (int precision = 6; precision < std::numeric_limits<NT>::max_digits10; ++precision) {
        std::snprintf(buffer, sizeof(buffer), "%.*g", precision, double(value));
        if (NT(std::strtod(buffer, nullptr)) == value) return std::string(buffer);
      }
      std::snprintf(buffer, sizeof(buffer), "%.*g", std::numeric_limits<NT>::max_digits10, double(value));
      return std::string(buffer);
    }
  } // namespace internal
 
  template <typename... Coordinates> class CoordinatePackND
  {
    // Assert that all coordinates have the same ctype
    static_assert((std::is_same<typename Coordinates::ctype,
                                typename device::tuple_element<0, device::tuple<Coordinates...>>::type::ctype>::value &&
                   ...),
                  "CoordinatePackND: all coordinate systems must use the same ctype");
    static_assert(sizeof...(Coordinates) > 0, "CoordinatePackND requires at least one coordinate system");
    // Assert that all coordinates have the dim 1
    static_assert(((Coordinates::dim == 1) && ...), "CoordinatePackND requires all coordinates to have dim 1");
 
  public:
    using ctype = typename device::tuple_element<0, device::tuple<Coordinates...>>::type::ctype;
    static constexpr size_t dim = sizeof...(Coordinates);
 
    CoordinatePackND(Coordinates... coordinates) : coordinates(coordinates...) {}
 
    template <typename... Coordinates2>
      requires(sizeof...(Coordinates2) == sizeof...(Coordinates) &&
               !std::is_same_v<CoordinatePackND<Coordinates2...>, CoordinatePackND> &&
               (std::is_constructible_v<Coordinates, const Coordinates2 &> && ...))
    CoordinatePackND(const CoordinatePackND<Coordinates2...> &other)
        : CoordinatePackND(other, std::index_sequence_for<Coordinates...>{})
    {
    }
 
    template <typename... I>
      requires(std::is_convertible_v<std::decay_t<I>, size_t> && ...)
    KOKKOS_FORCEINLINE_FUNCTION device::array<ctype, dim> forward(I &&...i) const
    {
      static_assert(sizeof...(I) == sizeof...(Coordinates));
      return forward_impl(device::make_integer_sequence<int, sizeof...(I)>(), device::forward<I>(i)...);
    }
    template <typename IT>
    KOKKOS_FORCEINLINE_FUNCTION device::array<ctype, dim> forward(const device::array<IT, dim> &coords) const
    {
      return device::apply([&](const auto &...iargs) { return forward(iargs...); }, coords);
    }
 
    template <typename... I, int... Is>
    device::array<ctype, dim> KOKKOS_FORCEINLINE_FUNCTION forward_impl(device::integer_sequence<int, Is...>,
                                                                       I &&...i) const
    {
      return {{device::get<Is>(coordinates).forward(device::get<Is>(device::tie(i...)))...}};
    }
 
    template <typename... I> KOKKOS_FORCEINLINE_FUNCTION device::array<ctype, sizeof...(I)> backward(I &&...i) const
    {
      static_assert(sizeof...(I) == sizeof...(Coordinates));
      return backward_impl(device::make_integer_sequence<int, sizeof...(I)>(), device::forward<I>(i)...);
    }
 
    template <typename... I, int... Is>
    device::array<ctype, sizeof...(I)> KOKKOS_FORCEINLINE_FUNCTION backward_impl(device::integer_sequence<int, Is...>,
                                                                                 I &&...i) const
    {
      return {{device::get<Is>(coordinates).backward(device::get<Is>(device::tie(i...)))...}};
    }
 
    template <size_t i> const auto &get_coordinates() const { return device::get<i>(coordinates); }
 
    template <size_t i>
    using coordinate_type = typename device::tuple_element<i, device::tuple<Coordinates...>>::type;
 
    size_t KOKKOS_FORCEINLINE_FUNCTION size() const
    {
      size_t size = 1;
      constexpr_for<0, sizeof...(Coordinates), 1>([&](auto i) { size *= device::get<i>(coordinates).size(); });
      return size;
    }
 
    device::array<size_t, sizeof...(Coordinates)> KOKKOS_FORCEINLINE_FUNCTION sizes() const
    {
      device::array<size_t, sizeof...(Coordinates)> sizes;
      constexpr_for<0, sizeof...(Coordinates), 1>([&](auto i) { sizes[i] = device::get<i>(coordinates).size(); });
      return sizes;
    }
 
    device::array<size_t, sizeof...(Coordinates)> KOKKOS_INLINE_FUNCTION from_linear_index(size_t s) const
    {
      device::array<size_t, sizeof...(Coordinates)> idx;
      // calculate the index for each coordinate system
      constexpr_for<0, sizeof...(Coordinates), 1>([&](auto i) {
        idx[sizeof...(Coordinates) - 1 - i] = s % device::get<sizeof...(Coordinates) - 1 - i>(coordinates).size();
        s = s / device::get<sizeof...(Coordinates) - 1 - i>(coordinates).size();
      });
 
      return idx;
    }
 
    template <typename... Coordinates2>
    friend bool operator==(const CoordinatePackND<Coordinates...> &r, const CoordinatePackND<Coordinates2...> &l)
    {
      if constexpr (!std::is_same_v<device::tuple<Coordinates...>, device::tuple<Coordinates2...>>)
        return false; // different axes
      else {
        // Check if all coordinates are equal
        bool equal = true;
        constexpr_for<0, sizeof...(Coordinates), 1>(
            [&](auto i) { equal &= (device::get<i>(r.coordinates) == device::get<i>(l.coordinates)); });
        return equal;
      }
    }
 
    std::string to_string() const
    {
      std::string ret = "CoordinatePackND(";
      constexpr_for<0, dim, 1>([&](auto i) {
        ret += device::get<i>(coordinates).to_string();
        if (i + 1 < dim) ret += ", ";
      });
      ret += ")";
      return ret;
    }
 
  protected:
    const device::tuple<Coordinates...> coordinates;
 
  private:
    template <typename Other, size_t... Is>
    CoordinatePackND(const Other &other, std::index_sequence<Is...>)
        : coordinates(Coordinates(other.template get_coordinates<Is>())...)
    {
    }
  };
 
  template <typename Base> class SubCoordinates : public Base
  {
  public:
    using ctype = typename Base::ctype;
    static constexpr size_t dim = Base::dim;
 
    SubCoordinates(const Base &base, size_t offset, size_t size) : Base(base), m_offset(offset), m_size(size)
    {
      if (size == 0) throw std::runtime_error("SubCoordinates: size must be > 0");
      if (offset + size > base.size()) throw std::runtime_error("SubCoordinates: offset + size must be <= base.size()");
    }
 
    device::array<size_t, dim> KOKKOS_FORCEINLINE_FUNCTION sizes() const
    {
      static_assert(dim == 1, "SubCoordinates is a linear index window, which is a per-axis box only for dim == 1.");
      return device::array<size_t, dim>{{m_size}};
    }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION size() const { return m_size; }
 
    device::array<size_t, dim> KOKKOS_INLINE_FUNCTION from_linear_index(size_t s) const
    {
      return Base::from_linear_index(m_offset + s);
    }
 
    template <typename... I>
      requires(std::is_convertible_v<std::decay_t<I>, size_t> && ...)
    KOKKOS_FORCEINLINE_FUNCTION device::array<ctype, dim> forward(I &&...i) const
    {
      static_assert(sizeof...(I) == dim);
      return forward(device::array<size_t, dim>{{static_cast<size_t>(i)...}});
    }
 
    template <typename IT> KOKKOS_INLINE_FUNCTION device::array<ctype, dim> forward(device::array<IT, dim> i) const
    {
      return Base::forward(i);
    }
 
    // backward() is Base's.
 
    std::string to_string() const
    {
      return "SubCoordinates(" + Base::to_string() + ", " + std::to_string(m_offset) + ", " + std::to_string(m_size) +
             ")";
    }
 
  private:
    size_t m_offset;
    size_t m_size;
  };
 
  template <typename NT = double>
    requires std::is_floating_point_v<NT>
  class LinearCoordinates1D
  {
  public:
    using ctype = NT;
    static constexpr size_t dim = 1;
 
    LinearCoordinates1D(size_t grid_extent, double start, double stop)
        : start(start), stop(stop), grid_extent(grid_extent), parent_extent(grid_extent)
    {
      if (grid_extent < 2) throw std::runtime_error("LinearCoordinates1D: grid_extent must be > 1");
      a = (stop - start) / (grid_extent - 1.);
    }
 
    template <typename NT2>
    LinearCoordinates1D(const LinearCoordinates1D<NT2> &other)
        : LinearCoordinates1D(LinearCoordinates1D(other.parent_size(), other.start, other.stop), other.size())
    {
    }
 
    LinearCoordinates1D prefix(size_t n) const
    {
      if (n < 1 || n > grid_extent) throw std::runtime_error("LinearCoordinates1D::prefix: need 1 <= n <= size()");
      return LinearCoordinates1D(*this, n);
    }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION parent_size() const { return parent_extent; }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION from_linear_index(size_t i) const
    {
      return device::array<size_t, 1>{i};
    }
 
    template <typename IT> NT KOKKOS_FORCEINLINE_FUNCTION forward(const IT &x) const { return start + a * x; }
 
    template <typename IT> device::array<NT, 1> KOKKOS_FORCEINLINE_FUNCTION forward(const device::array<IT, 1> &x) const
    {
      return {forward(x[0])};
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION backward(const NT &y) const { return (y - start) / a; }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION size() const { return grid_extent; }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION sizes() const { return {grid_extent}; }
 
    const NT start, stop;
 
    template <typename NT2>
    friend bool operator==(const LinearCoordinates1D<NT> &lhs, const LinearCoordinates1D<NT2> &rhs)
    {
      if constexpr (!std::is_same_v<NT, NT2>) return false; // Different types, cannot be equal
      return lhs.start == rhs.start && lhs.stop == rhs.stop && lhs.size() == rhs.size() &&
             lhs.parent_size() == rhs.parent_size();
    }
 
    std::string to_string() const
    {
      return "LinearCoordinates1D(" + std::to_string(parent_extent) + ", " + internal::round_trip(start) + ", " +
             internal::round_trip(stop) + ")" +
             (grid_extent == parent_extent ? std::string() : ".prefix(" + std::to_string(grid_extent) + ")");
    }
 
  private:
    // prefix(): the parent's parameters and spacing, fewer points
    LinearCoordinates1D(const LinearCoordinates1D &parent, size_t n)
        : start(parent.start), stop(parent.stop), grid_extent(n), parent_extent(parent.parent_extent), a(parent.a)
    {
    }
 
    const size_t grid_extent;
    size_t parent_extent;
    NT a;
  };
 
  template <typename NT = double>
    requires std::is_floating_point_v<NT>
  class LinearPeriodicCoordinates1D
  {
  public:
    using ctype = NT;
    static constexpr size_t dim = 1;
    static constexpr bool periodic = true;
 
    LinearPeriodicCoordinates1D(size_t grid_extent, double start, double stop)
        : start(start), stop(stop), grid_extent(grid_extent), extent(static_cast<NT>(grid_extent))
    {
      if (grid_extent == 0) throw std::runtime_error("LinearPeriodicCoordinates1D: grid_extent must be > 0");
      a = (stop - start) / NT(grid_extent);
    }
 
    template <typename NT2>
    LinearPeriodicCoordinates1D(const LinearPeriodicCoordinates1D<NT2> &other)
        : LinearPeriodicCoordinates1D(other.size(), other.start, other.stop)
    {
    }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION from_linear_index(size_t i) const
    {
      return device::array<size_t, 1>{i};
    }
 
    template <typename IT> NT KOKKOS_FORCEINLINE_FUNCTION forward(const IT &x) const { return start + a * x; }
 
    template <typename IT> device::array<NT, 1> KOKKOS_FORCEINLINE_FUNCTION forward(const device::array<IT, 1> &x) const
    {
      return {forward(x[0])};
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION backward(const NT &y) const
    {
      using Kokkos::floor;
      // fmod would keep the sign of idx, so fold explicitly into [0, grid_extent)
      NT idx = (y - start) / a;
      idx -= floor(idx / extent) * extent;
      // right at the seam the division can round such that the fold leaves idx just outside the interval
      if (idx >= extent) idx -= extent;
      if (idx < NT(0)) idx = NT(0);
      return idx;
    }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION size() const { return grid_extent; }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION sizes() const { return {grid_extent}; }
 
    const NT start, stop;
 
    template <typename NT2>
    friend bool operator==(const LinearPeriodicCoordinates1D<NT> &lhs, const LinearPeriodicCoordinates1D<NT2> &rhs)
    {
      if constexpr (!std::is_same_v<NT, NT2>) return false; // Different types, cannot be equal
      return lhs.start == rhs.start && lhs.stop == rhs.stop && lhs.grid_extent == rhs.grid_extent;
    }
 
    std::string to_string() const
    {
      return "LinearPeriodicCoordinates1D(" + std::to_string(grid_extent) + ", " + internal::round_trip(start) + ", " +
             internal::round_trip(stop) + ")";
    }
 
  private:
    const size_t grid_extent;
    const NT extent;
    NT a;
  };
 
  template <typename NT = double>
    requires std::is_floating_point_v<NT>
  class LogarithmicCoordinates1D
  {
  public:
    using ctype = NT;
    static constexpr size_t dim = 1;
 
    LogarithmicCoordinates1D(size_t grid_extent, double start, double stop, double bias)
        : start(start), stop(stop), bias(bias), grid_extent(grid_extent), parent_extent(grid_extent),
          gem1(grid_extent - 1.), gem1inv(1. / (grid_extent - 1.)), a(bias), b((stop - start) / std::expm1(bias)),
          c(start)
    {
      if (grid_extent < 2) throw std::runtime_error("LogarithmicCoordinates1D: grid_extent must be > 1");
      if (bias == 0.) throw std::runtime_error("LogarithmicCoordinates1D: bias must be != 0");
    }
 
    template <typename NT2>
    LogarithmicCoordinates1D(const LogarithmicCoordinates1D<NT2> &other)
        : LogarithmicCoordinates1D(LogarithmicCoordinates1D(other.parent_size(), other.start, other.stop, other.bias),
                                   other.size())
    {
    }
 
    LogarithmicCoordinates1D prefix(size_t n) const
    {
      if (n < 1 || n > grid_extent) throw std::runtime_error("LogarithmicCoordinates1D::prefix: need 1 <= n <= size()");
      return LogarithmicCoordinates1D(*this, n);
    }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION parent_size() const { return parent_extent; }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION from_linear_index(size_t i) const
    {
      return device::array<size_t, 1>{i};
    }
 
    template <typename IT> NT KOKKOS_FORCEINLINE_FUNCTION forward(const IT &x) const
    {
      using Kokkos::expm1;
      return b * expm1(a * x * gem1inv) + c;
    }
 
    template <typename IT> device::array<NT, 1> KOKKOS_FORCEINLINE_FUNCTION forward(const device::array<IT, 1> &x) const
    {
      return {forward(x[0])};
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION backward(const NT &y) const
    {
      using Kokkos::log1p;
      return log1p((y - c) / b) * gem1 / a;
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION backward_derivative(const NT &y) const { return gem1 / (a * (b + y - c)); }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION size() const { return grid_extent; }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION sizes() const { return {grid_extent}; }
 
    const NT start, stop, bias;
 
    template <typename NT2>
    friend bool operator==(const LogarithmicCoordinates1D<NT> &lhs, const LogarithmicCoordinates1D<NT2> &rhs)
    {
      if constexpr (!std::is_same_v<NT, NT2>) return false; // Different types, cannot be equal
      return lhs.start == rhs.start && lhs.stop == rhs.stop && lhs.bias == rhs.bias && lhs.size() == rhs.size() &&
             lhs.parent_size() == rhs.parent_size();
    }
 
    std::string to_string() const
    {
      return "LogarithmicCoordinates1D(" + internal::round_trip(start) + ", " + internal::round_trip(stop) + ", " +
             internal::round_trip(bias) + ", " + std::to_string(parent_extent) + ")" +
             (grid_extent == parent_extent ? std::string() : ".prefix(" + std::to_string(grid_extent) + ")");
    }
 
  private:
    // prefix(): the parent's parameters and spacing (gem1 is the parent's N - 1), fewer points
    LogarithmicCoordinates1D(const LogarithmicCoordinates1D &parent, size_t n)
        : start(parent.start), stop(parent.stop), bias(parent.bias), grid_extent(n),
          parent_extent(parent.parent_extent), gem1(parent.gem1), gem1inv(parent.gem1inv), a(parent.a), b(parent.b),
          c(parent.c)
    {
    }
 
    const size_t grid_extent;
    size_t parent_extent;
    const NT gem1, gem1inv;
    NT a, b, c;
  };
 
  template <typename NT = double>
    requires std::is_floating_point_v<NT>
  class FocusedLogCoordinates1D
  {
  public:
    using ctype = NT;
    static constexpr size_t dim = 1;
 
    FocusedLogCoordinates1D(size_t grid_extent, double start, double stop, double center, double focus)
        : start(start), stop(stop), center(center), focus(focus), grid_extent(grid_extent), parent_extent(grid_extent)
    {
      if (grid_extent < 2) throw std::runtime_error("FocusedLogCoordinates1D: grid_extent must be > 1");
      if (!(start > 0.)) throw std::runtime_error("FocusedLogCoordinates1D: start must be > 0");
      if (!(stop > start)) throw std::runtime_error("FocusedLogCoordinates1D: stop must be > start");
      if (!(center > 0.)) throw std::runtime_error("FocusedLogCoordinates1D: center must be > 0");
      if (!(focus >= 0.)) throw std::runtime_error("FocusedLogCoordinates1D: focus must be >= 0");
 
      const double u_start = std::log(start);
      const double u_stop = std::log(stop);
      const double u0_d = std::log(center);
      // sinh(s)/c is 0/0 at c == 0; below this threshold the map is a pure logarithmic grid to
      // within round-off anyway, so switch to it explicitly.
      pure_log = !(focus > 1e-8);
 
      double s_min_d, a_d;
      if (pure_log) {
        // in this branch s *is* u, so forward() degenerates to start * exp(u - log(start))
        s_min_d = u_start;
        a_d = (u_stop - u_start) / double(grid_extent - 1);
      } else {
        s_min_d = s_of_u(u_start, focus, u0_d);
        a_d = (s_of_u(u_stop, focus, u0_d) - s_min_d) / double(grid_extent - 1);
      }
      u0 = u0_d;
      c = focus;
      c_inv = pure_log ? NT(1) : NT(1. / focus);
      s_min = s_min_d;
      a = a_d;
      a_inv = NT(1. / a_d);
      // Store sinh(s_min), NOT sinh(s_min) * c_inv: forward() subtracts this BEFORE scaling by
      // c_inv, which is what makes x == 0 cancel exactly. See the note in forward(). For the same
      // reason it is evaluated in NT on the ROUNDED s_min, i.e. on exactly what forward() sees.
      using Kokkos::sinh;
      g_min = pure_log ? s_min : sinh(s_min);
    }
 
    template <typename NT2>
    FocusedLogCoordinates1D(const FocusedLogCoordinates1D<NT2> &other)
        : FocusedLogCoordinates1D(
              FocusedLogCoordinates1D(other.parent_size(), other.start, other.stop, other.center, other.focus),
              other.size())
    {
    }
 
    FocusedLogCoordinates1D prefix(size_t n) const
    {
      if (n < 1 || n > grid_extent) throw std::runtime_error("FocusedLogCoordinates1D::prefix: need 1 <= n <= size()");
      return FocusedLogCoordinates1D(*this, n);
    }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION parent_size() const { return parent_extent; }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION from_linear_index(size_t i) const
    {
      return device::array<size_t, 1>{i};
    }
 
    template <typename IT> NT KOKKOS_FORCEINLINE_FUNCTION forward(const IT &x) const
    {
      using Kokkos::exp, Kokkos::sinh;
      const NT s = s_min + a * x;
      // start * exp(...) and subtracting g_min before scaling by c_inv (no FMA contraction) make x == 0 give start
      // exactly.
      const NT d = pure_log ? (s - g_min) : (sinh(s) - g_min) * c_inv;
      return start * exp(d);
    }
 
    template <typename IT> device::array<NT, 1> KOKKOS_FORCEINLINE_FUNCTION forward(const device::array<IT, 1> &x) const
    {
      return {forward(x[0])};
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION backward(const NT &y) const
    {
      using Kokkos::log;
      const NT u = log(y);
      return ((pure_log ? u : s_of_u(u, c, u0)) - s_min) * a_inv;
    }
 
    NT KOKKOS_FORCEINLINE_FUNCTION backward_derivative(const NT &y) const
    {
      using Kokkos::log, Kokkos::sqrt;
      if (pure_log) return a_inv / y;
      const NT v = c * (log(y) - u0);
      return a_inv * c / (sqrt(v * v + NT(1)) * y);
    }
 
    size_t KOKKOS_FORCEINLINE_FUNCTION size() const { return grid_extent; }
 
    device::array<size_t, 1> KOKKOS_FORCEINLINE_FUNCTION sizes() const { return {grid_extent}; }
 
    const NT start, stop, center, focus;
 
    template <typename NT2>
    friend bool operator==(const FocusedLogCoordinates1D<NT> &lhs, const FocusedLogCoordinates1D<NT2> &rhs)
    {
      if constexpr (!std::is_same_v<NT, NT2>) return false; // Different types, cannot be equal
      return lhs.start == rhs.start && lhs.stop == rhs.stop && lhs.center == rhs.center && lhs.focus == rhs.focus &&
             lhs.size() == rhs.size() && lhs.parent_size() == rhs.parent_size();
    }
 
    std::string to_string() const
    {
      return "FocusedLogCoordinates1D(" + std::to_string(parent_extent) + ", " + internal::round_trip(start) + ", " +
             internal::round_trip(stop) + ", " + internal::round_trip(center) + ", " + internal::round_trip(focus) +
             ")" + (grid_extent == parent_extent ? std::string() : ".prefix(" + std::to_string(grid_extent) + ")");
    }
 
  private:
    template <typename T> static KOKKOS_FORCEINLINE_FUNCTION T s_of_u(const T u, const T c, const T u0)
    {
      using Kokkos::log, Kokkos::sqrt;
      const T v = c * (u - u0);
      const T D = sqrt(v * v + T(1));
      return v >= T(0) ? log(v + D) : -log(D - v);
    }
 
    // prefix(): the parent's parameters and spacing, fewer points
    FocusedLogCoordinates1D(const FocusedLogCoordinates1D &parent, size_t n)
        : start(parent.start), stop(parent.stop), center(parent.center), focus(parent.focus), grid_extent(n),
          parent_extent(parent.parent_extent), u0(parent.u0), s_min(parent.s_min), a(parent.a), a_inv(parent.a_inv),
          c(parent.c), c_inv(parent.c_inv), g_min(parent.g_min), pure_log(parent.pure_log)
    {
    }
 
    const size_t grid_extent;
    size_t parent_extent;
    NT u0, s_min, a, a_inv, c, c_inv, g_min;
    bool pure_log;
  };
 
  template <typename Coordinates> auto make_grid(const Coordinates &coordinates)
  {
    using ctype = typename Coordinates::ctype;
    using cortype = device::array<ctype, Coordinates::dim>;
    std::vector<cortype> grid(coordinates.size());
    for (size_t i = 0; i < coordinates.size(); ++i) {
      const auto forwarded = coordinates.forward(coordinates.from_linear_index(i));
      for (size_t j = 0; j < Coordinates::dim; ++j) {
        grid[i][j] = forwarded[j];
      }
    }
    return grid;
  }
 
  template <typename Coordinates> std::vector<typename Coordinates::ctype> dump_grid(const Coordinates &coordinates)
  {
    using ctype = typename Coordinates::ctype;
    std::vector<typename Coordinates::ctype> grid(coordinates.size() * Coordinates::dim);
    for (size_t i = 0; i < coordinates.size(); ++i) {
      for (size_t j = 0; j < Coordinates::dim; ++j) {
        const auto forwarded = coordinates.forward(coordinates.from_linear_index(i));
        grid[i * coordinates.dim + j] = forwarded[j];
      }
    }
    return grid;
  }
 
  template <typename Coordinates, typename NT> struct rebind_ctype;
  template <typename NT0, typename NT> struct rebind_ctype<LinearCoordinates1D<NT0>, NT> {
    using type = LinearCoordinates1D<NT>;
  };
  template <typename NT0, typename NT> struct rebind_ctype<LinearPeriodicCoordinates1D<NT0>, NT> {
    using type = LinearPeriodicCoordinates1D<NT>;
  };
  template <typename NT0, typename NT> struct rebind_ctype<LogarithmicCoordinates1D<NT0>, NT> {
    using type = LogarithmicCoordinates1D<NT>;
  };
  template <typename NT0, typename NT> struct rebind_ctype<FocusedLogCoordinates1D<NT0>, NT> {
    using type = FocusedLogCoordinates1D<NT>;
  };
  template <typename... Cs, typename NT>
    requires(requires { typename rebind_ctype<Cs, NT>::type; } && ...)
  struct rebind_ctype<CoordinatePackND<Cs...>, NT> {
    using type = CoordinatePackND<typename rebind_ctype<Cs, NT>::type...>;
  };
  template <typename Coordinates, typename NT> using rebind_ctype_t = typename rebind_ctype<Coordinates, NT>::type;
 
  // Common coordinate systems (double precision)
  using LogCoordinates = LogarithmicCoordinates1D<double>;
  using LinCoordinates = LinearCoordinates1D<double>;
  using LogLogCoordinates = CoordinatePackND<LogarithmicCoordinates1D<double>, LogarithmicCoordinates1D<double>>;
  using LogLinCoordinates = CoordinatePackND<LogarithmicCoordinates1D<double>, LinearCoordinates1D<double>>;
  using LinLogCoordinates = CoordinatePackND<LinearCoordinates1D<double>, LogarithmicCoordinates1D<double>>;
  using LinLinCoordinates = CoordinatePackND<LinearCoordinates1D<double>, LinearCoordinates1D<double>>;
  using LogLogLinCoordinates =
      CoordinatePackND<LogarithmicCoordinates1D<double>, LogarithmicCoordinates1D<double>, LinearCoordinates1D<double>>;
  using LogLinLinCoordinates =
      CoordinatePackND<LogarithmicCoordinates1D<double>, LinearCoordinates1D<double>, LinearCoordinates1D<double>>;
  using LinLinLinCoordinates =
      CoordinatePackND<LinearCoordinates1D<double>, LinearCoordinates1D<double>, LinearCoordinates1D<double>>;
 
  // Logarithmic coordinates focused on an interior scale, see FocusedLogCoordinates1D
  using FocusedLogCoordinates = FocusedLogCoordinates1D<double>;
  using FocusedLogLinCoordinates = CoordinatePackND<FocusedLogCoordinates1D<double>, LinearCoordinates1D<double>>;
  using FocusedLogLinLinCoordinates = CoordinatePackND<FocusedLogCoordinates1D<double>, LinearCoordinates1D<double>,
                                                       LinearCoordinates1D<double>>;
 
  // Periodic (angular) coordinates
  using LinPeriodicCoordinates = LinearPeriodicCoordinates1D<double>;
  using LogLinPeriodicCoordinates =
      CoordinatePackND<LogarithmicCoordinates1D<double>, LinearPeriodicCoordinates1D<double>>;
  using LogLinLinPeriodicCoordinates = CoordinatePackND<LogarithmicCoordinates1D<double>, LinearCoordinates1D<double>,
                                                        LinearPeriodicCoordinates1D<double>>;
  using FocusedLogLinPeriodicCoordinates =
      CoordinatePackND<FocusedLogCoordinates1D<double>, LinearPeriodicCoordinates1D<double>>;
  using FocusedLogLinLinPeriodicCoordinates =
      CoordinatePackND<FocusedLogCoordinates1D<double>, LinearCoordinates1D<double>,
                       LinearPeriodicCoordinates1D<double>>;
} // namespace DiFfRG