fastor_vector.hpp

// This file is part of the ACTS project.
//
// Copyright (C) 2016 CERN for the benefit of the ACTS project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
#pragma once
 
// Project include(s).
#include "detray/algebra/common/math.hpp"
#include "detray/algebra/concepts.hpp"
#include "detray/definitions/detail/qualifiers.hpp"
 
// Fastor include(s).
#ifdef _MSC_VER
#pragma warning(disable : 4244 4701 4702)
#endif  // MSVC
#include <Fastor/Fastor.h>
#ifdef _MSC_VER
#pragma warning(default : 4244 4701 4702)
#endif  // MSVC
 
namespace detray::algebra::fastor::math {
 
// Note that for all `Fastor::AbstractTensors`, the `N` template parameters
// refers to the number of dimensions, not the number of elements.

template <typename Derived, auto N>
DETRAY_HOST_DEVICE constexpr auto phi(
    const Fastor::AbstractTensor<Derived, N> &a) {
  // we first force evaluation of whatever was passed in.
  auto v = Fastor::evaluate(a);
  // using `auto` relieves us from having to extract the exact dimension of
  // the vector from somewhere. For all intents and purposes, we can consider
  // the type to be `Fastor::Tensor<scalar_t, SIZE>`.
 
  // we use the cmath version of `atan2` because Fastor's `atan2` works on
  // `Fastor::Tensor`s element-wise, which we don't want.
  return algebra::math::atan2(v[1], v[0]);
}

template <typename Derived, auto N>
DETRAY_HOST constexpr auto theta(
    const Fastor::AbstractTensor<Derived, N> &a) noexcept {
  // we first force evaluation of whatever was passed in.
  auto v = Fastor::evaluate(a);
  // using `auto` relieves us from having to extract the exact dimension of
  // the vector from somewhere. For all intents and purposes, we can consider
  // the type to be `Fastor::Tensor<scalar_t, SIZE>`.
 
  // we use the cmath version of `atan2` because Fastor's `atan2` works on
  // `Fastor::Tensor`s element-wise, which we don't want.
  return algebra::math::atan2(Fastor::norm(v(Fastor::fseq<0, 2>())), v[2]);
}

template <typename Derived, auto N>
DETRAY_HOST constexpr auto perp(
    const Fastor::AbstractTensor<Derived, N> &a) noexcept {
  // we first force evaluation of whatever was passed in.
  // using `auto` relieves us from having to extract the exact dimension of
  // the vector from somewhere. For all intents and purposes, we can consider
  // the type to be `Fastor::Tensor<scalar_t, SIZE>`.
  auto v = Fastor::evaluate(a);
 
  // we use the cmath version of `sqrt` because Fastor's `sqrt` works on
  // `Fastor::Tensor`s element-wise, which we don't want.
  return algebra::math::sqrt(
      Fastor::inner(v(Fastor::fseq<0, 2>()), v(Fastor::fseq<0, 2>())));
}

template <typename Derived, auto N>
DETRAY_HOST constexpr auto norm(const Fastor::AbstractTensor<Derived, N> &v) {
  return Fastor::norm(v);
}

template <typename Derived, auto N>
DETRAY_HOST constexpr auto eta(
    const Fastor::AbstractTensor<Derived, N> &a) noexcept {
  auto v = Fastor::evaluate(a);
 
  // we use the cmath version of `atanh` because Fastor's `atanh` works on
  // `Fastor::Tensor`s element-wise, which we don't want.
  return algebra::math::atanh(v[2] / Fastor::norm(v));
}

template <typename Derived, auto N>
DETRAY_HOST constexpr auto normalize(
    const Fastor::AbstractTensor<Derived, N> &v) {
  return v / Fastor::norm(v);
}

template <typename Derived0, auto N0, typename Derived1, auto N1>
DETRAY_HOST constexpr auto dot(const Fastor::AbstractTensor<Derived0, N0> &a,
                               const Fastor::AbstractTensor<Derived1, N1> &b) {
  return Fastor::inner(a, b);
}

template <concepts::scalar scalar_t, auto N>
DETRAY_HOST_DEVICE constexpr scalar_t dot(
    const Fastor::Tensor<scalar_t, N> &a,
    const Fastor::Tensor<scalar_t, N> &b) {
  return Fastor::inner(a, b);
}

template <concepts::scalar scalar_t, auto N>
DETRAY_HOST constexpr scalar_t dot(const Fastor::Tensor<scalar_t, N> &a,
                                   const Fastor::Tensor<scalar_t, N, 1> &b) {
  // We need to specify the type of the Tensor slice because Fastor by default
  // is lazy, so it returns an intermediate type which does not play well with
  // the Fastor::inner function.
  return Fastor::inner(a,
                       Fastor::Tensor<scalar_t, N>(b(Fastor::fseq<0, N>(), 0)));
}

template <concepts::scalar scalar_t, auto N>
DETRAY_HOST constexpr scalar_t dot(const Fastor::Tensor<scalar_t, N, 1> &a,
                                   const Fastor::Tensor<scalar_t, N> &b) {
  return Fastor::inner(Fastor::Tensor<scalar_t, N>(a(Fastor::fseq<0, N>(), 0)),
                       b);
}

template <concepts::scalar scalar_t, auto N>
DETRAY_HOST constexpr scalar_t dot(const Fastor::Tensor<scalar_t, N, 1> &a,
                                   const Fastor::Tensor<scalar_t, N, 1> &b) {
  return Fastor::inner(Fastor::Tensor<scalar_t, N>(a(Fastor::fseq<0, 3>(), 0)),
                       Fastor::Tensor<scalar_t, N>(b(Fastor::fseq<0, 3>(), 0)));
}

template <typename Derived0, auto N0, typename Derived1, auto N1>
DETRAY_HOST constexpr auto cross(
    const Fastor::AbstractTensor<Derived0, N0> &a,
    const Fastor::AbstractTensor<Derived1, N1> &b) {
  return Fastor::cross(a, b);
}

template <concepts::scalar scalar_t>
DETRAY_HOST constexpr auto cross(const Fastor::Tensor<scalar_t, 3> &a,
                                 const Fastor::Tensor<scalar_t, 3> &b) {
  return Fastor::cross(a, b);
}

template <concepts::scalar scalar_t>
DETRAY_HOST constexpr auto cross(const Fastor::Tensor<scalar_t, 3> &a,
                                 const Fastor::Tensor<scalar_t, 3, 1> &b) {
  // We need to specify the type of the Tensor slice because Fastor by default
  // is lazy, so it returns an intermediate type which does not play well with
  // the Fastor::cross function.
  return Fastor::cross(a,
                       Fastor::Tensor<scalar_t, 3>(b(Fastor::fseq<0, 3>(), 0)));
}

template <concepts::scalar scalar_t>
DETRAY_HOST constexpr auto cross(const Fastor::Tensor<scalar_t, 3, 1> &a,
                                 const Fastor::Tensor<scalar_t, 3> &b) {
  return Fastor::cross(Fastor::Tensor<scalar_t, 3>(a(Fastor::fseq<0, 3>(), 0)),
                       b);
}

template <concepts::scalar scalar_t>
DETRAY_HOST constexpr auto cross(const Fastor::Tensor<scalar_t, 3, 1> &a,
                                 const Fastor::Tensor<scalar_t, 3, 1> &b) {
  return Fastor::cross(Fastor::Tensor<scalar_t, 3>(a(Fastor::fseq<0, 3>(), 0)),
                       Fastor::Tensor<scalar_t, 3>(b(Fastor::fseq<0, 3>(), 0)));
}
 
}  // namespace detray::algebra::fastor::math