// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr> // // 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/. #include <limits.h>#include "main.h" #include "../Eigen/SpecialFunctions" // Hack to allow "implicit" conversions from double to Scalar via comma-initialization. template <typename Derived>operator <<(Eigen::DenseBase<Derived>& dense, double v) {return (dense << static_cast <typename Derived::Scalar>(v));template <typename XprType>operator ,(Eigen::CommaInitializer<XprType>& ci, double v) {return (ci, static_cast <typename XprType::Scalar>(v));template <typename X, typename Y>void verify_component_wise(const X& x, const Y& y)for (Index i=0; i<x.size(); ++i)if ((numext::isfinite)(y(i)))else if ((numext::isnan)(y(i)))else template <typename ArrayType> void array_special_functions()using std::abs;using std::sqrt;typedef typename ArrayType::Scalar Scalar;typedef typename NumTraits<Scalar>::Real RealScalar;// API #if EIGEN_HAS_C99_MATH#endif // EIGEN_HAS_C99_MATH #if EIGEN_HAS_C99_MATH// check special functions (comparing against numpy implementation) if (!NumTraits<Scalar>::IsComplex)// Test various propreties of igamma & igammac. These are normalized // gamma integrals where // igammac(a, x) = Gamma(a, x) / Gamma(a) // igamma(a, x) = gamma(a, x) / Gamma(a) // where Gamma and gamma are considered the standard unnormalized // upper and lower incomplete gamma functions, respectively. // Gamma(a, 0) == Gamma(a) // Gamma(a, x) + gamma(a, x) == Gamma(a) // Gamma(a, x) == (a - 1) * Gamma(a-1, x) + x^(a-1) * exp(-x) p());// gamma(a, x) == (a - 1) * gamma(a-1, x) - x^(a-1) * exp(-x) p());// Verify for large a and x that values are between 0 and 1. int max_exponent = std::numeric_limits<Scalar>::max_exponent10;for (int i = 0; i < a.size(); ++i) {// Check exact values of igamma and igammac against a third party calculation. 0.5)};0.5)};// location i*6+j corresponds to a_s[i], x_s[j]. for (int i = 0; i < 6; ++i) {for (int j = 0; j < 6; ++j) {if ((std::isnan)(igamma_s[i][j])) {else {if ((std::isnan)(igammac_s[i][j])) {else {#endif // EIGEN_HAS_C99_MATH // Check the ndtri function against scipy.special.ndtri 004, 2.3263478740408408, -2.3263478740408408;// ndtri(normal_cdf(x)) ~= x using std::sqrt;// Check the zeta function against scipy.special.zeta 4097, plusinf, nan, plusinf, nan, plusinf;// digamma 10340372392849, nan, nan, nan, nan;#if EIGEN_HAS_C99_MATH6, -2.47810300902e-07, -8.29668781082e-09, -0.434562276666, 0.567742190178, -0.0108615497927, nan, nan, plusinf, nan, plusinf;if (sizeof (RealScalar)>=8) { // double // Reason for commented line: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1232 // CALL_SUBTEST( res = x.polygamma(n); verify_component_wise(res, ref); ); else {// CALL_SUBTEST( res = x.polygamma(n); verify_component_wise(res.head(8), ref.head(8)); ); ;#endif #if EIGEN_HAS_C99_MATH// Inputs and ground truth generated with scipy via: // a = np.logspace(-3, 3, 5) - 1e-3 // b = np.logspace(-3, 3, 5) - 1e-3 // x = np.linspace(-0.1, 1.1, 5) // (full_a, full_b, full_x) = np.vectorize(lambda a, b, x: (a, b, x))(*np.ix_(a, b, x)) // full_a = full_a.flatten().tolist() # same for full_b, full_x // v = scipy.special.betainc(full_a, full_b, full_x).flatten().tolist() // // Note in Eigen, we call betainc with arguments in the order (x, a, b). // Test various properties of betainc const Scalar eps = std::numeric_limits<Scalar>::epsilon();// betainc(a, 1, x) == x**a // betainc(1, b, x) == 1 - (1 - x)**b // betainc(a, b, x) == 1 - betainc(b, a, 1-x) // betainc(a+1, b, x) = betainc(a, b, x) - x**a * (1 - x)**b / (a * beta(a, b)) // Add eps to rhs and lhs so that component-wise test doesn't result in // nans when both outputs are zeros. if (sizeof (Scalar) >= 8) { // double else {// Reason for limited test: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1232 // betainc(a, b+1, x) = betainc(a, b, x) + x**a * (1 - x)**b / (b * beta(a, b)) // Add eps to rhs and lhs so that component-wise test doesn't result in // nans when both outputs are zeros. #endif // EIGEN_HAS_C99_MATH /* Code to generate the data for the following two test cases. N = 5 np.random.seed(3) a = np.logspace(-2, 3, 6) a = np.ravel(np.tile(np.reshape(a, [-1, 1]), [1, N])) x = np.random.gamma(a, 1.0) x = np.maximum(x, np.finfo(np.float32).tiny) def igamma(a, x): return mpmath.gammainc(a, 0, x, regularized=True) def igamma_der_a(a, x): res = mpmath.diff(lambda a_prime: igamma(a_prime, x), a) return np.float64(res) def gamma_sample_der_alpha(a, x): igamma_x = igamma(a, x) def igammainv_of_igamma(a_prime): return mpmath.findroot(lambda x_prime: igamma(a_prime, x_prime) - igamma_x, x, solver='newton') return np.float64(mpmath.diff(igammainv_of_igamma, a)) v_igamma_der_a = np.vectorize(igamma_der_a)(a, x) v_gamma_sample_der_alpha = np.vectorize(gamma_sample_der_alpha)(a, x) #if EIGEN_HAS_C99_MATH// Test igamma_der_a // Test gamma_sample_der_alpha #endif // EIGEN_HAS_C99_MATH // TODO(cantonios): half/bfloat16 don't have enough precision to reproduce results above. // CALL_SUBTEST_3(array_special_functions<ArrayX<Eigen::half>>()); // CALL_SUBTEST_4(array_special_functions<ArrayX<Eigen::bfloat16>>()); quality 85%
¤ Dauer der Verarbeitung: 0.14 Sekunden
(vorverarbeitet)
¤ *© Formatika GbR, Deutschland
