GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

merge

Browse Source
This commit is contained in:
Gael Guennebaud 2016-01-28 13:21:48 +01:00
commit ddf64babde
9 changed files with 405 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
Eigen/src/Core/GenericPacketMath.h
View File

@ -75,6 +75,7 @@ struct default_packet_traits
HasCosh = 0,
HasTanh = 0,
HasLGamma = 0,
HasDiGamma = 0,
HasErf = 0,
HasErfc = 0,
@ -439,6 +440,10 @@ Packet pceil(const Packet& a) { using numext::ceil; return ceil(a); }
template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
Packet plgamma(const Packet& a) { using numext::lgamma; return lgamma(a); }
/** \internal \returns the derivative of lgamma, psi(\a a) (coeff-wise) */
template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
Packet pdigamma(const Packet& a) { using numext::digamma; return digamma(a); }
/** \internal \returns the erf(\a a) (coeff-wise) */
template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
Packet perf(const Packet& a) { using numext::erf; return erf(a); }

1
Eigen/src/Core/GlobalFunctions.h
View File

@ -50,6 +50,7 @@ namespace Eigen
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(cosh,scalar_cosh_op)
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(tanh,scalar_tanh_op)
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(lgamma,scalar_lgamma_op)
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(digamma,scalar_digamma_op)
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(erf,scalar_erf_op)
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(erfc,scalar_erfc_op)
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(exp,scalar_exp_op)

385
Eigen/src/Core/SpecialFunctions.h
View File

@ -13,79 +13,349 @@
namespace Eigen {
namespace internal {
// Parts of this code are based on the Cephes Math Library.
//
// Cephes Math Library Release 2.8: June, 2000
// Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
//
// Permission has been kindly provided by the original author
// to incorporate the Cephes software into the Eigen codebase:
//
// From: Stephen Moshier
// To: Eugene Brevdo
// Subject: Re: Permission to wrap several cephes functions in Eigen
//
// Hello Eugene,
//
// Thank you for writing.
//
// If your licensing is similar to BSD, the formal way that has been
// handled is simply to add a statement to the effect that you are incorporating
// the Cephes software by permission of the author.
//
// Good luck with your project,
// Steve
namespace cephes {
/* polevl (modified for Eigen)
*
* Evaluate polynomial
*
*
*
* SYNOPSIS:
*
* int N;
* Scalar x, y, coef[N+1];
*
* y = polevl<decltype(x), N>( x, coef);
*
*
*
* DESCRIPTION:
*
* Evaluates polynomial of degree N:
*
* 2 N
* y = C + C x + C x +...+ C x
* 0 1 2 N
*
* Coefficients are stored in reverse order:
*
* coef[0] = C , ..., coef[N] = C .
* N 0
*
* The function p1evl() assumes that coef[N] = 1.0 and is
* omitted from the array. Its calling arguments are
* otherwise the same as polevl().
*
*
* The Eigen implementation is templatized. For best speed, store
* coef as a const array (constexpr), e.g.
*
* const double coef[] = {1.0, 2.0, 3.0, ...};
*
*/
template <typename Scalar, int N>
struct polevl {
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
static Scalar run(const Scalar x, const Scalar coef[]) {
EIGEN_STATIC_ASSERT(N > 0, YOU_MADE_A_PROGRAMMING_MISTAKE);
return polevl<Scalar, N - 1>::run(x, coef) * x + coef[N];
}
};
template <typename Scalar>
struct polevl<Scalar, 0> {
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
static Scalar run(const Scalar, const Scalar coef[]) {
return coef[0];
}
};
} // end namespace cephes
/****************************************************************************
* Implementation of lgamma *
****************************************************************************/
template<typename Scalar>
struct lgamma_impl
{
template <typename Scalar>
struct lgamma_impl {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE Scalar run(const Scalar&)
{
static EIGEN_STRONG_INLINE Scalar run(const Scalar) {
EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
THIS_TYPE_IS_NOT_SUPPORTED);
return Scalar(0);
}
};
template<typename Scalar>
struct lgamma_retval
{
template <typename Scalar>
struct lgamma_retval {
typedef Scalar type;
};
#ifdef EIGEN_HAS_C99_MATH
template<>
struct lgamma_impl<float>
{
template <>
struct lgamma_impl<float> {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE double run(const float& x) { return ::lgammaf(x); }
static EIGEN_STRONG_INLINE float run(float x) { return ::lgammaf(x); }
};
template<>
struct lgamma_impl<double>
{
template <>
struct lgamma_impl<double> {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE double run(const double& x) { return ::lgamma(x); }
static EIGEN_STRONG_INLINE double run(double x) { return ::lgamma(x); }
};
#endif
/****************************************************************************
* Implementation of digamma (psi) *
****************************************************************************/
#ifdef EIGEN_HAS_C99_MATH
/*
*
* Polynomial evaluation helper for the Psi (digamma) function.
*
* digamma_impl_maybe_poly::run(s) evaluates the asymptotic Psi expansion for
* input Scalar s, assuming s is above 10.0.
*
* If s is above a certain threshold for the given Scalar type, zero
* is returned. Otherwise the polynomial is evaluated with enough
* coefficients for results matching Scalar machine precision.
*
*
*/
template <typename Scalar>
struct digamma_impl_maybe_poly {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE Scalar run(const Scalar) {
EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
THIS_TYPE_IS_NOT_SUPPORTED);
return Scalar(0);
}
};
template <>
struct digamma_impl_maybe_poly<float> {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE float run(const float s) {
const float A[] = {
-4.16666666666666666667E-3,
3.96825396825396825397E-3,
-8.33333333333333333333E-3,
8.33333333333333333333E-2
};
float z;
if (s < 1.0e8f) {
z = 1.0f / (s * s);
return z * cephes::polevl<float, 3>::run(z, A);
} else return 0.0f;
}
};
template <>
struct digamma_impl_maybe_poly<double> {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE double run(const double s) {
const double A[] = {
8.33333333333333333333E-2,
-2.10927960927960927961E-2,
7.57575757575757575758E-3,
-4.16666666666666666667E-3,
3.96825396825396825397E-3,
-8.33333333333333333333E-3,
8.33333333333333333333E-2
};
double z;
if (s < 1.0e17) {
z = 1.0 / (s * s);
return z * cephes::polevl<double, 6>::run(z, A);
}
else return 0.0;
}
};
#endif // EIGEN_HAS_C99_MATH
template <typename Scalar>
struct digamma_retval {
typedef Scalar type;
};
#ifdef EIGEN_HAS_C99_MATH
template <typename Scalar>
struct digamma_impl {
EIGEN_DEVICE_FUNC
static Scalar run(Scalar x) {
/*
*
* Psi (digamma) function (modified for Eigen)
*
*
* SYNOPSIS:
*
* double x, y, psi();
*
* y = psi( x );
*
*
* DESCRIPTION:
*
* d -
* psi(x) = -- ln | (x)
* dx
*
* is the logarithmic derivative of the gamma function.
* For integer x,
* n-1
* -
* psi(n) = -EUL + > 1/k.
* -
* k=1
*
* If x is negative, it is transformed to a positive argument by the
* reflection formula psi(1-x) = psi(x) + pi cot(pi x).
* For general positive x, the argument is made greater than 10
* using the recurrence psi(x+1) = psi(x) + 1/x.
* Then the following asymptotic expansion is applied:
*
* inf. B
* - 2k
* psi(x) = log(x) - 1/2x - > -------
* - 2k
* k=1 2k x
*
* where the B2k are Bernoulli numbers.
*
* ACCURACY (float):
* Relative error (except absolute when |psi| < 1):
* arithmetic domain # trials peak rms
* IEEE 0,30 30000 1.3e-15 1.4e-16
* IEEE -30,0 40000 1.5e-15 2.2e-16
*
* ACCURACY (double):
* Absolute error, relative when |psi| > 1 :
* arithmetic domain # trials peak rms
* IEEE -33,0 30000 8.2e-7 1.2e-7
* IEEE 0,33 100000 7.3e-7 7.7e-8
*
* ERROR MESSAGES:
* message condition value returned
* psi singularity x integer <=0 INFINITY
*/
Scalar p, q, nz, s, w, y;
bool negative;
const Scalar maxnum = std::numeric_limits<Scalar>::infinity();
const Scalar m_pi = 3.14159265358979323846;
negative = 0;
nz = 0.0;
const Scalar zero = 0.0;
const Scalar one = 1.0;
const Scalar half = 0.5;
if (x <= zero) {
negative = one;
q = x;
p = ::floor(q);
if (p == q) {
return maxnum;
}
/* Remove the zeros of tan(m_pi x)
* by subtracting the nearest integer from x
*/
nz = q - p;
if (nz != half) {
if (nz > half) {
p += one;
nz = q - p;
}
nz = m_pi / ::tan(m_pi * nz);
}
else {
nz = zero;
}
x = one - x;
}
/* use the recurrence psi(x+1) = psi(x) + 1/x. */
s = x;
w = zero;
while (s < Scalar(10)) {
w += one / s;
s += one;
}
y = digamma_impl_maybe_poly<Scalar>::run(s);
y = ::log(s) - (half / s) - y - w;
return (negative) ? y - nz : y;
}
};
#endif // EIGEN_HAS_C99_MATH
/****************************************************************************
* Implementation of erf *
****************************************************************************/
template<typename Scalar>
struct erf_impl
{
template <typename Scalar>
struct erf_impl {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE Scalar run(const Scalar&)
{
static EIGEN_STRONG_INLINE Scalar run(const Scalar) {
EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
THIS_TYPE_IS_NOT_SUPPORTED);
return Scalar(0);
}
};
template<typename Scalar>
struct erf_retval
{
template <typename Scalar>
struct erf_retval {
typedef Scalar type;
};
#ifdef EIGEN_HAS_C99_MATH
template<>
struct erf_impl<float>
{
template <>
struct erf_impl<float> {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE float run(const float& x) { return ::erff(x); }
static EIGEN_STRONG_INLINE float run(float x) { return ::erff(x); }
};
template<>
struct erf_impl<double>
{
template <>
struct erf_impl<double> {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE double run(const double& x) { return ::erf(x); }
static EIGEN_STRONG_INLINE double run(double x) { return ::erf(x); }
};
#endif // EIGEN_HAS_C99_MATH
@ -93,35 +363,30 @@ struct erf_impl<double>
* Implementation of erfc *
****************************************************************************/
template<typename Scalar>
struct erfc_impl
{
template <typename Scalar>
struct erfc_impl {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE Scalar run(const Scalar&)
{
static EIGEN_STRONG_INLINE Scalar run(const Scalar) {
EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
THIS_TYPE_IS_NOT_SUPPORTED);
return Scalar(0);
}
};
template<typename Scalar>
struct erfc_retval
{
template <typename Scalar>
struct erfc_retval {
typedef Scalar type;
};
#ifdef EIGEN_HAS_C99_MATH
template<>
struct erfc_impl<float>
{
template <>
struct erfc_impl<float> {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE float run(const float x) { return ::erfcf(x); }
};
template<>
struct erfc_impl<double>
{
template <>
struct erfc_impl<double> {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE double run(const double x) { return ::erfc(x); }
};
@ -129,27 +394,29 @@ struct erfc_impl<double>
} // end namespace internal
namespace numext {
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(lgamma, Scalar) lgamma(const Scalar& x)
{
template <typename Scalar>
EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(lgamma, Scalar)
lgamma(const Scalar& x) {
return EIGEN_MATHFUNC_IMPL(lgamma, Scalar)::run(x);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(erf, Scalar) erf(const Scalar& x)
{
template <typename Scalar>
EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(digamma, Scalar)
digamma(const Scalar& x) {
return EIGEN_MATHFUNC_IMPL(digamma, Scalar)::run(x);
}
template <typename Scalar>
EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erf, Scalar)
erf(const Scalar& x) {
return EIGEN_MATHFUNC_IMPL(erf, Scalar)::run(x);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(erfc, Scalar) erfc(const Scalar& x)
{
template <typename Scalar>
EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erfc, Scalar)
erfc(const Scalar& x) {
return EIGEN_MATHFUNC_IMPL(erfc, Scalar)::run(x);
}

14
Eigen/src/Core/arch/CUDA/MathFunctions.h
View File

@ -78,6 +78,20 @@ double2 plgamma<double2>(const double2& a)
return make_double2(lgamma(a.x), lgamma(a.y));
}
template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
float4 pdigamma<float4>(const float4& a)
{
using numext::digamma;
return make_float4(digamma(a.x), digamma(a.y), digamma(a.z), digamma(a.w));
}
template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
double2 pdigamma<double2>(const double2& a)
{
using numext::digamma;
return make_double2(digamma(a.x), digamma(a.y));
}
template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
float4 perf<float4>(const float4& a)
{

2
Eigen/src/Core/arch/CUDA/PacketMath.h
View File

@ -40,6 +40,7 @@ template<> struct packet_traits<float> : default_packet_traits
HasSqrt = 1,
HasRsqrt = 1,
HasLGamma = 1,
HasDiGamma = 1,
HasErf = 1,
HasErfc = 1,
@ -63,6 +64,7 @@ template<> struct packet_traits<double> : default_packet_traits
HasSqrt = 1,
HasRsqrt = 1,
HasLGamma = 1,
HasDiGamma = 1,
HasErf = 1,
HasErfc = 1,

22
Eigen/src/Core/functors/UnaryFunctors.h
View File

@ -427,6 +427,28 @@ struct functor_traits<scalar_lgamma_op<Scalar> >
};
};
/** \internal
* \brief Template functor to compute psi, the derivative of lgamma of a scalar.
* \sa class CwiseUnaryOp, Cwise::digamma()
*/
template<typename Scalar> struct scalar_digamma_op {
EIGEN_EMPTY_STRUCT_CTOR(scalar_digamma_op)
EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& a) const {
using numext::digamma; return digamma(a);
}
typedef typename packet_traits<Scalar>::type Packet;
EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& a) const { return internal::pdigamma(a); }
};
template<typename Scalar>
struct functor_traits<scalar_digamma_op<Scalar> >
{
enum {
// Guesstimate
Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
PacketAccess = packet_traits<Scalar>::HasDiGamma
};
};
/** \internal
* \brief Template functor to compute the Gauss error function of a
* scalar

11
Eigen/src/plugins/ArrayCwiseUnaryOps.h
View File

@ -22,6 +22,7 @@ typedef CwiseUnaryOp<internal::scalar_tanh_op<Scalar>, const Derived> TanhReturn
typedef CwiseUnaryOp<internal::scalar_sinh_op<Scalar>, const Derived> SinhReturnType;
typedef CwiseUnaryOp<internal::scalar_cosh_op<Scalar>, const Derived> CoshReturnType;
typedef CwiseUnaryOp<internal::scalar_lgamma_op<Scalar>, const Derived> LgammaReturnType;
typedef CwiseUnaryOp<internal::scalar_digamma_op<Scalar>, const Derived> DigammaReturnType;
typedef CwiseUnaryOp<internal::scalar_erf_op<Scalar>, const Derived> ErfReturnType;
typedef CwiseUnaryOp<internal::scalar_erfc_op<Scalar>, const Derived> ErfcReturnType;
typedef CwiseUnaryOp<internal::scalar_pow_op<Scalar>, const Derived> PowReturnType;
@ -318,6 +319,16 @@ lgamma() const
return LgammaReturnType(derived());
}
/** \returns an expression of the coefficient-wise digamma (psi, derivative of lgamma).
*
* \sa cos(), sin(), tan()
*/
inline const DigammaReturnType
digamma() const
{
return DigammaReturnType(derived());
}
/** \returns an expression of the coefficient-wise Gauss error
* function of *this.
*

19
test/array.cpp
View File

@ -219,6 +219,7 @@ template<typename ArrayType> void array_real(const ArrayType& m)
VERIFY_IS_APPROX(m1.tanh(), tanh(m1));
#ifdef EIGEN_HAS_C99_MATH
VERIFY_IS_APPROX(m1.lgamma(), lgamma(m1));
VERIFY_IS_APPROX(m1.digamma(), digamma(m1));
VERIFY_IS_APPROX(m1.erf(), erf(m1));
VERIFY_IS_APPROX(m1.erfc(), erfc(m1));
#endif // EIGEN_HAS_C99_MATH
@ -310,6 +311,21 @@ template<typename ArrayType> void array_real(const ArrayType& m)
m1 += ArrayType::Constant(rows,cols,Scalar(tiny));
VERIFY_IS_APPROX(s1/m1, s1 * m1.inverse());
// check special functions (comparing against numpy implementation)
#ifdef EIGEN_HAS_C99_MATH
if (!NumTraits<Scalar>::IsComplex) {
VERIFY_IS_APPROX(numext::digamma(Scalar(1)), RealScalar(-0.5772156649015329));
VERIFY_IS_APPROX(numext::digamma(Scalar(1.5)), RealScalar(0.03648997397857645));
VERIFY_IS_APPROX(numext::digamma(Scalar(4)), RealScalar(1.2561176684318));
VERIFY_IS_APPROX(numext::digamma(Scalar(-10.5)), RealScalar(2.398239129535781));
VERIFY_IS_APPROX(numext::digamma(Scalar(10000.5)), RealScalar(9.210340372392849));
VERIFY_IS_EQUAL(numext::digamma(Scalar(0)),
std::numeric_limits<RealScalar>::infinity());
VERIFY_IS_EQUAL(numext::digamma(Scalar(-1)),
std::numeric_limits<RealScalar>::infinity());
}
#endif // EIGEN_HAS_C99_MATH
// check inplace transpose
m3 = m1;
m3.transposeInPlace();
@ -336,8 +352,6 @@ template<typename ArrayType> void array_complex(const ArrayType& m)
Array<RealScalar, -1, -1> m3(rows, cols);
Scalar s1 = internal::random<Scalar>();
for (Index i = 0; i < m.rows(); ++i)
for (Index j = 0; j < m.cols(); ++j)
m2(i,j) = sqrt(m1(i,j));
@ -410,6 +424,7 @@ template<typename ArrayType> void array_complex(const ArrayType& m)
VERIFY_IS_APPROX( m1.sign() * m1.abs(), m1);
// scalar by array division
Scalar s1 = internal::random<Scalar>();
const RealScalar tiny = sqrt(std::numeric_limits<RealScalar>::epsilon());
s1 += Scalar(tiny);
m1 += ArrayType::Constant(rows,cols,Scalar(tiny));

6
unsupported/Eigen/CXX11/src/Tensor/TensorBase.h
View File

@ -128,6 +128,12 @@ class TensorBase<Derived, ReadOnlyAccessors>
return unaryExpr(internal::scalar_lgamma_op<Scalar>());
}
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE const TensorCwiseUnaryOp<internal::scalar_digamma_op<Scalar>, const Derived>
digamma() const {
return unaryExpr(internal::scalar_digamma_op<Scalar>());
}
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE const TensorCwiseUnaryOp<internal::scalar_erf_op<Scalar>, const Derived>
erf() const {