Functions
special_functions.cpp File Reference
#include <limits.h>
#include "main.h"
#include "../Eigen/SpecialFunctions"

Functions

template<typename Derived >
Eigen::CommaInitializer< Derived > operator<< (Eigen::DenseBase< Derived > &dense, double v)
 
template<typename XprType >
Eigen::CommaInitializer< XprType > & operator, (Eigen::CommaInitializer< XprType > &ci, double v)
 
template<typename X , typename Y >
void verify_component_wise (const X &x, const Y &y)
 
template<typename ArrayType >
void array_special_functions ()
 
 EIGEN_DECLARE_TEST (special_functions)
 

Function Documentation

◆ array_special_functions()

template<typename ArrayType >
void array_special_functions ( )
38  {
39  using std::abs;
40  using std::sqrt;
41  typedef typename ArrayType::Scalar Scalar;
42  typedef typename NumTraits<Scalar>::Real RealScalar;
43 
44  Scalar plusinf = std::numeric_limits<Scalar>::infinity();
45  Scalar nan = std::numeric_limits<Scalar>::quiet_NaN();
46 
47  Index rows = internal::random<Index>(1, 30);
48  Index cols = 1;
49 
50  // API
51  {
52  ArrayType m1 = ArrayType::Random(rows, cols);
53 #if EIGEN_HAS_C99_MATH
54  VERIFY_IS_APPROX(m1.lgamma(), lgamma(m1));
55  VERIFY_IS_APPROX(m1.digamma(), digamma(m1));
56  VERIFY_IS_APPROX(m1.erf(), erf(m1));
57  VERIFY_IS_APPROX(m1.erfc(), erfc(m1));
58 #endif // EIGEN_HAS_C99_MATH
59  }
60 
61 #if EIGEN_HAS_C99_MATH
62  // check special functions (comparing against numpy implementation)
63  if (!NumTraits<Scalar>::IsComplex) {
64  {
65  ArrayType m1 = ArrayType::Random(rows, cols);
66  ArrayType m2 = ArrayType::Random(rows, cols);
67 
68  // Test various propreties of igamma & igammac. These are normalized
69  // gamma integrals where
70  // igammac(a, x) = Gamma(a, x) / Gamma(a)
71  // igamma(a, x) = gamma(a, x) / Gamma(a)
72  // where Gamma and gamma are considered the standard unnormalized
73  // upper and lower incomplete gamma functions, respectively.
74  ArrayType a = m1.abs() + Scalar(2);
75  ArrayType x = m2.abs() + Scalar(2);
76  ArrayType zero = ArrayType::Zero(rows, cols);
77  ArrayType one = ArrayType::Constant(rows, cols, Scalar(1.0));
78  ArrayType a_m1 = a - one;
79  ArrayType Gamma_a_x = Eigen::igammac(a, x) * a.lgamma().exp();
80  ArrayType Gamma_a_m1_x = Eigen::igammac(a_m1, x) * a_m1.lgamma().exp();
81  ArrayType gamma_a_x = Eigen::igamma(a, x) * a.lgamma().exp();
82  ArrayType gamma_a_m1_x = Eigen::igamma(a_m1, x) * a_m1.lgamma().exp();
83 
84  // Gamma(a, 0) == Gamma(a)
85  VERIFY_IS_APPROX(Eigen::igammac(a, zero), one);
86 
87  // Gamma(a, x) + gamma(a, x) == Gamma(a)
88  VERIFY_IS_APPROX(Gamma_a_x + gamma_a_x, a.lgamma().exp());
89 
90  // Gamma(a, x) == (a - 1) * Gamma(a-1, x) + x^(a-1) * exp(-x)
91  VERIFY_IS_APPROX(Gamma_a_x, (a - Scalar(1)) * Gamma_a_m1_x + x.pow(a - Scalar(1)) * (-x).exp());
92 
93  // gamma(a, x) == (a - 1) * gamma(a-1, x) - x^(a-1) * exp(-x)
94  VERIFY_IS_APPROX(gamma_a_x, (a - Scalar(1)) * gamma_a_m1_x - x.pow(a - Scalar(1)) * (-x).exp());
95  }
96  {
97  // Verify for large a and x that values are between 0 and 1.
98  ArrayType m1 = ArrayType::Random(rows, cols);
99  ArrayType m2 = ArrayType::Random(rows, cols);
100  int max_exponent = std::numeric_limits<Scalar>::max_exponent10;
101  ArrayType a = m1.abs() * Scalar(pow(10., max_exponent - 1));
102  ArrayType x = m2.abs() * Scalar(pow(10., max_exponent - 1));
103  for (int i = 0; i < a.size(); ++i) {
104  Scalar igam = numext::igamma(a(i), x(i));
105  VERIFY(0 <= igam);
106  VERIFY(igam <= 1);
107  }
108  }
109 
110  {
111  // Check exact values of igamma and igammac against a third party calculation.
112  Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)};
113  Scalar x_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)};
114 
115  // location i*6+j corresponds to a_s[i], x_s[j].
116  Scalar igamma_s[][6] = {
117  {Scalar(0.0), nan, nan, nan, nan, nan},
118  {Scalar(0.0), Scalar(0.6321205588285578), Scalar(0.7768698398515702), Scalar(0.9816843611112658),
119  Scalar(9.999500016666262e-05), Scalar(1.0)},
120  {Scalar(0.0), Scalar(0.4275932955291202), Scalar(0.608374823728911), Scalar(0.9539882943107686),
121  Scalar(7.522076445089201e-07), Scalar(1.0)},
122  {Scalar(0.0), Scalar(0.01898815687615381), Scalar(0.06564245437845008), Scalar(0.5665298796332909),
123  Scalar(4.166333347221828e-18), Scalar(1.0)},
124  {Scalar(0.0), Scalar(0.9999780593618628), Scalar(0.9999899967080838), Scalar(0.9999996219837988),
125  Scalar(0.9991370418689945), Scalar(1.0)},
126  {Scalar(0.0), Scalar(0.0), Scalar(0.0), Scalar(0.0), Scalar(0.0), Scalar(0.5042041932513908)}};
127  Scalar igammac_s[][6] = {
128  {nan, nan, nan, nan, nan, nan},
129  {Scalar(1.0), Scalar(0.36787944117144233), Scalar(0.22313016014842982), Scalar(0.018315638888734182),
130  Scalar(0.9999000049998333), Scalar(0.0)},
131  {Scalar(1.0), Scalar(0.5724067044708798), Scalar(0.3916251762710878), Scalar(0.04601170568923136),
132  Scalar(0.9999992477923555), Scalar(0.0)},
133  {Scalar(1.0), Scalar(0.9810118431238462), Scalar(0.9343575456215499), Scalar(0.4334701203667089), Scalar(1.0),
134  Scalar(0.0)},
135  {Scalar(1.0), Scalar(2.1940638138146658e-05), Scalar(1.0003291916285e-05), Scalar(3.7801620118431334e-07),
136  Scalar(0.0008629581310054535), Scalar(0.0)},
137  {Scalar(1.0), Scalar(1.0), Scalar(1.0), Scalar(1.0), Scalar(1.0), Scalar(0.49579580674813944)}};
138 
139  for (int i = 0; i < 6; ++i) {
140  for (int j = 0; j < 6; ++j) {
141  if ((std::isnan)(igamma_s[i][j])) {
142  VERIFY((std::isnan)(numext::igamma(a_s[i], x_s[j])));
143  } else {
144  VERIFY_IS_APPROX(numext::igamma(a_s[i], x_s[j]), igamma_s[i][j]);
145  }
146 
147  if ((std::isnan)(igammac_s[i][j])) {
148  VERIFY((std::isnan)(numext::igammac(a_s[i], x_s[j])));
149  } else {
150  VERIFY_IS_APPROX(numext::igammac(a_s[i], x_s[j]), igammac_s[i][j]);
151  }
152  }
153  }
154  }
155  }
156 #endif // EIGEN_HAS_C99_MATH
157 
158  // Check the ndtri function against scipy.special.ndtri
159  {
160  ArrayType x(11), res(11), ref(11);
161  x << 0.5, 0.2, 0.8, 0.9, 0.1, 0.99, 0.01, 0, 1, -0.01, 1.01;
162  ref << 0., -0.8416212335729142, 0.8416212335729142, 1.2815515655446004, -1.2815515655446004, 2.3263478740408408,
163  -2.3263478740408408, -plusinf, plusinf, nan, nan;
164  CALL_SUBTEST(verify_component_wise(ref, ref););
165  CALL_SUBTEST(res = x.ndtri(); verify_component_wise(res, ref););
166  CALL_SUBTEST(res = ndtri(x); verify_component_wise(res, ref););
167 
168  // ndtri(normal_cdf(x)) ~= x
169  CALL_SUBTEST(ArrayType m1 = ArrayType::Random(32); using std::sqrt;
170 
171  ArrayType cdf_val = (m1 / Scalar(sqrt(2.))).erf(); cdf_val = (cdf_val + Scalar(1)) / Scalar(2);
172  verify_component_wise(cdf_val.ndtri(), m1););
173  }
174 
175  // Check the zeta function against scipy.special.zeta
176  {
177  ArrayType x(11), q(11), res(11), ref(11);
178  x << 1.5, 4, 10.5, 10000.5, 3, 1, 0.9, 2, 3, 4, 2000;
179  q << 2, 1.5, 3, 1.0001, -2.5, 1.2345, 1.2345, -1, -2, -3, 2000;
180  ref << 1.61237534869, 0.234848505667, 1.03086757337e-5, 0.367879440865, 0.054102025820864097, plusinf, nan, plusinf,
181  nan, plusinf, 0;
182  CALL_SUBTEST(verify_component_wise(ref, ref););
183  CALL_SUBTEST(res = x.zeta(q); verify_component_wise(res, ref););
184  CALL_SUBTEST(res = zeta(x, q); verify_component_wise(res, ref););
185  }
186 
187  // digamma
188  {
189  ArrayType x(9), res(9), ref(9);
190  x << 1, 1.5, 4, -10.5, 10000.5, 0, -1, -2, -3;
191  ref << -0.5772156649015329, 0.03648997397857645, 1.2561176684318, 2.398239129535781, 9.210340372392849, nan, nan,
192  nan, nan;
193  CALL_SUBTEST(verify_component_wise(ref, ref););
194 
195  CALL_SUBTEST(res = x.digamma(); verify_component_wise(res, ref););
196  CALL_SUBTEST(res = digamma(x); verify_component_wise(res, ref););
197  }
198 
199 #if EIGEN_HAS_C99_MATH
200  {
201  ArrayType n(16), x(16), res(16), ref(16);
202  n << 1, 1, 1, 1.5, 17, 31, 28, 8, 42, 147, 170, -1, 0, 1, 2, 3;
203  x << 2, 3, 25.5, 1.5, 4.7, 11.8, 17.7, 30.2, 15.8, 54.1, 64, -1, -2, -3, -4, -5;
204  ref << 0.644934066848, 0.394934066848, 0.0399946696496, nan, 293.334565435, 0.445487887616, -2.47810300902e-07,
205  -8.29668781082e-09, -0.434562276666, 0.567742190178, -0.0108615497927, nan, nan, plusinf, nan, plusinf;
206  CALL_SUBTEST(verify_component_wise(ref, ref););
207 
208  if (sizeof(RealScalar) >= 8) { // double
209  // Reason for commented line: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1232
210  // CALL_SUBTEST( res = x.polygamma(n); verify_component_wise(res, ref); );
211  CALL_SUBTEST(res = polygamma(n, x); verify_component_wise(res, ref););
212  } else {
213  // CALL_SUBTEST( res = x.polygamma(n); verify_component_wise(res.head(8), ref.head(8)); );
214  CALL_SUBTEST(res = polygamma(n, x); verify_component_wise(res.head(8), ref.head(8)););
215  }
216  }
217 #endif
218 
219 #if EIGEN_HAS_C99_MATH
220  {
221  // Inputs and ground truth generated with scipy via:
222  // a = np.logspace(-3, 3, 5) - 1e-3
223  // b = np.logspace(-3, 3, 5) - 1e-3
224  // x = np.linspace(-0.1, 1.1, 5)
225  // (full_a, full_b, full_x) = np.vectorize(lambda a, b, x: (a, b, x))(*np.ix_(a, b, x))
226  // full_a = full_a.flatten().tolist() # same for full_b, full_x
227  // v = scipy.special.betainc(full_a, full_b, full_x).flatten().tolist()
228  //
229  // Note in Eigen, we call betainc with arguments in the order (x, a, b).
230  ArrayType a(125);
231  ArrayType b(125);
232  ArrayType x(125);
233  ArrayType v(125);
234  ArrayType res(125);
235 
236  a << 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
237  0.0, 0.0, 0.0, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
238  0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
239  0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
240  0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
241  0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
242  0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999,
243  0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 31.62177660168379,
244  31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379,
245  31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379,
246  31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379,
247  31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379,
248  31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, 999.999,
249  999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999,
250  999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999;
251 
252  b << 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
253  0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, 31.62177660168379, 31.62177660168379, 31.62177660168379,
254  31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999, 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0,
255  0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.999,
256  0.999, 0.999, 0.999, 0.999, 31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379,
257  31.62177660168379, 999.999, 999.999, 999.999, 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379,
258  0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.999, 0.999, 0.999, 0.999,
259  0.999, 31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379, 999.999,
260  999.999, 999.999, 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, 0.03062277660168379,
261  0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999,
262  31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379, 999.999, 999.999,
263  999.999, 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, 0.03062277660168379,
264  0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999,
265  31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379, 31.62177660168379, 999.999, 999.999,
266  999.999, 999.999, 999.999;
267 
268  x << -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1,
269  0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2,
270  0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5,
271  0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8,
272  1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1,
273  -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1;
274 
275  v << nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
276  nan, nan, nan, nan, nan, nan, nan, nan, nan, 0.47972119876364683, 0.5, 0.5202788012363533, nan, nan,
277  0.9518683957740043, 0.9789663010413743, 0.9931729188073435, nan, nan, 0.999995949033062, 0.9999999999993698,
278  0.9999999999999999, nan, nan, 0.9999999999999999, 0.9999999999999999, 0.9999999999999999, nan, nan, nan, nan,
279  nan, nan, nan, 0.006827081192655869, 0.0210336989586256, 0.04813160422599567, nan, nan, 0.20014344256217678,
280  0.5000000000000001, 0.7998565574378232, nan, nan, 0.9991401428435834, 0.999999999698403, 0.9999999999999999,
281  nan, nan, 0.9999999999999999, 0.9999999999999999, 0.9999999999999999, nan, nan, nan, nan, nan, nan, nan,
282  1.0646600232370887e-25, 6.301722877826246e-13, 4.050966937974938e-06, nan, nan, 7.864342668429763e-23,
283  3.015969667594166e-10, 0.0008598571564165444, nan, nan, 6.031987710123844e-08, 0.5000000000000007,
284  0.9999999396801229, nan, nan, 0.9999999999999999, 0.9999999999999999, 0.9999999999999999, nan, nan, nan, nan,
285  nan, nan, nan, 0.0, 7.029920380986636e-306, 2.2450728208591345e-101, nan, nan, 0.0, 9.275871147869727e-302,
286  1.2232913026152827e-97, nan, nan, 0.0, 3.0891393081932924e-252, 2.9303043666183996e-60, nan, nan,
287  2.248913486879199e-196, 0.5000000000004947, 0.9999999999999999, nan;
288 
289  CALL_SUBTEST(res = betainc(a, b, x); verify_component_wise(res, v););
290  }
291 
292  // Test various properties of betainc
293  {
294  ArrayType m1 = ArrayType::Random(32);
295  ArrayType m2 = ArrayType::Random(32);
296  ArrayType m3 = ArrayType::Random(32);
297  ArrayType one = ArrayType::Constant(32, Scalar(1.0));
298  const Scalar eps = std::numeric_limits<Scalar>::epsilon();
299  ArrayType a = (m1 * Scalar(4)).exp();
300  ArrayType b = (m2 * Scalar(4)).exp();
301  ArrayType x = m3.abs();
302 
303  // betainc(a, 1, x) == x**a
304  CALL_SUBTEST(ArrayType test = betainc(a, one, x); ArrayType expected = x.pow(a);
305  verify_component_wise(test, expected););
306 
307  // betainc(1, b, x) == 1 - (1 - x)**b
308  CALL_SUBTEST(ArrayType test = betainc(one, b, x); ArrayType expected = one - (one - x).pow(b);
309  verify_component_wise(test, expected););
310 
311  // betainc(a, b, x) == 1 - betainc(b, a, 1-x)
312  CALL_SUBTEST(ArrayType test = betainc(a, b, x) + betainc(b, a, one - x); ArrayType expected = one;
313  verify_component_wise(test, expected););
314 
315  // betainc(a+1, b, x) = betainc(a, b, x) - x**a * (1 - x)**b / (a * beta(a, b))
316  CALL_SUBTEST(
317  ArrayType num = x.pow(a) * (one - x).pow(b);
318  ArrayType denom = a * (a.lgamma() + b.lgamma() - (a + b).lgamma()).exp();
319  // Add eps to rhs and lhs so that component-wise test doesn't result in
320  // nans when both outputs are zeros.
321  ArrayType expected = betainc(a, b, x) - num / denom + eps;
322  ArrayType test = betainc(a + one, b, x) + eps; if (sizeof(Scalar) >= 8) { // double
323  verify_component_wise(test, expected);
324  } else {
325  // Reason for limited test: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1232
326  verify_component_wise(test.head(8), expected.head(8));
327  });
328 
329  // betainc(a, b+1, x) = betainc(a, b, x) + x**a * (1 - x)**b / (b * beta(a, b))
330  CALL_SUBTEST(
331  // Add eps to rhs and lhs so that component-wise test doesn't result in
332  // nans when both outputs are zeros.
333  ArrayType num = x.pow(a) * (one - x).pow(b);
334  ArrayType denom = b * (a.lgamma() + b.lgamma() - (a + b).lgamma()).exp();
335  ArrayType expected = betainc(a, b, x) + num / denom + eps; ArrayType test = betainc(a, b + one, x) + eps;
336  verify_component_wise(test, expected););
337  }
338 #endif // EIGEN_HAS_C99_MATH
339 
340  /* Code to generate the data for the following two test cases.
341  N = 5
342  np.random.seed(3)
343 
344  a = np.logspace(-2, 3, 6)
345  a = np.ravel(np.tile(np.reshape(a, [-1, 1]), [1, N]))
346  x = np.random.gamma(a, 1.0)
347  x = np.maximum(x, np.finfo(np.float32).tiny)
348 
349  def igamma(a, x):
350  return mpmath.gammainc(a, 0, x, regularized=True)
351 
352  def igamma_der_a(a, x):
353  res = mpmath.diff(lambda a_prime: igamma(a_prime, x), a)
354  return np.float64(res)
355 
356  def gamma_sample_der_alpha(a, x):
357  igamma_x = igamma(a, x)
358  def igammainv_of_igamma(a_prime):
359  return mpmath.findroot(lambda x_prime: igamma(a_prime, x_prime) -
360  igamma_x, x, solver='newton')
361  return np.float64(mpmath.diff(igammainv_of_igamma, a))
362 
363  v_igamma_der_a = np.vectorize(igamma_der_a)(a, x)
364  v_gamma_sample_der_alpha = np.vectorize(gamma_sample_der_alpha)(a, x)
365 */
366 
367 #if EIGEN_HAS_C99_MATH
368  // Test igamma_der_a
369  {
370  ArrayType a(30);
371  ArrayType x(30);
372  ArrayType res(30);
373  ArrayType v(30);
374 
375  a << 0.01, 0.01, 0.01, 0.01, 0.01, 0.1, 0.1, 0.1, 0.1, 0.1, 1.0, 1.0, 1.0, 1.0, 1.0, 10.0, 10.0, 10.0, 10.0, 10.0,
376  100.0, 100.0, 100.0, 100.0, 100.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0;
377 
378  x << 1.25668890405e-26, 1.17549435082e-38, 1.20938905072e-05, 1.17549435082e-38, 1.17549435082e-38,
379  5.66572070696e-16, 0.0132865061065, 0.0200034203853, 6.29263709118e-17, 1.37160367764e-06, 0.333412038288,
380  1.18135687766, 0.580629033777, 0.170631439426, 0.786686768458, 7.63873279537, 13.1944344379, 11.896042354,
381  10.5830172417, 10.5020942233, 92.8918587747, 95.003720371, 86.3715926467, 96.0330217672, 82.6389930677,
382  968.702906754, 969.463546828, 1001.79726022, 955.047416547, 1044.27458568;
383 
384  v << -32.7256441441, -36.4394150514, -9.66467612263, -36.4394150514, -36.4394150514, -1.0891900302, -2.66351229645,
385  -2.48666868596, -0.929700494428, -3.56327722764, -0.455320135314, -0.391437214323, -0.491352055991,
386  -0.350454834292, -0.471773162921, -0.104084440522, -0.0723646747909, -0.0992828975532, -0.121638215446,
387  -0.122619605294, -0.0317670267286, -0.0359974812869, -0.0154359225363, -0.0375775365921, -0.00794899153653,
388  -0.00777303219211, -0.00796085782042, -0.0125850719397, -0.00455500206958, -0.00476436993148;
389 
390  CALL_SUBTEST(res = igamma_der_a(a, x); verify_component_wise(res, v););
391  }
392 
393  // Test gamma_sample_der_alpha
394  {
395  ArrayType alpha(30);
396  ArrayType sample(30);
397  ArrayType res(30);
398  ArrayType v(30);
399 
400  alpha << 0.01, 0.01, 0.01, 0.01, 0.01, 0.1, 0.1, 0.1, 0.1, 0.1, 1.0, 1.0, 1.0, 1.0, 1.0, 10.0, 10.0, 10.0, 10.0,
401  10.0, 100.0, 100.0, 100.0, 100.0, 100.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0;
402 
403  sample << 1.25668890405e-26, 1.17549435082e-38, 1.20938905072e-05, 1.17549435082e-38, 1.17549435082e-38,
404  5.66572070696e-16, 0.0132865061065, 0.0200034203853, 6.29263709118e-17, 1.37160367764e-06, 0.333412038288,
405  1.18135687766, 0.580629033777, 0.170631439426, 0.786686768458, 7.63873279537, 13.1944344379, 11.896042354,
406  10.5830172417, 10.5020942233, 92.8918587747, 95.003720371, 86.3715926467, 96.0330217672, 82.6389930677,
407  968.702906754, 969.463546828, 1001.79726022, 955.047416547, 1044.27458568;
408 
409  v << 7.42424742367e-23, 1.02004297287e-34, 0.0130155240738, 1.02004297287e-34, 1.02004297287e-34, 1.96505168277e-13,
410  0.525575786243, 0.713903991771, 2.32077561808e-14, 0.000179348049886, 0.635500453302, 1.27561284917,
411  0.878125852156, 0.41565819538, 1.03606488534, 0.885964824887, 1.16424049334, 1.10764479598, 1.04590810812,
412  1.04193666963, 0.965193152414, 0.976217589464, 0.93008035061, 0.98153216096, 0.909196397698, 0.98434963993,
413  0.984738050206, 1.00106492525, 0.97734200649, 1.02198794179;
414 
415  CALL_SUBTEST(res = gamma_sample_der_alpha(alpha, sample); verify_component_wise(res, v););
416  }
417 #endif // EIGEN_HAS_C99_MATH
418 }
abs
AnnoyingScalar abs(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:135
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
m1
Matrix3d m1
Definition: IOFormat.cpp:2
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
m2
MatrixType m2(n_dims)
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
b
Scalar * b
Definition: benchVecAdd.cpp:17
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:45
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:46
Eigen::ArrayBase::pow
EIGEN_DEVICE_FUNC const GlobalUnaryPowReturnType< Derived, ScalarExponent > pow(const Eigen::ArrayBase< Derived > &x, const ScalarExponent &exponent)
Definition: GlobalFunctions.h:137
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:13
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:151
a
const Scalar * a
Definition: level2_cplx_impl.h:32
VERIFY
#define VERIFY(a)
Definition: main.h:362
CALL_SUBTEST
#define CALL_SUBTEST(FUNC)
Definition: main.h:382
isnan
#define isnan(X)
Definition: main.h:109
CRBond_Bessel::eps
double eps
Definition: crbond_bessel.cc:24
Eigen::bfloat16_impl::exp
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 exp(const bfloat16 &a)
Definition: BFloat16.h:615
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:2019
Eigen::zeta
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:152
Eigen::igammac
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_igammac_op< typename Derived::Scalar >, const Derived, const ExponentDerived > igammac(const Eigen::ArrayBase< Derived > &a, const Eigen::ArrayBase< ExponentDerived > &x)
Definition: SpecialFunctionsArrayAPI.h:93
Eigen::igamma
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_igamma_op< typename Derived::Scalar >, const Derived, const ExponentDerived > igamma(const Eigen::ArrayBase< Derived > &a, const Eigen::ArrayBase< ExponentDerived > &x)
Definition: SpecialFunctionsArrayAPI.h:31
Eigen::polygamma
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_polygamma_op< typename DerivedX::Scalar >, const DerivedN, const DerivedX > polygamma(const Eigen::ArrayBase< DerivedN > &n, const Eigen::ArrayBase< DerivedX > &x)
Definition: SpecialFunctionsArrayAPI.h:113
Eigen::betainc
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const TensorCwiseTernaryOp< internal::scalar_betainc_op< typename XDerived::Scalar >, const ADerived, const BDerived, const XDerived > betainc(const Eigen::TensorBase< ADerived, ReadOnlyAccessors > &a, const Eigen::TensorBase< BDerived, ReadOnlyAccessors > &b, const Eigen::TensorBase< XDerived, ReadOnlyAccessors > &x)
Definition: TensorGlobalFunctions.h:26
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:83
Eigen::gamma_sample_der_alpha
EIGEN_STRONG_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)
Definition: SpecialFunctionsArrayAPI.h:75
Eigen::igamma_der_a
EIGEN_STRONG_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)
Definition: SpecialFunctionsArrayAPI.h:52
oomph::PseudoSolidHelper::Zero
double Zero
Definition: pseudosolid_node_update_elements.cc:35
oomph::SarahBL::epsilon
double epsilon
Definition: osc_ring_sarah_asymptotics.h:43
plotDoE.x
list x
Definition: plotDoE.py:28
test
Definition: indexed_view.cpp:20
verify_component_wise
void verify_component_wise(const X &x, const Y &y)
Definition: special_functions.cpp:26
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:217
zero
EIGEN_DONT_INLINE Scalar zero()
Definition: svd_common.h:232
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References a, abs(), alpha, b, Eigen::betainc(), CALL_SUBTEST, cols, CRBond_Bessel::eps, oomph::SarahBL::epsilon, Eigen::bfloat16_impl::exp(), Eigen::gamma_sample_der_alpha(), i, Eigen::igamma(), Eigen::igamma_der_a(), Eigen::igammac(), isnan, j, m1, m2(), n, Eigen::polygamma(), Eigen::ArrayBase< Derived >::pow(), Eigen::numext::q, res, rows, sqrt(), v, VERIFY, verify_component_wise(), VERIFY_IS_APPROX, plotDoE::x, zero(), oomph::PseudoSolidHelper::Zero, and Eigen::zeta().

◆ EIGEN_DECLARE_TEST()

EIGEN_DECLARE_TEST ( special_functions  )
420  {
421  CALL_SUBTEST_1(array_special_functions<ArrayXf>());
422  CALL_SUBTEST_2(array_special_functions<ArrayXd>());
423  // TODO(cantonios): half/bfloat16 don't have enough precision to reproduce results above.
424  // CALL_SUBTEST_3(array_special_functions<ArrayX<Eigen::half>>());
425  // CALL_SUBTEST_4(array_special_functions<ArrayX<Eigen::bfloat16>>());
426 }
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10

References CALL_SUBTEST_1, and CALL_SUBTEST_2.

◆ operator,()

template<typename XprType >
Eigen::CommaInitializer<XprType>& operator, ( Eigen::CommaInitializer< XprType > &  ci,
double  v 
)
21  {
22  return (ci, static_cast<typename XprType::Scalar>(v));
23 }

References v.

◆ operator<<()

template<typename Derived >
Eigen::CommaInitializer<Derived> operator<< ( Eigen::DenseBase< Derived > &  dense,
double  v 
)
16  {
17  return (dense << static_cast<typename Derived::Scalar>(v));
18 }

References v.

◆ verify_component_wise()

template<typename X , typename Y >
void verify_component_wise ( const X &  x,
const Y &  y 
)
26  {
27  for (Index i = 0; i < x.size(); ++i) {
28  if ((numext::isfinite)(y(i)))
29  VERIFY_IS_APPROX(x(i), y(i));
30  else if ((numext::isnan)(y(i)))
31  VERIFY((numext::isnan)(x(i)));
32  else
33  VERIFY_IS_EQUAL(x(i), y(i));
34  }
35 }
y
Scalar * y
Definition: level1_cplx_impl.h:128
isfinite
#define isfinite(X)
Definition: main.h:111
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:367

References i, isfinite, isnan, VERIFY, VERIFY_IS_APPROX, VERIFY_IS_EQUAL, plotDoE::x, and y.

Referenced by array_special_functions().