diff options
author | Deven Desai <deven.desai.amd@gmail.com> | 2018-06-13 12:09:52 -0400 |
---|---|---|
committer | Deven Desai <deven.desai.amd@gmail.com> | 2018-06-13 12:09:52 -0400 |
commit | d1d22ef0f4af42f58bdd9d78b22bf912852a6bf4 (patch) | |
tree | d137f1e11d54028c241eee61bf8cd5fe6441f602 /unsupported/Eigen/src | |
parent | 8fbd47052bcafea612b8ae2841c1de5db738f042 (diff) | |
parent | d3a380af4d17513ab71630b59f390589fa7c207b (diff) |
syncing this fork with upstream
Diffstat (limited to 'unsupported/Eigen/src')
23 files changed, 574 insertions, 264 deletions
diff --git a/unsupported/Eigen/src/BVH/KdBVH.h b/unsupported/Eigen/src/BVH/KdBVH.h index 1b8d75865..13f792cd0 100644 --- a/unsupported/Eigen/src/BVH/KdBVH.h +++ b/unsupported/Eigen/src/BVH/KdBVH.h @@ -170,7 +170,7 @@ private: typedef internal::vector_int_pair<Scalar, Dim> VIPair; typedef std::vector<VIPair, aligned_allocator<VIPair> > VIPairList; typedef Matrix<Scalar, Dim, 1> VectorType; - struct VectorComparator //compares vectors, or, more specificall, VIPairs along a particular dimension + struct VectorComparator //compares vectors, or more specifically, VIPairs along a particular dimension { VectorComparator(int inDim) : dim(inDim) {} inline bool operator()(const VIPair &v1, const VIPair &v2) const { return v1.first[dim] < v2.first[dim]; } diff --git a/unsupported/Eigen/src/Eigenvalues/ArpackSelfAdjointEigenSolver.h b/unsupported/Eigen/src/Eigenvalues/ArpackSelfAdjointEigenSolver.h index 866a8a460..9f7bff764 100644 --- a/unsupported/Eigen/src/Eigenvalues/ArpackSelfAdjointEigenSolver.h +++ b/unsupported/Eigen/src/Eigenvalues/ArpackSelfAdjointEigenSolver.h @@ -300,7 +300,7 @@ public: /** \brief Reports whether previous computation was successful. * - * \returns \c Success if computation was succesful, \c NoConvergence otherwise. + * \returns \c Success if computation was successful, \c NoConvergence otherwise. */ ComputationInfo info() const { diff --git a/unsupported/Eigen/src/EulerAngles/EulerSystem.h b/unsupported/Eigen/src/EulerAngles/EulerSystem.h index 28f52da61..65c2e94c7 100644 --- a/unsupported/Eigen/src/EulerAngles/EulerSystem.h +++ b/unsupported/Eigen/src/EulerAngles/EulerSystem.h @@ -12,7 +12,7 @@ namespace Eigen { - // Forward declerations + // Forward declarations template <typename _Scalar, class _System> class EulerAngles; diff --git a/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h b/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h index dc0093eb9..37d5b4c6c 100644 --- a/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h +++ b/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h @@ -99,7 +99,7 @@ void pseudo_inverse(const CMatrix &C, CINVMatrix &CINV) /** \ingroup IterativeSolvers_Module * Constrained conjugate gradient * - * Computes the minimum of \f$ 1/2((Ax).x) - bx \f$ under the contraint \f$ Cx \le f \f$ + * Computes the minimum of \f$ 1/2((Ax).x) - bx \f$ under the constraint \f$ Cx \le f \f$ */ template<typename TMatrix, typename CMatrix, typename VectorX, typename VectorB, typename VectorF> diff --git a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h index d603ba336..be039e07f 100644 --- a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h @@ -39,7 +39,6 @@ template <typename VectorType, typename IndexType> void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::Scalar& ncut) { eigen_assert(vec.size() == perm.size()); - typedef typename IndexType::Scalar Index; bool flag; for (Index k = 0; k < ncut; k++) { @@ -112,7 +111,6 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > using Base::_solve_impl; typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::Index Index; typedef typename MatrixType::StorageIndex StorageIndex; typedef typename MatrixType::RealScalar RealScalar; typedef _Preconditioner Preconditioner; @@ -146,7 +144,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > void _solve_with_guess_impl(const Rhs& b, Dest& x) const { bool failed = false; - for(int j=0; j<b.cols(); ++j) + for(Index j=0; j<b.cols(); ++j) { m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; @@ -170,17 +168,17 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > /** * Get the restart value */ - int restart() { return m_restart; } + Index restart() { return m_restart; } /** * Set the restart value (default is 30) */ - void set_restart(const int restart) { m_restart=restart; } + void set_restart(const Index restart) { m_restart=restart; } /** * Set the number of eigenvalues to deflate at each restart */ - void setEigenv(const int neig) + void setEigenv(const Index neig) { m_neig = neig; if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates @@ -189,12 +187,12 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > /** * Get the size of the deflation subspace size */ - int deflSize() {return m_r; } + Index deflSize() {return m_r; } /** * Set the maximum size of the deflation subspace */ - void setMaxEigenv(const int maxNeig) { m_maxNeig = maxNeig; } + void setMaxEigenv(const Index maxNeig) { m_maxNeig = maxNeig; } protected: // DGMRES algorithm @@ -202,27 +200,27 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const; // Perform one cycle of GMRES template<typename Dest> - int dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const; + Index dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const; // Compute data to use for deflation - int dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const; + Index dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const; // Apply deflation to a vector template<typename RhsType, typename DestType> - int dgmresApplyDeflation(const RhsType& In, DestType& Out) const; + Index dgmresApplyDeflation(const RhsType& In, DestType& Out) const; ComplexVector schurValues(const ComplexSchur<DenseMatrix>& schurofH) const; ComplexVector schurValues(const RealSchur<DenseMatrix>& schurofH) const; // Init data for deflation void dgmresInitDeflation(Index& rows) const; mutable DenseMatrix m_V; // Krylov basis vectors mutable DenseMatrix m_H; // Hessenberg matrix - mutable DenseMatrix m_Hes; // Initial hessenberg matrix wihout Givens rotations applied + mutable DenseMatrix m_Hes; // Initial hessenberg matrix without Givens rotations applied mutable Index m_restart; // Maximum size of the Krylov subspace mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles) mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */ mutable PartialPivLU<DenseMatrix> m_luT; // LU factorization of m_T mutable StorageIndex m_neig; //Number of eigenvalues to extract at each restart - mutable int m_r; // Current number of deflated eigenvalues, size of m_U - mutable int m_maxNeig; // Maximum number of eigenvalues to deflate + mutable Index m_r; // Current number of deflated eigenvalues, size of m_U + mutable Index m_maxNeig; // Maximum number of eigenvalues to deflate mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A mutable bool m_isDeflAllocated; mutable bool m_isDeflInitialized; @@ -244,13 +242,13 @@ void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rh const Preconditioner& precond) const { //Initialization - int n = mat.rows(); + Index n = mat.rows(); DenseVector r0(n); - int nbIts = 0; + Index nbIts = 0; m_H.resize(m_restart+1, m_restart); m_Hes.resize(m_restart, m_restart); m_V.resize(n,m_restart+1); - //Initial residual vector and intial norm + //Initial residual vector and initial norm x = precond.solve(x); r0 = rhs - mat * x; RealScalar beta = r0.norm(); @@ -284,7 +282,7 @@ void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rh */ template< typename _MatrixType, typename _Preconditioner> template<typename Dest> -int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const +Index DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const { //Initialization DenseVector g(m_restart+1); // Right hand side of the least square problem @@ -293,8 +291,8 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con m_V.col(0) = r0/beta; m_info = NoConvergence; std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations - int it = 0; // Number of inner iterations - int n = mat.rows(); + Index it = 0; // Number of inner iterations + Index n = mat.rows(); DenseVector tv1(n), tv2(n); //Temporary vectors while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations) { @@ -312,7 +310,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con // Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt Scalar coef; - for (int i = 0; i <= it; ++i) + for (Index i = 0; i <= it; ++i) { coef = tv1.dot(m_V.col(i)); tv1 = tv1 - coef * m_V.col(i); @@ -328,7 +326,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con // FIXME Check for happy breakdown // Update Hessenberg matrix with Givens rotations - for (int i = 1; i <= it; ++i) + for (Index i = 1; i <= it; ++i) { m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint()); } @@ -418,7 +416,7 @@ inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_Matr } template< typename _MatrixType, typename _Preconditioner> -int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const +Index DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const { // First, find the Schur form of the Hessenberg matrix H typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::type schurofH; @@ -433,8 +431,8 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri // Reorder the absolute values of Schur values DenseRealVector modulEig(it); - for (int j=0; j<it; ++j) modulEig(j) = std::abs(eig(j)); - perm.setLinSpaced(it,0,it-1); + for (Index j=0; j<it; ++j) modulEig(j) = std::abs(eig(j)); + perm.setLinSpaced(it,0,internal::convert_index<StorageIndex>(it-1)); internal::sortWithPermutation(modulEig, perm, neig); if (!m_lambdaN) @@ -442,7 +440,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN); } //Count the real number of extracted eigenvalues (with complex conjugates) - int nbrEig = 0; + Index nbrEig = 0; while (nbrEig < neig) { if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++; @@ -451,7 +449,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri // Extract the Schur vectors corresponding to the smallest Ritz values DenseMatrix Sr(it, nbrEig); Sr.setZero(); - for (int j = 0; j < nbrEig; j++) + for (Index j = 0; j < nbrEig; j++) { Sr.col(j) = schurofH.matrixU().col(perm(it-j-1)); } @@ -462,8 +460,8 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri if (m_r) { // Orthogonalize X against m_U using modified Gram-Schmidt - for (int j = 0; j < nbrEig; j++) - for (int k =0; k < m_r; k++) + for (Index j = 0; j < nbrEig; j++) + for (Index k =0; k < m_r; k++) X.col(j) = X.col(j) - (m_U.col(k).dot(X.col(j)))*m_U.col(k); } @@ -473,7 +471,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri dgmresInitDeflation(m); DenseMatrix MX(m, nbrEig); DenseVector tv1(m); - for (int j = 0; j < nbrEig; j++) + for (Index j = 0; j < nbrEig; j++) { tv1 = mat * X.col(j); MX.col(j) = precond.solve(tv1); @@ -488,8 +486,8 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri } // Save X into m_U and m_MX in m_MU - for (int j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j); - for (int j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j); + for (Index j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j); + for (Index j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j); // Increase the size of the invariant subspace m_r += nbrEig; @@ -502,7 +500,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri } template<typename _MatrixType, typename _Preconditioner> template<typename RhsType, typename DestType> -int DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const +Index DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const { DenseVector x1 = m_U.leftCols(m_r).transpose() * x; y = x + m_U.leftCols(m_r) * ( m_lambdaN * m_luT.solve(x1) - x1); diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LMqrsolv.h b/unsupported/Eigen/src/LevenbergMarquardt/LMqrsolv.h index ae9d793b1..123485817 100644 --- a/unsupported/Eigen/src/LevenbergMarquardt/LMqrsolv.h +++ b/unsupported/Eigen/src/LevenbergMarquardt/LMqrsolv.h @@ -73,7 +73,7 @@ void lmqrsolv( qtbpj = -givens.s() * wa[k] + givens.c() * qtbpj; wa[k] = temp; - /* accumulate the tranformation in the row of s. */ + /* accumulate the transformation in the row of s. */ for (i = k+1; i<n; ++i) { temp = givens.c() * s(i,k) + givens.s() * sdiag[i]; sdiag[i] = -givens.s() * s(i,k) + givens.c() * sdiag[i]; diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h index 995427978..5f64501be 100644 --- a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h +++ b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h @@ -233,9 +233,9 @@ class LevenbergMarquardt : internal::no_assignment_operator /** * \brief Reports whether the minimization was successful - * \returns \c Success if the minimization was succesful, + * \returns \c Success if the minimization was successful, * \c NumericalIssue if a numerical problem arises during the - * minimization process, for exemple during the QR factorization + * minimization process, for example during the QR factorization * \c NoConvergence if the minimization did not converge after * the maximum number of function evaluation allowed * \c InvalidInput if the input matrix is invalid diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index 85ab3d97c..03356998b 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -313,7 +313,7 @@ struct matrix_exp_computeUV<MatrixType, long double> matrix_exp_pade17(A, U, V); } -#elif LDBL_MANT_DIG <= 112 // quadruple precison +#elif LDBL_MANT_DIG <= 112 // quadruple precision if (l1norm < 1.639394610288918690547467954466970e-005L) { matrix_exp_pade3(arg, U, V); diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index ebc433d89..33609aea9 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -81,7 +81,7 @@ class MatrixPowerParenthesesReturnValue : public ReturnByValue< MatrixPowerParen * * \note Currently this class is only used by MatrixPower. One may * insist that this be nested into MatrixPower. This class is here to - * faciliate future development of triangular matrix functions. + * facilitate future development of triangular matrix functions. */ template<typename MatrixType> class MatrixPowerAtomic : internal::noncopyable diff --git a/unsupported/Eigen/src/NonLinearOptimization/qrsolv.h b/unsupported/Eigen/src/NonLinearOptimization/qrsolv.h index feafd62a8..4f2f560b3 100644 --- a/unsupported/Eigen/src/NonLinearOptimization/qrsolv.h +++ b/unsupported/Eigen/src/NonLinearOptimization/qrsolv.h @@ -61,7 +61,7 @@ void qrsolv( qtbpj = -givens.s() * wa[k] + givens.c() * qtbpj; wa[k] = temp; - /* accumulate the tranformation in the row of s. */ + /* accumulate the transformation in the row of s. */ for (i = k+1; i<n; ++i) { temp = givens.c() * s(i,k) + givens.s() * sdiag[i]; sdiag[i] = -givens.s() * s(i,k) + givens.c() * sdiag[i]; diff --git a/unsupported/Eigen/src/NonLinearOptimization/r1updt.h b/unsupported/Eigen/src/NonLinearOptimization/r1updt.h index f28766061..09fc65255 100644 --- a/unsupported/Eigen/src/NonLinearOptimization/r1updt.h +++ b/unsupported/Eigen/src/NonLinearOptimization/r1updt.h @@ -22,7 +22,7 @@ void r1updt( Scalar temp; JacobiRotation<Scalar> givens; - // r1updt had a broader usecase, but we dont use it here. And, more + // r1updt had a broader usecase, but we don't use it here. And, more // importantly, we can not test it. eigen_assert(m==n); eigen_assert(u.size()==m); diff --git a/unsupported/Eigen/src/Polynomials/Companion.h b/unsupported/Eigen/src/Polynomials/Companion.h index e0af6ebe4..41a4efc2f 100644 --- a/unsupported/Eigen/src/Polynomials/Companion.h +++ b/unsupported/Eigen/src/Polynomials/Companion.h @@ -104,7 +104,7 @@ class companion /** Helper function for the balancing algorithm. * \returns true if the row and the column, having colNorm and rowNorm * as norms, are balanced, false otherwise. - * colB and rowB are repectively the multipliers for + * colB and rowB are respectively the multipliers for * the column and the row in order to balance them. * */ bool balanced( RealScalar colNorm, RealScalar rowNorm, @@ -113,7 +113,7 @@ class companion /** Helper function for the balancing algorithm. * \returns true if the row and the column, having colNorm and rowNorm * as norms, are balanced, false otherwise. - * colB and rowB are repectively the multipliers for + * colB and rowB are respectively the multipliers for * the column and the row in order to balance them. * */ bool balancedR( RealScalar colNorm, RealScalar rowNorm, diff --git a/unsupported/Eigen/src/Skyline/SkylineInplaceLU.h b/unsupported/Eigen/src/Skyline/SkylineInplaceLU.h index a1f54ed35..bda057a85 100644 --- a/unsupported/Eigen/src/Skyline/SkylineInplaceLU.h +++ b/unsupported/Eigen/src/Skyline/SkylineInplaceLU.h @@ -41,7 +41,7 @@ public: /** Sets the relative threshold value used to prune zero coefficients during the decomposition. * - * Setting a value greater than zero speeds up computation, and yields to an imcomplete + * Setting a value greater than zero speeds up computation, and yields to an incomplete * factorization with fewer non zero coefficients. Such approximate factors are especially * useful to initialize an iterative solver. * diff --git a/unsupported/Eigen/src/Skyline/SkylineMatrix.h b/unsupported/Eigen/src/Skyline/SkylineMatrix.h index a2a8933ca..f77d79a04 100644 --- a/unsupported/Eigen/src/Skyline/SkylineMatrix.h +++ b/unsupported/Eigen/src/Skyline/SkylineMatrix.h @@ -206,26 +206,26 @@ public: if (col > row) //upper matrix { const Index minOuterIndex = inner - m_data.upperProfile(inner); - eigen_assert(outer >= minOuterIndex && "you try to acces a coeff that do not exist in the storage"); + eigen_assert(outer >= minOuterIndex && "You tried to access a coeff that does not exist in the storage"); return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner))); } if (col < row) //lower matrix { const Index minInnerIndex = outer - m_data.lowerProfile(outer); - eigen_assert(inner >= minInnerIndex && "you try to acces a coeff that do not exist in the storage"); + eigen_assert(inner >= minInnerIndex && "You tried to access a coeff that does not exist in the storage"); return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer))); } } else { if (outer > inner) //upper matrix { const Index maxOuterIndex = inner + m_data.upperProfile(inner); - eigen_assert(outer <= maxOuterIndex && "you try to acces a coeff that do not exist in the storage"); + eigen_assert(outer <= maxOuterIndex && "You tried to access a coeff that does not exist in the storage"); return this->m_data.upper(m_colStartIndex[inner] + (outer - inner)); } if (outer < inner) //lower matrix { const Index maxInnerIndex = outer + m_data.lowerProfile(outer); - eigen_assert(inner <= maxInnerIndex && "you try to acces a coeff that do not exist in the storage"); + eigen_assert(inner <= maxInnerIndex && "You tried to access a coeff that does not exist in the storage"); return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer)); } } @@ -300,11 +300,11 @@ public: if (IsRowMajor) { const Index minInnerIndex = outer - m_data.lowerProfile(outer); - eigen_assert(inner >= minInnerIndex && "you try to acces a coeff that do not exist in the storage"); + eigen_assert(inner >= minInnerIndex && "You tried to access a coeff that does not exist in the storage"); return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer))); } else { const Index maxInnerIndex = outer + m_data.lowerProfile(outer); - eigen_assert(inner <= maxInnerIndex && "you try to acces a coeff that do not exist in the storage"); + eigen_assert(inner <= maxInnerIndex && "You tried to access a coeff that does not exist in the storage"); return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer)); } } @@ -336,11 +336,11 @@ public: if (IsRowMajor) { const Index minOuterIndex = inner - m_data.upperProfile(inner); - eigen_assert(outer >= minOuterIndex && "you try to acces a coeff that do not exist in the storage"); + eigen_assert(outer >= minOuterIndex && "You tried to access a coeff that does not exist in the storage"); return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner))); } else { const Index maxOuterIndex = inner + m_data.upperProfile(inner); - eigen_assert(outer <= maxOuterIndex && "you try to acces a coeff that do not exist in the storage"); + eigen_assert(outer <= maxOuterIndex && "You tried to access a coeff that does not exist in the storage"); return this->m_data.upper(m_colStartIndex[inner] + (outer - inner)); } } diff --git a/unsupported/Eigen/src/SparseExtra/DynamicSparseMatrix.h b/unsupported/Eigen/src/SparseExtra/DynamicSparseMatrix.h index 037a13f86..3f1ff14ad 100644 --- a/unsupported/Eigen/src/SparseExtra/DynamicSparseMatrix.h +++ b/unsupported/Eigen/src/SparseExtra/DynamicSparseMatrix.h @@ -187,7 +187,7 @@ template<typename _Scalar, int _Options, typename _StorageIndex> /** Does nothing: provided for compatibility with SparseMatrix */ inline void finalize() {} - /** Suppress all nonzeros which are smaller than \a reference under the tolerence \a epsilon */ + /** Suppress all nonzeros which are smaller than \a reference under the tolerance \a epsilon */ void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision()) { for (Index j=0; j<outerSize(); ++j) @@ -224,21 +224,21 @@ template<typename _Scalar, int _Options, typename _StorageIndex> } } - /** The class DynamicSparseMatrix is deprectaed */ + /** The class DynamicSparseMatrix is deprecated */ EIGEN_DEPRECATED inline DynamicSparseMatrix() : m_innerSize(0), m_data(0) { eigen_assert(innerSize()==0 && outerSize()==0); } - /** The class DynamicSparseMatrix is deprectaed */ + /** The class DynamicSparseMatrix is deprecated */ EIGEN_DEPRECATED inline DynamicSparseMatrix(Index rows, Index cols) : m_innerSize(0) { resize(rows, cols); } - /** The class DynamicSparseMatrix is deprectaed */ + /** The class DynamicSparseMatrix is deprecated */ template<typename OtherDerived> EIGEN_DEPRECATED explicit inline DynamicSparseMatrix(const SparseMatrixBase<OtherDerived>& other) : m_innerSize(0) diff --git a/unsupported/Eigen/src/SparseExtra/MarketIO.h b/unsupported/Eigen/src/SparseExtra/MarketIO.h index 6d57ab2e9..1618b09a8 100644 --- a/unsupported/Eigen/src/SparseExtra/MarketIO.h +++ b/unsupported/Eigen/src/SparseExtra/MarketIO.h @@ -104,7 +104,7 @@ namespace internal out << value.real << " " << value.imag()<< "\n"; } -} // end namepsace internal +} // end namespace internal inline bool getMarketHeader(const std::string& filename, int& sym, bool& iscomplex, bool& isvector) { diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h index b7a9d035b..30cdf4751 100644 --- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h @@ -33,6 +33,48 @@ igamma(const Eigen::ArrayBase<Derived>& a, const Eigen::ArrayBase<ExponentDerive ); } +/** \cpp11 \returns an expression of the coefficient-wise igamma_der_a(\a a, \a x) to the given arrays. + * + * This function computes the coefficient-wise derivative of the incomplete + * gamma function with respect to the parameter a. + * + * \note This function supports only float and double scalar types in c++11 + * mode. To support other scalar types, + * or float/double in non c++11 mode, the user has to provide implementations + * of igamma_der_a(T,T) for any scalar + * type T to be supported. + * + * \sa Eigen::igamma(), Eigen::lgamma() + */ +template <typename Derived, typename ExponentDerived> +inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_igamma_der_a_op<typename Derived::Scalar>, const Derived, const ExponentDerived> +igamma_der_a(const Eigen::ArrayBase<Derived>& a, const Eigen::ArrayBase<ExponentDerived>& x) { + return Eigen::CwiseBinaryOp<Eigen::internal::scalar_igamma_der_a_op<typename Derived::Scalar>, const Derived, const ExponentDerived>( + a.derived(), + x.derived()); +} + +/** \cpp11 \returns an expression of the coefficient-wise gamma_sample_der_alpha(\a alpha, \a sample) to the given arrays. + * + * This function computes the coefficient-wise derivative of the sample + * of a Gamma(alpha, 1) random variable with respect to the parameter alpha. + * + * \note This function supports only float and double scalar types in c++11 + * mode. To support other scalar types, + * or float/double in non c++11 mode, the user has to provide implementations + * of gamma_sample_der_alpha(T,T) for any scalar + * type T to be supported. + * + * \sa Eigen::igamma(), Eigen::lgamma() + */ +template <typename AlphaDerived, typename SampleDerived> +inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_gamma_sample_der_alpha_op<typename AlphaDerived::Scalar>, const AlphaDerived, const SampleDerived> +gamma_sample_der_alpha(const Eigen::ArrayBase<AlphaDerived>& alpha, const Eigen::ArrayBase<SampleDerived>& sample) { + return Eigen::CwiseBinaryOp<Eigen::internal::scalar_gamma_sample_der_alpha_op<typename AlphaDerived::Scalar>, const AlphaDerived, const SampleDerived>( + alpha.derived(), + sample.derived()); +} + /** \cpp11 \returns an expression of the coefficient-wise igammac(\a a, \a x) to the given arrays. * * This function computes the coefficient-wise complementary incomplete gamma function. diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h index 8420f0174..3a63dcdd6 100644 --- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h @@ -41,6 +41,60 @@ struct functor_traits<scalar_igamma_op<Scalar> > { }; }; +/** \internal + * \brief Template functor to compute the derivative of the incomplete gamma + * function igamma_der_a(a, x) + * + * \sa class CwiseBinaryOp, Cwise::igamma_der_a + */ +template <typename Scalar> +struct scalar_igamma_der_a_op { + EIGEN_EMPTY_STRUCT_CTOR(scalar_igamma_der_a_op) + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& a, const Scalar& x) const { + using numext::igamma_der_a; + return igamma_der_a(a, x); + } + template <typename Packet> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const { + return internal::pigamma_der_a(a, x); + } +}; +template <typename Scalar> +struct functor_traits<scalar_igamma_der_a_op<Scalar> > { + enum { + // 2x the cost of igamma + Cost = 40 * NumTraits<Scalar>::MulCost + 20 * NumTraits<Scalar>::AddCost, + PacketAccess = packet_traits<Scalar>::HasIGammaDerA + }; +}; + +/** \internal + * \brief Template functor to compute the derivative of the sample + * of a Gamma(alpha, 1) random variable with respect to the parameter alpha + * gamma_sample_der_alpha(alpha, sample) + * + * \sa class CwiseBinaryOp, Cwise::gamma_sample_der_alpha + */ +template <typename Scalar> +struct scalar_gamma_sample_der_alpha_op { + EIGEN_EMPTY_STRUCT_CTOR(scalar_gamma_sample_der_alpha_op) + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& alpha, const Scalar& sample) const { + using numext::gamma_sample_der_alpha; + return gamma_sample_der_alpha(alpha, sample); + } + template <typename Packet> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& alpha, const Packet& sample) const { + return internal::pgamma_sample_der_alpha(alpha, sample); + } +}; +template <typename Scalar> +struct functor_traits<scalar_gamma_sample_der_alpha_op<Scalar> > { + enum { + // 2x the cost of igamma, minus the lgamma cost (the lgamma cancels out) + Cost = 30 * NumTraits<Scalar>::MulCost + 15 * NumTraits<Scalar>::AddCost, + PacketAccess = packet_traits<Scalar>::HasGammaSampleDerAlpha + }; +}; /** \internal * \brief Template functor to compute the complementary incomplete gamma function igammac(a, x) diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h index c5867002e..fbdfd299e 100644 --- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h @@ -33,6 +33,14 @@ template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half erfc(const Eigen::h template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igamma(const Eigen::half& a, const Eigen::half& x) { return Eigen::half(Eigen::numext::igamma(static_cast<float>(a), static_cast<float>(x))); } +template <> +EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igamma_der_a(const Eigen::half& a, const Eigen::half& x) { + return Eigen::half(Eigen::numext::igamma_der_a(static_cast<float>(a), static_cast<float>(x))); +} +template <> +EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half gamma_sample_der_alpha(const Eigen::half& alpha, const Eigen::half& sample) { + return Eigen::half(Eigen::numext::gamma_sample_der_alpha(static_cast<float>(alpha), static_cast<float>(sample))); +} template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igammac(const Eigen::half& a, const Eigen::half& x) { return Eigen::half(Eigen::numext::igammac(static_cast<float>(a), static_cast<float>(x))); } diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h index 5ab67dc1f..9a719e3c0 100644 --- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h @@ -521,6 +521,197 @@ struct cephes_helper<double> { } }; +enum IgammaComputationMode { VALUE, DERIVATIVE, SAMPLE_DERIVATIVE }; + +template <typename Scalar, IgammaComputationMode mode> +EIGEN_DEVICE_FUNC +int igamma_num_iterations() { + /* Returns the maximum number of internal iterations for igamma computation. + */ + if (mode == VALUE) { + return 2000; + } + + if (internal::is_same<Scalar, float>::value) { + return 200; + } else if (internal::is_same<Scalar, double>::value) { + return 500; + } else { + return 2000; + } +} + +template <typename Scalar, IgammaComputationMode mode> +struct igammac_cf_impl { + /* Computes igamc(a, x) or derivative (depending on the mode) + * using the continued fraction expansion of the complementary + * incomplete Gamma function. + * + * Preconditions: + * a > 0 + * x >= 1 + * x >= a + */ + EIGEN_DEVICE_FUNC + static Scalar run(Scalar a, Scalar x) { + const Scalar zero = 0; + const Scalar one = 1; + const Scalar two = 2; + const Scalar machep = cephes_helper<Scalar>::machep(); + const Scalar big = cephes_helper<Scalar>::big(); + const Scalar biginv = cephes_helper<Scalar>::biginv(); + + if ((numext::isinf)(x)) { + return zero; + } + + // continued fraction + Scalar y = one - a; + Scalar z = x + y + one; + Scalar c = zero; + Scalar pkm2 = one; + Scalar qkm2 = x; + Scalar pkm1 = x + one; + Scalar qkm1 = z * x; + Scalar ans = pkm1 / qkm1; + + Scalar dpkm2_da = zero; + Scalar dqkm2_da = zero; + Scalar dpkm1_da = zero; + Scalar dqkm1_da = -x; + Scalar dans_da = (dpkm1_da - ans * dqkm1_da) / qkm1; + + for (int i = 0; i < igamma_num_iterations<Scalar, mode>(); i++) { + c += one; + y += one; + z += two; + + Scalar yc = y * c; + Scalar pk = pkm1 * z - pkm2 * yc; + Scalar qk = qkm1 * z - qkm2 * yc; + + Scalar dpk_da = dpkm1_da * z - pkm1 - dpkm2_da * yc + pkm2 * c; + Scalar dqk_da = dqkm1_da * z - qkm1 - dqkm2_da * yc + qkm2 * c; + + if (qk != zero) { + Scalar ans_prev = ans; + ans = pk / qk; + + Scalar dans_da_prev = dans_da; + dans_da = (dpk_da - ans * dqk_da) / qk; + + if (mode == VALUE) { + if (numext::abs(ans_prev - ans) <= machep * numext::abs(ans)) { + break; + } + } else { + if (numext::abs(dans_da - dans_da_prev) <= machep) { + break; + } + } + } + + pkm2 = pkm1; + pkm1 = pk; + qkm2 = qkm1; + qkm1 = qk; + + dpkm2_da = dpkm1_da; + dpkm1_da = dpk_da; + dqkm2_da = dqkm1_da; + dqkm1_da = dqk_da; + + if (numext::abs(pk) > big) { + pkm2 *= biginv; + pkm1 *= biginv; + qkm2 *= biginv; + qkm1 *= biginv; + + dpkm2_da *= biginv; + dpkm1_da *= biginv; + dqkm2_da *= biginv; + dqkm1_da *= biginv; + } + } + + /* Compute x**a * exp(-x) / gamma(a) */ + Scalar logax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a); + Scalar dlogax_da = numext::log(x) - digamma_impl<Scalar>::run(a); + Scalar ax = numext::exp(logax); + Scalar dax_da = ax * dlogax_da; + + switch (mode) { + case VALUE: + return ans * ax; + case DERIVATIVE: + return ans * dax_da + dans_da * ax; + case SAMPLE_DERIVATIVE: + return -(dans_da + ans * dlogax_da) * x; + } + } +}; + +template <typename Scalar, IgammaComputationMode mode> +struct igamma_series_impl { + /* Computes igam(a, x) or its derivative (depending on the mode) + * using the series expansion of the incomplete Gamma function. + * + * Preconditions: + * x > 0 + * a > 0 + * !(x > 1 && x > a) + */ + EIGEN_DEVICE_FUNC + static Scalar run(Scalar a, Scalar x) { + const Scalar zero = 0; + const Scalar one = 1; + const Scalar machep = cephes_helper<Scalar>::machep(); + + /* power series */ + Scalar r = a; + Scalar c = one; + Scalar ans = one; + + Scalar dc_da = zero; + Scalar dans_da = zero; + + for (int i = 0; i < igamma_num_iterations<Scalar, mode>(); i++) { + r += one; + Scalar term = x / r; + Scalar dterm_da = -x / (r * r); + dc_da = term * dc_da + dterm_da * c; + dans_da += dc_da; + c *= term; + ans += c; + + if (mode == VALUE) { + if (c <= machep * ans) { + break; + } + } else { + if (numext::abs(dc_da) <= machep * numext::abs(dans_da)) { + break; + } + } + } + + /* Compute x**a * exp(-x) / gamma(a + 1) */ + Scalar logax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a + one); + Scalar dlogax_da = numext::log(x) - digamma_impl<Scalar>::run(a + one); + Scalar ax = numext::exp(logax); + Scalar dax_da = ax * dlogax_da; + + switch (mode) { + case VALUE: + return ans * ax; + case DERIVATIVE: + return ans * dax_da + dans_da * ax; + case SAMPLE_DERIVATIVE: + return -(dans_da + ans * dlogax_da) * x / a; + } + } +}; + #if !EIGEN_HAS_C99_MATH template <typename Scalar> @@ -535,8 +726,6 @@ struct igammac_impl { #else -template <typename Scalar> struct igamma_impl; // predeclare igamma_impl - template <typename Scalar> struct igammac_impl { EIGEN_DEVICE_FUNC @@ -604,97 +793,15 @@ struct igammac_impl { return nan; } - if ((numext::isnan)(a) || (numext::isnan)(x)) { // propagate nans + if ((numext::isnan)(a) || (numext::isnan)(x)) { // propagate nans return nan; } if ((x < one) || (x < a)) { - /* The checks above ensure that we meet the preconditions for - * igamma_impl::Impl(), so call it, rather than igamma_impl::Run(). - * Calling Run() would also work, but in that case the compiler may not be - * able to prove that igammac_impl::Run and igamma_impl::Run are not - * mutually recursive. This leads to worse code, particularly on - * platforms like nvptx, where recursion is allowed only begrudgingly. - */ - return (one - igamma_impl<Scalar>::Impl(a, x)); - } - - return Impl(a, x); - } - - private: - /* igamma_impl calls igammac_impl::Impl. */ - friend struct igamma_impl<Scalar>; - - /* Actually computes igamc(a, x). - * - * Preconditions: - * a > 0 - * x >= 1 - * x >= a - */ - EIGEN_DEVICE_FUNC static Scalar Impl(Scalar a, Scalar x) { - const Scalar zero = 0; - const Scalar one = 1; - const Scalar two = 2; - const Scalar machep = cephes_helper<Scalar>::machep(); - const Scalar maxlog = numext::log(NumTraits<Scalar>::highest()); - const Scalar big = cephes_helper<Scalar>::big(); - const Scalar biginv = cephes_helper<Scalar>::biginv(); - const Scalar inf = NumTraits<Scalar>::infinity(); - - Scalar ans, ax, c, yc, r, t, y, z; - Scalar pk, pkm1, pkm2, qk, qkm1, qkm2; - - if (x == inf) return zero; // std::isinf crashes on CUDA - - /* Compute x**a * exp(-x) / gamma(a) */ - ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a); - if (ax < -maxlog) { // underflow - return zero; - } - ax = numext::exp(ax); - - // continued fraction - y = one - a; - z = x + y + one; - c = zero; - pkm2 = one; - qkm2 = x; - pkm1 = x + one; - qkm1 = z * x; - ans = pkm1 / qkm1; - - for (int i = 0; i < 2000; i++) { - c += one; - y += one; - z += two; - yc = y * c; - pk = pkm1 * z - pkm2 * yc; - qk = qkm1 * z - qkm2 * yc; - if (qk != zero) { - r = pk / qk; - t = numext::abs((ans - r) / r); - ans = r; - } else { - t = one; - } - pkm2 = pkm1; - pkm1 = pk; - qkm2 = qkm1; - qkm1 = qk; - if (numext::abs(pk) > big) { - pkm2 *= biginv; - pkm1 *= biginv; - qkm2 *= biginv; - qkm1 *= biginv; - } - if (t <= machep) { - break; - } + return (one - igamma_series_impl<Scalar, VALUE>::run(a, x)); } - return (ans * ax); + return igammac_cf_impl<Scalar, VALUE>::run(a, x); } }; @@ -704,15 +811,10 @@ struct igammac_impl { * Implementation of igamma (incomplete gamma integral), based on Cephes but requires C++11/C99 * ************************************************************************************************/ -template <typename Scalar> -struct igamma_retval { - typedef Scalar type; -}; - #if !EIGEN_HAS_C99_MATH -template <typename Scalar> -struct igamma_impl { +template <typename Scalar, IgammaComputationMode mode> +struct igamma_generic_impl { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar x) { EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), @@ -723,69 +825,17 @@ struct igamma_impl { #else -template <typename Scalar> -struct igamma_impl { +template <typename Scalar, IgammaComputationMode mode> +struct igamma_generic_impl { EIGEN_DEVICE_FUNC static Scalar run(Scalar a, Scalar x) { - /* igam() - * Incomplete gamma integral - * - * - * - * SYNOPSIS: - * - * double a, x, y, igam(); - * - * y = igam( a, x ); - * - * DESCRIPTION: - * - * The function is defined by - * - * x - * - - * 1 | | -t a-1 - * igam(a,x) = ----- | e t dt. - * - | | - * | (a) - - * 0 - * - * - * In this implementation both arguments must be positive. - * The integral is evaluated by either a power series or - * continued fraction expansion, depending on the relative - * values of a and x. - * - * ACCURACY (double): - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE 0,30 200000 3.6e-14 2.9e-15 - * IEEE 0,100 300000 9.9e-14 1.5e-14 - * - * - * ACCURACY (float): - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE 0,30 20000 7.8e-6 5.9e-7 - * - */ - /* - Cephes Math Library Release 2.2: June, 1992 - Copyright 1985, 1987, 1992 by Stephen L. Moshier - Direct inquiries to 30 Frost Street, Cambridge, MA 02140 - */ - - - /* left tail of incomplete gamma function: - * - * inf. k - * a -x - x - * x e > ---------- - * - - - * k=0 | (a+k+1) + /* Depending on the mode, returns + * - VALUE: incomplete Gamma function igamma(a, x) + * - DERIVATIVE: derivative of incomplete Gamma function d/da igamma(a, x) + * - SAMPLE_DERIVATIVE: implicit derivative of a Gamma random variable + * x ~ Gamma(x | a, 1), dx/da = -1 / Gamma(x | a, 1) * d igamma(a, x) / dx * + * Derivatives are implemented by forward-mode differentiation. */ const Scalar zero = 0; const Scalar one = 1; @@ -797,71 +847,167 @@ struct igamma_impl { return nan; } - if ((numext::isnan)(a) || (numext::isnan)(x)) { // propagate nans + if ((numext::isnan)(a) || (numext::isnan)(x)) { // propagate nans return nan; } if ((x > one) && (x > a)) { - /* The checks above ensure that we meet the preconditions for - * igammac_impl::Impl(), so call it, rather than igammac_impl::Run(). - * Calling Run() would also work, but in that case the compiler may not be - * able to prove that igammac_impl::Run and igamma_impl::Run are not - * mutually recursive. This leads to worse code, particularly on - * platforms like nvptx, where recursion is allowed only begrudgingly. - */ - return (one - igammac_impl<Scalar>::Impl(a, x)); + Scalar ret = igammac_cf_impl<Scalar, mode>::run(a, x); + if (mode == VALUE) { + return one - ret; + } else { + return -ret; + } } - return Impl(a, x); + return igamma_series_impl<Scalar, mode>::run(a, x); } +}; + +#endif // EIGEN_HAS_C99_MATH - private: - /* igammac_impl calls igamma_impl::Impl. */ - friend struct igammac_impl<Scalar>; +template <typename Scalar> +struct igamma_retval { + typedef Scalar type; +}; - /* Actually computes igam(a, x). +template <typename Scalar> +struct igamma_impl : igamma_generic_impl<Scalar, VALUE> { + /* igam() + * Incomplete gamma integral. + * + * The CDF of Gamma(a, 1) random variable at the point x. + * + * Accuracy estimation. For each a in [10^-2, 10^-1...10^3] we sample + * 50 Gamma random variables x ~ Gamma(x | a, 1), a total of 300 points. + * The ground truth is computed by mpmath. Mean absolute error: + * float: 1.26713e-05 + * double: 2.33606e-12 + * + * Cephes documentation below. + * + * SYNOPSIS: + * + * double a, x, y, igam(); + * + * y = igam( a, x ); + * + * DESCRIPTION: + * + * The function is defined by + * + * x + * - + * 1 | | -t a-1 + * igam(a,x) = ----- | e t dt. + * - | | + * | (a) - + * 0 + * + * + * In this implementation both arguments must be positive. + * The integral is evaluated by either a power series or + * continued fraction expansion, depending on the relative + * values of a and x. + * + * ACCURACY (double): + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,30 200000 3.6e-14 2.9e-15 + * IEEE 0,100 300000 9.9e-14 1.5e-14 + * + * + * ACCURACY (float): + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,30 20000 7.8e-6 5.9e-7 * - * Preconditions: - * x > 0 - * a > 0 - * !(x > 1 && x > a) */ - EIGEN_DEVICE_FUNC static Scalar Impl(Scalar a, Scalar x) { - const Scalar zero = 0; - const Scalar one = 1; - const Scalar machep = cephes_helper<Scalar>::machep(); - const Scalar maxlog = numext::log(NumTraits<Scalar>::highest()); + /* + Cephes Math Library Release 2.2: June, 1992 + Copyright 1985, 1987, 1992 by Stephen L. Moshier + Direct inquiries to 30 Frost Street, Cambridge, MA 02140 + */ - Scalar ans, ax, c, r; + /* left tail of incomplete gamma function: + * + * inf. k + * a -x - x + * x e > ---------- + * - - + * k=0 | (a+k+1) + * + */ +}; - /* Compute x**a * exp(-x) / gamma(a) */ - ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a); - if (ax < -maxlog) { - // underflow - return zero; - } - ax = numext::exp(ax); +template <typename Scalar> +struct igamma_der_a_retval : igamma_retval<Scalar> {}; - /* power series */ - r = a; - c = one; - ans = one; +template <typename Scalar> +struct igamma_der_a_impl : igamma_generic_impl<Scalar, DERIVATIVE> { + /* Derivative of the incomplete Gamma function with respect to a. + * + * Computes d/da igamma(a, x) by forward differentiation of the igamma code. + * + * Accuracy estimation. For each a in [10^-2, 10^-1...10^3] we sample + * 50 Gamma random variables x ~ Gamma(x | a, 1), a total of 300 points. + * The ground truth is computed by mpmath. Mean absolute error: + * float: 6.17992e-07 + * double: 4.60453e-12 + * + * Reference: + * R. Moore. "Algorithm AS 187: Derivatives of the incomplete gamma + * integral". Journal of the Royal Statistical Society. 1982 + */ +}; - for (int i = 0; i < 2000; i++) { - r += one; - c *= x/r; - ans += c; - if (c/ans <= machep) { - break; - } - } +template <typename Scalar> +struct gamma_sample_der_alpha_retval : igamma_retval<Scalar> {}; - return (ans * ax / a); - } +template <typename Scalar> +struct gamma_sample_der_alpha_impl + : igamma_generic_impl<Scalar, SAMPLE_DERIVATIVE> { + /* Derivative of a Gamma random variable sample with respect to alpha. + * + * Consider a sample of a Gamma random variable with the concentration + * parameter alpha: sample ~ Gamma(alpha, 1). The reparameterization + * derivative that we want to compute is dsample / dalpha = + * d igammainv(alpha, u) / dalpha, where u = igamma(alpha, sample). + * However, this formula is numerically unstable and expensive, so instead + * we use implicit differentiation: + * + * igamma(alpha, sample) = u, where u ~ Uniform(0, 1). + * Apply d / dalpha to both sides: + * d igamma(alpha, sample) / dalpha + * + d igamma(alpha, sample) / dsample * dsample/dalpha = 0 + * d igamma(alpha, sample) / dalpha + * + Gamma(sample | alpha, 1) dsample / dalpha = 0 + * dsample/dalpha = - (d igamma(alpha, sample) / dalpha) + * / Gamma(sample | alpha, 1) + * + * Here Gamma(sample | alpha, 1) is the PDF of the Gamma distribution + * (note that the derivative of the CDF w.r.t. sample is the PDF). + * See the reference below for more details. + * + * The derivative of igamma(alpha, sample) is computed by forward + * differentiation of the igamma code. Division by the Gamma PDF is performed + * in the same code, increasing the accuracy and speed due to cancellation + * of some terms. + * + * Accuracy estimation. For each alpha in [10^-2, 10^-1...10^3] we sample + * 50 Gamma random variables sample ~ Gamma(sample | alpha, 1), a total of 300 + * points. The ground truth is computed by mpmath. Mean absolute error: + * float: 2.1686e-06 + * double: 1.4774e-12 + * + * Reference: + * M. Figurnov, S. Mohamed, A. Mnih "Implicit Reparameterization Gradients". + * 2018 + */ }; -#endif // EIGEN_HAS_C99_MATH - /***************************************************************************** * Implementation of Riemann zeta function of two arguments, based on Cephes * *****************************************************************************/ @@ -1574,6 +1720,8 @@ struct betainc_impl<double> { } }; +#endif // EIGEN_HAS_C99_MATH + /**************************************************************************** * Implementation of Bessel function, based on Cephes * ****************************************************************************/ @@ -1902,8 +2050,6 @@ struct i1e_impl<double> { } }; -#endif // EIGEN_HAS_C99_MATH - } // end namespace internal namespace numext { @@ -1951,6 +2097,18 @@ EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igamma, Scalar) } template <typename Scalar> +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igamma_der_a, Scalar) + igamma_der_a(const Scalar& a, const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(igamma_der_a, Scalar)::run(a, x); +} + +template <typename Scalar> +EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(gamma_sample_der_alpha, Scalar) + gamma_sample_der_alpha(const Scalar& a, const Scalar& x) { + return EIGEN_MATHFUNC_IMPL(gamma_sample_der_alpha, Scalar)::run(a, x); +} + +template <typename Scalar> EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igammac, Scalar) igammac(const Scalar& a, const Scalar& x) { return EIGEN_MATHFUNC_IMPL(igammac, Scalar)::run(a, x); diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h index 4c176716b..465f41d54 100644 --- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h @@ -42,6 +42,21 @@ Packet perfc(const Packet& a) { using numext::erfc; return erfc(a); } template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pigamma(const Packet& a, const Packet& x) { using numext::igamma; return igamma(a, x); } +/** \internal \returns the derivative of the incomplete gamma function + * igamma_der_a(\a a, \a x) */ +template <typename Packet> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pigamma_der_a(const Packet& a, const Packet& x) { + using numext::igamma_der_a; return igamma_der_a(a, x); +} + +/** \internal \returns compute the derivative of the sample + * of Gamma(alpha, 1) random variable with respect to the parameter a + * gamma_sample_der_alpha(\a alpha, \a sample) */ +template <typename Packet> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pgamma_sample_der_alpha(const Packet& alpha, const Packet& sample) { + using numext::gamma_sample_der_alpha; return gamma_sample_der_alpha(alpha, sample); +} + /** \internal \returns the complementary incomplete gamma function igammac(\a a, \a x) */ template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pigammac(const Packet& a, const Packet& x) { using numext::igammac; return igammac(a, x); } diff --git a/unsupported/Eigen/src/SpecialFunctions/arch/CUDA/CudaSpecialFunctions.h b/unsupported/Eigen/src/SpecialFunctions/arch/CUDA/CudaSpecialFunctions.h index c25fea0b3..020ac1b62 100644 --- a/unsupported/Eigen/src/SpecialFunctions/arch/CUDA/CudaSpecialFunctions.h +++ b/unsupported/Eigen/src/SpecialFunctions/arch/CUDA/CudaSpecialFunctions.h @@ -120,6 +120,41 @@ double2 pigamma<double2>(const double2& a, const double2& x) return make_double2(igamma(a.x, x.x), igamma(a.y, x.y)); } +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pigamma_der_a<float4>( + const float4& a, const float4& x) { + using numext::igamma_der_a; + return make_float4(igamma_der_a(a.x, x.x), igamma_der_a(a.y, x.y), + igamma_der_a(a.z, x.z), igamma_der_a(a.w, x.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pigamma_der_a<double2>(const double2& a, const double2& x) { + using numext::igamma_der_a; + return make_double2(igamma_der_a(a.x, x.x), igamma_der_a(a.y, x.y)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pgamma_sample_der_alpha<float4>( + const float4& alpha, const float4& sample) { + using numext::gamma_sample_der_alpha; + return make_float4( + gamma_sample_der_alpha(alpha.x, sample.x), + gamma_sample_der_alpha(alpha.y, sample.y), + gamma_sample_der_alpha(alpha.z, sample.z), + gamma_sample_der_alpha(alpha.w, sample.w)); +} + +template <> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2 +pgamma_sample_der_alpha<double2>(const double2& alpha, const double2& sample) { + using numext::gamma_sample_der_alpha; + return make_double2( + gamma_sample_der_alpha(alpha.x, sample.x), + gamma_sample_der_alpha(alpha.y, sample.y)); +} + template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pigammac<float4>(const float4& a, const float4& x) { diff --git a/unsupported/Eigen/src/Splines/SplineFitting.h b/unsupported/Eigen/src/Splines/SplineFitting.h index c761a9b3d..1a4b80a2f 100644 --- a/unsupported/Eigen/src/Splines/SplineFitting.h +++ b/unsupported/Eigen/src/Splines/SplineFitting.h @@ -181,7 +181,7 @@ namespace Eigen * \ingroup Splines_Module * * \param[in] pts The data points to which a spline should be fit. - * \param[out] chord_lengths The resulting chord lenggth vector. + * \param[out] chord_lengths The resulting chord length vector. * * \sa Les Piegl and Wayne Tiller, The NURBS book (2nd ed.), 1997, 9.2.1 Global Curve Interpolation to Point Data **/ |