contrib/Eigen/src/Tensor/TensorFunctors.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2014 Benoit Steiner <benoit.steiner.goog@gmail.com>
//
// 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 http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0

#ifndef EIGEN_TENSOR_TENSOR_FUNCTORS_H
#define EIGEN_TENSOR_TENSOR_FUNCTORS_H

// IWYU pragma: private
#include "./InternalHeaderCheck.h"

namespace Eigen {
namespace internal {

/** \internal
 * \brief Template functor to compute the modulo between an array and a scalar.
 */
template <typename Scalar>
struct scalar_mod_op {
  EIGEN_DEVICE_FUNC scalar_mod_op(const Scalar& divisor) : m_divisor(divisor) {}
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar operator()(const Scalar& a) const { return a % m_divisor; }
  const Scalar m_divisor;
};
template <typename Scalar>
struct functor_traits<scalar_mod_op<Scalar>> {
  enum { Cost = scalar_div_cost<Scalar, false>::value, PacketAccess = false };
};

/** \internal
 * \brief Template functor to compute the modulo between 2 arrays.
 */
template <typename Scalar>
struct scalar_mod2_op {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar operator()(const Scalar& a, const Scalar& b) const { return a % b; }
};
template <typename Scalar>
struct functor_traits<scalar_mod2_op<Scalar>> {
  enum { Cost = scalar_div_cost<Scalar, false>::value, PacketAccess = false };
};

template <typename Scalar>
struct scalar_fmod_op {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar operator()(const Scalar& a, const Scalar& b) const {
    return numext::fmod(a, b);
  }
};
template <typename Scalar>
struct functor_traits<scalar_fmod_op<Scalar>> {
  enum {
    Cost = 13,  // Reciprocal throughput of FPREM on Haswell.
    PacketAccess = false
  };
};

template <typename Reducer, typename Device>
struct reducer_traits {
  enum { Cost = 1, PacketAccess = false, IsStateful = false, IsExactlyAssociative = true };
};

// Marks reducers whose partial accumulators may be accumulated independently and combined in any
// order. That requires reduce(value, accum) to be a pure combine -- so a partial accumulator may
// be fed back through it -- and that combine to be associative and commutative with initialize()
// as its identity, since the accumulators receive interleaved, not contiguous, operands.
//
// Neither property is implied by PacketAccess: a reducer may support packets yet transform each
// value it accepts (e.g. a sum-of-squares reducer), in which case merging through reduce() is
// wrong; and a pure combine may still be order-dependent (e.g. composition). Conservative
// default; built-in combines opt in below, gated on scalar categories whose semantics Eigen
// controls.
template <typename Reducer>
struct reducer_can_reorder_accumulators : std::false_type {};

// Standard reduction functors
template <typename T>
struct SumReducer {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const T t, T* accum) const {
    internal::scalar_sum_op<T> sum_op;
    *accum = sum_op(*accum, t);
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reducePacket(const Packet& p, Packet* accum) const {
    (*accum) = padd<Packet>(*accum, p);
  }

  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T initialize() const {
    internal::scalar_cast_op<int, T> conv;
    return conv(0);
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet initializePacket() const {
    return pset1<Packet>(initialize());
  }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T finalize(const T accum) const { return accum; }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet finalizePacket(const Packet& vaccum) const {
    return vaccum;
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T finalizeBoth(const T saccum, const Packet& vaccum) const {
    internal::scalar_sum_op<T> sum_op;
    return sum_op(saccum, predux(vaccum));
  }
};

template <typename T, typename Device>
struct reducer_traits<SumReducer<T>, Device> {
  enum {
    Cost = NumTraits<T>::AddCost,
    PacketAccess = PacketType<T, Device>::HasAdd,
    IsStateful = false,
    IsExactlyAssociative = NumTraits<T>::IsInteger
  };
};

template <typename T>
struct MeanReducer {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE MeanReducer() = default;

  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const T t, T* accum) {
    internal::scalar_sum_op<T> sum_op;
    *accum = sum_op(*accum, t);
    scalarCount_++;
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reducePacket(const Packet& p, Packet* accum) {
    (*accum) = padd<Packet>(*accum, p);
    packetCount_++;
  }

  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T initialize() const {
    internal::scalar_cast_op<int, T> conv;
    return conv(0);
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet initializePacket() const {
    return pset1<Packet>(initialize());
  }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T finalize(const T accum) const {
    internal::scalar_quotient_op<T> quotient_op;
    return quotient_op(accum, T(scalarCount_));
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet finalizePacket(const Packet& vaccum) const {
    return pdiv(vaccum, pset1<Packet>(T(packetCount_)));
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T finalizeBoth(const T saccum, const Packet& vaccum) const {
    internal::scalar_sum_op<T> sum_op;
    internal::scalar_quotient_op<T> quotient_op;
    return quotient_op(sum_op(saccum, predux(vaccum)), T(scalarCount_ + packetCount_ * unpacket_traits<Packet>::size));
  }

 protected:
  DenseIndex scalarCount_ = 0;
  DenseIndex packetCount_ = 0;
};

template <typename T, typename Device>
struct reducer_traits<MeanReducer<T>, Device> {
  enum {
    Cost = NumTraits<T>::AddCost,
    PacketAccess = PacketType<T, Device>::HasAdd && PacketType<T, Device>::HasDiv && !NumTraits<T>::IsInteger,
    IsStateful = true,
    IsExactlyAssociative = NumTraits<T>::IsInteger
  };
};

template <typename T, bool IsMax = true, bool IsInteger = true>
struct MinMaxBottomValue {
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T bottom_value() { return Eigen::NumTraits<T>::lowest(); }
};
template <typename T>
struct MinMaxBottomValue<T, true, false> {
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T bottom_value() { return -Eigen::NumTraits<T>::infinity(); }
};
template <typename T>
struct MinMaxBottomValue<T, false, true> {
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T bottom_value() { return Eigen::NumTraits<T>::highest(); }
};
template <typename T>
struct MinMaxBottomValue<T, false, false> {
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T bottom_value() { return Eigen::NumTraits<T>::infinity(); }
};

template <typename T, int NaNPropagation = PropagateFast>
struct MaxReducer {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const T t, T* accum) const {
    scalar_max_op<T, T, NaNPropagation> op;
    *accum = op(t, *accum);
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reducePacket(const Packet& p, Packet* accum) const {
    scalar_max_op<T, T, NaNPropagation> op;
    (*accum) = op.packetOp(*accum, p);
  }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T initialize() const {
    return MinMaxBottomValue<T, /*IsMax=*/true, Eigen::NumTraits<T>::IsInteger>::bottom_value();
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet initializePacket() const {
    return pset1<Packet>(initialize());
  }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T finalize(const T accum) const { return accum; }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet finalizePacket(const Packet& vaccum) const {
    return vaccum;
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T finalizeBoth(const T saccum, const Packet& vaccum) const {
    scalar_max_op<T, T, NaNPropagation> op;
    return op(saccum, op.predux(vaccum));
  }
};

template <typename T, typename Device, int NaNPropagation>
struct reducer_traits<MaxReducer<T, NaNPropagation>, Device> {
  enum {
    Cost = NumTraits<T>::AddCost,
    PacketAccess = PacketType<T, Device>::HasMax,
    IsStateful = false,
    IsExactlyAssociative = (NaNPropagation != PropagateFast)
  };
};

template <typename T, int NaNPropagation = PropagateFast>
struct MinReducer {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const T t, T* accum) const {
    scalar_min_op<T, T, NaNPropagation> op;
    *accum = op(t, *accum);
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reducePacket(const Packet& p, Packet* accum) const {
    scalar_min_op<T, T, NaNPropagation> op;
    (*accum) = op.packetOp(*accum, p);
  }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T initialize() const {
    return MinMaxBottomValue<T, /*IsMax=*/false, Eigen::NumTraits<T>::IsInteger>::bottom_value();
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet initializePacket() const {
    return pset1<Packet>(initialize());
  }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T finalize(const T accum) const { return accum; }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet finalizePacket(const Packet& vaccum) const {
    return vaccum;
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T finalizeBoth(const T saccum, const Packet& vaccum) const {
    scalar_min_op<T, T, NaNPropagation> op;
    return op(saccum, op.predux(vaccum));
  }
};

template <typename T, typename Device, int NaNPropagation>
struct reducer_traits<MinReducer<T, NaNPropagation>, Device> {
  enum {
    Cost = NumTraits<T>::AddCost,
    PacketAccess = PacketType<T, Device>::HasMin,
    IsStateful = false,
    IsExactlyAssociative = (NaNPropagation != PropagateFast)
  };
};

template <typename T>
struct ProdReducer {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const T t, T* accum) const {
    internal::scalar_product_op<T> prod_op;
    (*accum) = prod_op(*accum, t);
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reducePacket(const Packet& p, Packet* accum) const {
    (*accum) = pmul<Packet>(*accum, p);
  }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T initialize() const {
    internal::scalar_cast_op<int, T> conv;
    return conv(1);
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet initializePacket() const {
    return pset1<Packet>(initialize());
  }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T finalize(const T accum) const { return accum; }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet finalizePacket(const Packet& vaccum) const {
    return vaccum;
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T finalizeBoth(const T saccum, const Packet& vaccum) const {
    internal::scalar_product_op<T> prod_op;
    return prod_op(saccum, predux_mul(vaccum));
  }
};

template <typename T, typename Device>
struct reducer_traits<ProdReducer<T>, Device> {
  enum {
    Cost = NumTraits<T>::MulCost,
    PacketAccess = PacketType<T, Device>::HasMul,
    IsStateful = false,
    IsExactlyAssociative = true
  };
};

// Inputs have already been evaluated when reduce() is called; eager bitwise operations keep scalar reduction loops
// branch-free.
struct AndReducer {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(bool t, bool* accum) const { *accum = *accum & t; }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool initialize() const { return true; }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool finalize(bool accum) const { return accum; }
};

template <typename Device>
struct reducer_traits<AndReducer, Device> {
  enum { Cost = 1, PacketAccess = false, IsStateful = false, IsExactlyAssociative = true };
};

struct OrReducer {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(bool t, bool* accum) const { *accum = *accum | t; }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool initialize() const { return false; }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool finalize(bool accum) const { return accum; }
};

template <typename Device>
struct reducer_traits<OrReducer, Device> {
  enum { Cost = 1, PacketAccess = false, IsStateful = false, IsExactlyAssociative = true };
};

// Sum/prod reorder for scalars whose + and * commute; min/max additionally exclude custom scalars,
// whose generic std::min/std::max keep the first operand when values compare equivalent.
template <typename T>
struct reducer_can_reorder_accumulators<SumReducer<T>>
    : bool_constant<internal::is_arithmetic<T>::value || NumTraits<T>::IsComplex> {};
template <typename T>
struct reducer_can_reorder_accumulators<ProdReducer<T>>
    : bool_constant<internal::is_arithmetic<T>::value || NumTraits<T>::IsComplex> {};
template <typename T, int NaNPropagation>
struct reducer_can_reorder_accumulators<MinReducer<T, NaNPropagation>> : internal::is_arithmetic<T> {};
template <typename T, int NaNPropagation>
struct reducer_can_reorder_accumulators<MaxReducer<T, NaNPropagation>> : internal::is_arithmetic<T> {};
template <>
struct reducer_can_reorder_accumulators<AndReducer> : std::true_type {};
template <>
struct reducer_can_reorder_accumulators<OrReducer> : std::true_type {};

// Argmin/Argmax reducers.  Returns the first occurrence if multiple locations
// contain the same min/max value.
template <typename T>
struct ArgMaxPairReducer {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const T t, T* accum) const {
    if (t.second < accum->second) {
      return;
    } else if (t.second > accum->second || accum->first > t.first) {
      *accum = t;
    }
  }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T initialize() const {
    return T(0, NumTraits<typename T::second_type>::lowest());
  }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T finalize(const T& accum) const { return accum; }
};

template <typename T, typename Device>
struct reducer_traits<ArgMaxPairReducer<T>, Device> {
  enum { Cost = NumTraits<T>::AddCost, PacketAccess = false, IsStateful = false, IsExactlyAssociative = true };
};

template <typename T>
struct ArgMinPairReducer {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const T& t, T* accum) const {
    if (t.second > accum->second) {
      return;
    } else if (t.second < accum->second || accum->first > t.first) {
      *accum = t;
    }
  }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T initialize() const {
    return T(0, NumTraits<typename T::second_type>::highest());
  }
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T finalize(const T& accum) const { return accum; }
};

template <typename T, typename Device>
struct reducer_traits<ArgMinPairReducer<T>, Device> {
  enum { Cost = NumTraits<T>::AddCost, PacketAccess = false, IsStateful = false, IsExactlyAssociative = true };
};

// Both select the extremum of the lexicographic order on (value, -index), which is a total order for
// arithmetic values once NaNs are excluded; a scan never admits a NaN, since the index tie-break
// favors the accumulator, whose index is never larger. Custom scalars stay conservative.
template <typename T>
struct reducer_can_reorder_accumulators<ArgMaxPairReducer<T>> : internal::is_arithmetic<typename T::second_type> {};
template <typename T>
struct reducer_can_reorder_accumulators<ArgMinPairReducer<T>> : internal::is_arithmetic<typename T::second_type> {};

template <typename T, typename Index, size_t NumDims>
class GaussianGenerator {
 public:
  static constexpr bool PacketAccess = false;

  EIGEN_DEVICE_FUNC GaussianGenerator(const array<T, NumDims>& means, const array<T, NumDims>& std_devs)
      : m_means(means) {
    EIGEN_UNROLL_LOOP
    for (size_t i = 0; i < NumDims; ++i) {
      m_two_sigmas[i] = std_devs[i] * std_devs[i] * 2;
    }
  }

  EIGEN_DEVICE_FUNC T operator()(const array<Index, NumDims>& coordinates) const {
    T tmp = T(0);
    EIGEN_UNROLL_LOOP
    for (size_t i = 0; i < NumDims; ++i) {
      T offset = coordinates[i] - m_means[i];
      tmp += offset * offset / m_two_sigmas[i];
    }
    return numext::exp(-tmp);
  }

 private:
  array<T, NumDims> m_means;
  array<T, NumDims> m_two_sigmas;
};

template <typename T, typename Index, size_t NumDims>
struct functor_traits<GaussianGenerator<T, Index, NumDims>> {
  enum {
    Cost =
        NumDims * (2 * NumTraits<T>::AddCost + NumTraits<T>::MulCost + functor_traits<scalar_quotient_op<T, T>>::Cost) +
        functor_traits<scalar_exp_op<T>>::Cost,
    PacketAccess = GaussianGenerator<T, Index, NumDims>::PacketAccess
  };
};

template <typename Scalar>
struct scalar_clamp_op {
  EIGEN_DEVICE_FUNC inline scalar_clamp_op(const Scalar& _min, const Scalar& _max) : m_min(_min), m_max(_max) {}
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
    return numext::mini(numext::maxi(x, m_min), m_max);
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& x) const {
    return internal::pmin(internal::pmax(x, pset1<Packet>(m_min)), pset1<Packet>(m_max));
  }
  const Scalar m_min;
  const Scalar m_max;
};
template <typename Scalar>
struct functor_traits<scalar_clamp_op<Scalar>> {
  enum {
    Cost = 2 * NumTraits<Scalar>::AddCost,
    PacketAccess = (packet_traits<Scalar>::HasMin && packet_traits<Scalar>::HasMax)
  };
};

}  // end namespace internal
}  // end namespace Eigen

#endif  // EIGEN_TENSOR_TENSOR_FUNCTORS_H