SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Integration.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (w) 2014 Wu Lin
8  * Written (W) 2013 Roman Votyakov
9  *
10  * The abscissae and weights for Gauss-Kronrod rules are taken form
11  * QUADPACK, which is in public domain.
12  * http://www.netlib.org/quadpack/
13  *
14  * See header file for which functions are adapted from GNU Octave,
15  * file quadgk.m: Copyright (C) 2008-2012 David Bateman under GPLv3
16  * http://www.gnu.org/software/octave/
17  *
18  * See header file for which functions are adapted from
19  * Gaussian Process Machine Learning Toolbox, file util/gauher.m,
20  * http://www.gaussianprocess.org/gpml/code/matlab/doc/
21  */
22 
24 
25 #ifdef HAVE_EIGEN3
26 
28 #include <shogun/lib/SGVector.h>
29 
30 using namespace shogun;
31 using namespace Eigen;
32 
33 namespace shogun
34 {
35 
36 #ifndef DOXYGEN_SHOULD_SKIP_THIS
37 
48 class CITransformFunction : public CFunction
49 {
50 public:
55  CITransformFunction(CFunction* f)
56  {
57  SG_REF(f);
58  m_f=f;
59  }
60 
61  virtual ~CITransformFunction() { SG_UNREF(m_f); }
62 
70  virtual float64_t operator() (float64_t x)
71  {
72  float64_t hx=1.0/(1.0-CMath::sq(x));
73  float64_t gx=x*hx;
74  float64_t dgx=(1.0+CMath::sq(x))*CMath::sq(hx);
75 
76  return (*m_f)(gx)*dgx;
77  }
78 
79 private:
81  CFunction* m_f;
82 };
83 
99 class CILTransformFunction : public CFunction
100 {
101 public:
107  CILTransformFunction(CFunction* f, float64_t b)
108  {
109  SG_REF(f);
110  m_f=f;
111  m_b=b;
112  }
113 
114  virtual ~CILTransformFunction() { SG_UNREF(m_f); }
115 
123  virtual float64_t operator() (float64_t x)
124  {
125  float64_t hx=1.0/(1.0+x);
126  float64_t gx=x*hx;
127  float64_t dgx=CMath::sq(hx);
128 
129  return -(*m_f)(m_b-CMath::sq(gx))*2*gx*dgx;
130  }
131 
132 private:
134  CFunction* m_f;
135 
137  float64_t m_b;
138 };
139 
155 class CIUTransformFunction : public CFunction
156 {
157 public:
163  CIUTransformFunction(CFunction* f, float64_t a)
164  {
165  SG_REF(f);
166  m_f=f;
167  m_a=a;
168  }
169 
170  virtual ~CIUTransformFunction() { SG_UNREF(m_f); }
171 
179  virtual float64_t operator() (float64_t x)
180  {
181  float64_t hx=1.0/(1.0-x);
182  float64_t gx=x*hx;
183  float64_t dgx=CMath::sq(hx);
184 
185  return (*m_f)(m_a+CMath::sq(gx))*2*gx*dgx;
186  }
187 
188 private:
190  CFunction* m_f;
191 
193  float64_t m_a;
194 };
195 
206 class CTransformFunction : public CFunction
207 {
208 public:
215  CTransformFunction(CFunction* f, float64_t a, float64_t b)
216  {
217  SG_REF(f);
218  m_f=f;
219  m_a=a;
220  m_b=b;
221  }
222 
223  virtual ~CTransformFunction() { SG_UNREF(m_f); }
224 
233  virtual float64_t operator() (float64_t x)
234  {
235  float64_t qw=(m_b-m_a)/4.0;
236  float64_t gx=qw*(x*(3.0-CMath::sq(x)))+(m_b+m_a)/2.0;
237  float64_t dgx=qw*3.0*(1.0-CMath::sq(x));
238 
239  return (*m_f)(gx)*dgx;
240  }
241 
242 private:
244  CFunction* m_f;
245 
247  float64_t m_a;
248 
250  float64_t m_b;
251 };
252 
253 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
254 
256  float64_t b, float64_t abs_tol, float64_t rel_tol, uint32_t max_iter,
257  index_t sn)
258 {
259  // check the parameters
260  REQUIRE(f, "Integrable function should not be NULL\n")
261  REQUIRE(abs_tol>0.0, "Absolute tolerance must be positive, but is %f\n",
262  abs_tol)
263  REQUIRE(rel_tol>0.0, "Relative tolerance must be positive, but is %f\n",
264  rel_tol)
265  REQUIRE(max_iter>0, "Maximum number of iterations must be greater than 0, "
266  "but is %d\n", max_iter)
267  REQUIRE(sn>0, "Initial number of subintervals must be greater than 0, "
268  "but is %d\n", sn)
269 
270  // integral evaluation function
271  typedef void TQuadGKEvaluationFunction(CFunction* f,
274 
275  TQuadGKEvaluationFunction* evaluate_quadgk;
276 
277  CFunction* tf;
278  float64_t ta;
279  float64_t tb;
280  float64_t q_sign;
281 
282  // negate integral value and swap a and b, if a>b
283  if (a>b)
284  {
285  ta=b;
286  tb=a;
287  q_sign=-1.0;
288  }
289  else
290  {
291  ta=a;
292  tb=b;
293  q_sign=1.0;
294  }
295 
296  // transform integrable function and domain of integration
297  if (a==-CMath::INFTY && b==CMath::INFTY)
298  {
299  tf=new CITransformFunction(f);
300  evaluate_quadgk=&evaluate_quadgk15;
301  ta=-1.0;
302  tb=1.0;
303  }
304  else if (a==-CMath::INFTY)
305  {
306  tf=new CILTransformFunction(f, b);
307  evaluate_quadgk=&evaluate_quadgk15;
308  ta=-1.0;
309  tb=0.0;
310  }
311  else if (b==CMath::INFTY)
312  {
313  tf=new CIUTransformFunction(f, a);
314  evaluate_quadgk=&evaluate_quadgk15;
315  ta=0.0;
316  tb=1.0;
317  }
318  else
319  {
320  tf=new CTransformFunction(f, a, b);
321  evaluate_quadgk=&evaluate_quadgk21;
322  ta=-1.0;
323  tb=1.0;
324  }
325 
326  // compute initial subintervals, by dividing domain [a, b] into sn
327  // parts
329 
330  // width of each subinterval
331  float64_t sw=(tb-ta)/sn;
332 
333  for (index_t i=0; i<sn; i++)
334  {
335  subs->push_back(ta+i*sw);
336  subs->push_back(ta+(i+1)*sw);
337  }
338 
339  // evaluate integrals on initial subintervals
342 
343  evaluate_quadgk(tf, subs, q_subs, err_subs);
344 
345  // compute value of integral and error on [a, b]
346  float64_t q=0.0;
347  float64_t err=0.0;
348 
349  for (index_t i=0; i<q_subs->get_num_elements(); i++)
350  q+=(*q_subs)[i];
351 
352  for (index_t i=0; i<err_subs->get_num_elements(); i++)
353  err+=(*err_subs)[i];
354 
355  // evaluate tolerance
356  float64_t tol=CMath::max(abs_tol, rel_tol*CMath::abs(q));
357 
358  // number of iterations
359  uint32_t iter=1;
360 
362 
363  while (err>tol && iter<max_iter)
364  {
365  // choose and bisect subintervals with estimated error, which
366  // is larger or equal to tolerance
367  for (index_t i=0; i<subs->get_num_elements()/2; i++)
368  {
369  if (CMath::abs((*err_subs)[i])>=tol*CMath::abs((*subs)[2*i+1]-
370  (*subs)[2*i])/(tb-ta))
371  {
372  // bisect subinterval
373  float64_t mid=((*subs)[2*i]+(*subs)[2*i+1])/2.0;
374 
375  new_subs->push_back((*subs)[2*i]);
376  new_subs->push_back(mid);
377  new_subs->push_back(mid);
378  new_subs->push_back((*subs)[2*i+1]);
379 
380  // subtract value of the integral and error on this
381  // subinterval from total value and error
382  q-=(*q_subs)[i];
383  err-=(*err_subs)[i];
384  }
385  }
386 
387  subs->set_array(new_subs->get_array(), new_subs->get_num_elements(),
388  new_subs->get_num_elements());
389 
390  new_subs->reset_array();
391 
392  // break if no new subintervals
393  if (!subs->get_num_elements())
394  break;
395 
396  // evaluate integrals on selected subintervals
397  evaluate_quadgk(tf, subs, q_subs, err_subs);
398 
399  for (index_t i=0; i<q_subs->get_num_elements(); i++)
400  q+=(*q_subs)[i];
401 
402  for (index_t i=0; i<err_subs->get_num_elements(); i++)
403  err+=(*err_subs)[i];
404 
405  // evaluate tolerance
406  tol=CMath::max(abs_tol, rel_tol*CMath::abs(q));
407 
408  iter++;
409  }
410 
411  SG_UNREF(new_subs);
412 
413  if (err>tol)
414  {
415  SG_SWARNING("Error tolerance not met. Estimated error is equal to %g "
416  "after %d iterations\n", err, iter)
417  }
418 
419  // clean up
420  SG_UNREF(subs);
421  SG_UNREF(q_subs);
422  SG_UNREF(err_subs);
423  SG_UNREF(tf);
424 
425  return q_sign*q;
426 }
427 
429 {
430  SG_REF(f);
431 
432  // evaluate integral using Gauss-Hermite 64-point rule
433  float64_t q=evaluate_quadgh64(f);
434 
435  SG_UNREF(f);
436 
437  return q;
438 }
439 
442 {
443  REQUIRE(xgh.vlen == wgh.vlen,
444  "The length of node array (%d) and weight array (%d) should be the same\n",
445  xgh.vlen, wgh.vlen);
446 
447  SG_REF(f);
448 
449  float64_t q=evaluate_quadgh(f, xgh.vlen, xgh.vector, wgh.vector);
450 
451  SG_UNREF(f);
452 
453  return q;
454 }
455 
456 void CIntegration::evaluate_quadgk(CFunction* f, CDynamicArray<float64_t>* subs,
458  float64_t* xgk, float64_t* wg, float64_t* wgk)
459 {
460  // check the parameters
461  REQUIRE(f, "Integrable function should not be NULL\n")
462  REQUIRE(subs, "Array of subintervals should not be NULL\n")
463  REQUIRE(!(subs->get_array_size()%2), "Size of the array of subintervals "
464  "should be even\n")
465  REQUIRE(q, "Array of values of integrals should not be NULL\n")
466  REQUIRE(err, "Array of errors should not be NULL\n")
467  REQUIRE(n%2, "Order of Gauss-Kronrod should be odd\n")
468  REQUIRE(xgk, "Gauss-Kronrod nodes should not be NULL\n")
469  REQUIRE(wgk, "Gauss-Kronrod weights should not be NULL\n")
470  REQUIRE(wg, "Gauss weights should not be NULL\n")
471 
472  // create eigen representation of subs, xgk, wg, wgk
473  Map<MatrixXd> eigen_subs(subs->get_array(), 2, subs->get_num_elements()/2);
474  Map<VectorXd> eigen_xgk(xgk, n);
475  Map<VectorXd> eigen_wg(wg, n/2);
476  Map<VectorXd> eigen_wgk(wgk, n);
477 
478  // compute half width and centers of each subinterval
479  VectorXd eigen_hw=(eigen_subs.row(1)-eigen_subs.row(0))/2.0;
480  VectorXd eigen_center=eigen_subs.colwise().sum()/2.0;
481 
482  // compute Gauss-Kronrod nodes x for each subinterval: x=hw*xgk+center
483  MatrixXd x=eigen_hw*eigen_xgk.adjoint()+eigen_center*
484  (VectorXd::Ones(n)).adjoint();
485 
486  // compute ygk=f(x)
487  MatrixXd ygk(x.rows(), x.cols());
488 
489  for (index_t i=0; i<ygk.rows(); i++)
490  for (index_t j=0; j<ygk.cols(); j++)
491  ygk(i,j)=(*f)(x(i,j));
492 
493  // compute value of definite integral on each subinterval
494  VectorXd eigen_q=((ygk*eigen_wgk.asDiagonal()).rowwise().sum()).cwiseProduct(
495  eigen_hw);
496  q->set_array(eigen_q.data(), eigen_q.size());
497 
498  // choose function values for Gauss nodes
499  MatrixXd yg(ygk.rows(), ygk.cols()/2);
500 
501  for (index_t i=1, j=0; i<ygk.cols(); i+=2, j++)
502  yg.col(j)=ygk.col(i);
503 
504  // compute error on each subinterval
505  VectorXd eigen_err=(((yg*eigen_wg.asDiagonal()).rowwise().sum()).cwiseProduct(
506  eigen_hw)-eigen_q).array().abs();
507  err->set_array(eigen_err.data(), eigen_err.size());
508 }
509 
510 void CIntegration::generate_gauher(SGVector<float64_t> xgh, SGVector<float64_t> wgh)
511 {
512  REQUIRE(xgh.vlen == wgh.vlen,
513  "The length of node array (%d) and weight array (%d) should be the same\n",
514  xgh.vlen, wgh.vlen);
515 
516  index_t n = xgh.vlen;
517 
518  if (n == 20)
519  {
520  generate_gauher20(xgh, wgh);
521  }
522  else
523  {
524  Map<VectorXd> eigen_xgh(xgh.vector, xgh.vlen);
525  Map<VectorXd> eigen_wgh(wgh.vector, wgh.vlen);
526 
527  eigen_xgh = MatrixXd::Zero(n,1);
528  eigen_wgh = MatrixXd::Ones(n,1);
529 
530  if (n > 1)
531  {
532  MatrixXd v = MatrixXd::Zero(n,n);
533 
534  //b = sqrt( (1:N-1)/2 )';
535  //[V,D] = eig( diag(b,1) + diag(b,-1) );
536  v.block(0, 1, n-1, n-1).diagonal() = (0.5*ArrayXd::LinSpaced(n-1,1,n-1)).sqrt();
537  v.block(1, 0, n-1, n-1).diagonal() = v.block(0, 1, n-1, n-1).diagonal();
538  EigenSolver<MatrixXd> eig(v);
539 
540  //w = V(1,:)'.^2
541  eigen_wgh = eig.eigenvectors().row(0).transpose().real().array().pow(2);
542 
543  //x = sqrt(2)*diag(D)
544  eigen_xgh = eig.eigenvalues().real()*sqrt(2.0);
545  }
546  }
547 }
548 
550 {
551  REQUIRE(xgh.vlen == wgh.vlen,
552  "The length of node array (%d) and weight array (%d) should be the same\n",
553  xgh.vlen, wgh.vlen);
554  REQUIRE(xgh.vlen == 20, "The length of xgh and wgh should be 20\n");
555 
556  static const index_t n = 20;
557  static float64_t wgh_pre[n]=
558  {
559  0.0000000000001257800672437920121938444754,
560  0.0000000002482062362315158465220083577413,
561  0.0000000612749025998290679114578012251502,
562  0.0000044021210902308611768963750310312832,
563  0.0001288262799619289543807260089991473251,
564  0.0018301031310804880686271545187082665507,
565  0.0139978374471010288959682554832397727296,
566  0.0615063720639768204967445797137770568952,
567  0.1617393339840000332507941038784338161349,
568  0.2607930634495548849471902030927594751120,
569  0.2607930634495547739248877405771054327488,
570  0.1617393339840003108065502601675689220428,
571  0.0615063720639767788633811562704067910090,
572  0.0139978374471010080792865437615546397865,
573  0.0018301031310804856833823750505985117343,
574  0.0001288262799619298488475183095403053812,
575  0.0000044021210902308865878847926600414553,
576  0.0000000612749025998294252534824241331057,
577  0.0000000002482062362315177593771748866178,
578  0.0000000000001257800672437921636551382778
579  };
580 
581  static float64_t xgh_pre[n]=
582  {
583  -7.6190485416797573137159815814811736345291,
584  -6.5105901570136559541879250900819897651672,
585  -5.5787388058932032564030123467091470956802,
586  -4.7345813340460569662582201999612152576447,
587  -3.9439673506573176275935566081898286938667,
588  -3.1890148165533904744961546384729444980621,
589  -2.4586636111723669806394809711491689085960,
590  -1.7452473208141270344384565760265104472637,
591  -1.0429453488027506935509336472023278474808,
592  -0.3469641570813560282893206476728664711118,
593  0.3469641570813561393116231101885205134749,
594  1.0429453488027513596847484222962521016598,
595  1.7452473208141265903492467259638942778111,
596  2.4586636111723669806394809711491689085960,
597  3.1890148165533904744961546384729444980621,
598  3.9439673506573162953259270580019801855087,
599  4.7345813340460569662582201999612152576447,
600  5.5787388058932014800461729464586824178696,
601  6.5105901570136532896526659897062927484512,
602  7.6190485416797573137159815814811736345291
603 
604  };
605 
606  for (index_t idx = 0; idx < n; idx++)
607  {
608  wgh[idx] = wgh_pre[idx];
609  xgh[idx] = xgh_pre[idx];
610  }
611 
612 }
613 
614 void CIntegration::evaluate_quadgk15(CFunction* f, CDynamicArray<float64_t>* subs,
616 {
617  static const index_t n=15;
618 
619  // Gauss-Kronrod nodes
620  static float64_t xgk[n]=
621  {
622  -0.991455371120812639206854697526329,
623  -0.949107912342758524526189684047851,
624  -0.864864423359769072789712788640926,
625  -0.741531185599394439863864773280788,
626  -0.586087235467691130294144838258730,
627  -0.405845151377397166906606412076961,
628  -0.207784955007898467600689403773245,
629  0.000000000000000000000000000000000,
630  0.207784955007898467600689403773245,
631  0.405845151377397166906606412076961,
632  0.586087235467691130294144838258730,
633  0.741531185599394439863864773280788,
634  0.864864423359769072789712788640926,
635  0.949107912342758524526189684047851,
636  0.991455371120812639206854697526329
637  };
638 
639  // Gauss weights
640  static float64_t wg[n/2]=
641  {
642  0.129484966168869693270611432679082,
643  0.279705391489276667901467771423780,
644  0.381830050505118944950369775488975,
645  0.417959183673469387755102040816327,
646  0.381830050505118944950369775488975,
647  0.279705391489276667901467771423780,
648  0.129484966168869693270611432679082
649  };
650 
651  // Gauss-Kronrod weights
652  static float64_t wgk[n]=
653  {
654  0.022935322010529224963732008058970,
655  0.063092092629978553290700663189204,
656  0.104790010322250183839876322541518,
657  0.140653259715525918745189590510238,
658  0.169004726639267902826583426598550,
659  0.190350578064785409913256402421014,
660  0.204432940075298892414161999234649,
661  0.209482141084727828012999174891714,
662  0.204432940075298892414161999234649,
663  0.190350578064785409913256402421014,
664  0.169004726639267902826583426598550,
665  0.140653259715525918745189590510238,
666  0.104790010322250183839876322541518,
667  0.063092092629978553290700663189204,
668  0.022935322010529224963732008058970
669  };
670 
671  // evaluate definite integral on each subinterval using Gauss-Kronrod rule
672  evaluate_quadgk(f, subs, q, err, n, xgk, wg, wgk);
673 }
674 
675 void CIntegration::evaluate_quadgk21(CFunction* f, CDynamicArray<float64_t>* subs,
677 {
678  static const index_t n=21;
679 
680  // Gauss-Kronrod nodes
681  static float64_t xgk[n]=
682  {
683  -0.995657163025808080735527280689003,
684  -0.973906528517171720077964012084452,
685  -0.930157491355708226001207180059508,
686  -0.865063366688984510732096688423493,
687  -0.780817726586416897063717578345042,
688  -0.679409568299024406234327365114874,
689  -0.562757134668604683339000099272694,
690  -0.433395394129247190799265943165784,
691  -0.294392862701460198131126603103866,
692  -0.148874338981631210884826001129720,
693  0.000000000000000000000000000000000,
694  0.148874338981631210884826001129720,
695  0.294392862701460198131126603103866,
696  0.433395394129247190799265943165784,
697  0.562757134668604683339000099272694,
698  0.679409568299024406234327365114874,
699  0.780817726586416897063717578345042,
700  0.865063366688984510732096688423493,
701  0.930157491355708226001207180059508,
702  0.973906528517171720077964012084452,
703  0.995657163025808080735527280689003
704  };
705 
706  // Gauss weights
707  static float64_t wg[n/2]=
708  {
709  0.066671344308688137593568809893332,
710  0.149451349150580593145776339657697,
711  0.219086362515982043995534934228163,
712  0.269266719309996355091226921569469,
713  0.295524224714752870173892994651338,
714  0.295524224714752870173892994651338,
715  0.269266719309996355091226921569469,
716  0.219086362515982043995534934228163,
717  0.149451349150580593145776339657697,
718  0.066671344308688137593568809893332
719  };
720 
721  // Gauss-Kronrod weights
722  static float64_t wgk[n]=
723  {
724  0.011694638867371874278064396062192,
725  0.032558162307964727478818972459390,
726  0.054755896574351996031381300244580,
727  0.075039674810919952767043140916190,
728  0.093125454583697605535065465083366,
729  0.109387158802297641899210590325805,
730  0.123491976262065851077958109831074,
731  0.134709217311473325928054001771707,
732  0.142775938577060080797094273138717,
733  0.147739104901338491374841515972068,
734  0.149445554002916905664936468389821,
735  0.147739104901338491374841515972068,
736  0.142775938577060080797094273138717,
737  0.134709217311473325928054001771707,
738  0.123491976262065851077958109831074,
739  0.109387158802297641899210590325805,
740  0.093125454583697605535065465083366,
741  0.075039674810919952767043140916190,
742  0.054755896574351996031381300244580,
743  0.032558162307964727478818972459390,
744  0.011694638867371874278064396062192
745  };
746 
747  evaluate_quadgk(f, subs, q, err, n, xgk, wg, wgk);
748 }
749 
750 float64_t CIntegration::evaluate_quadgh(CFunction* f, index_t n, float64_t* xgh,
751  float64_t* wgh)
752 {
753  // check the parameters
754  REQUIRE(f, "Integrable function should not be NULL\n");
755  REQUIRE(xgh, "Gauss-Hermite nodes should not be NULL\n");
756  REQUIRE(wgh, "Gauss-Hermite weights should not be NULL\n");
757 
758  float64_t q=0.0;
759 
760  for (index_t i=0; i<n; i++)
761  q+=wgh[i]*(*f)(xgh[i]);
762 
763  return q;
764 }
765 
766 float64_t CIntegration::evaluate_quadgh64(CFunction* f)
767 {
768  static const index_t n=64;
769 
770  // Gauss-Hermite nodes
771  static float64_t xgh[n]=
772  {
773  -10.52612316796054588332682628381528,
774  -9.895287586829539021204461477159608,
775  -9.373159549646721162545652439723862,
776  -8.907249099964769757295972885642943,
777  -8.477529083379863090564166344821916,
778  -8.073687285010225225858791140758144,
779  -7.68954016404049682844780422986949,
780  -7.321013032780949201189569363719477,
781  -6.965241120551107529242642193492688,
782  -6.620112262636027379036660108937914,
783  -6.284011228774828235418093195070243,
784  -5.955666326799486045344567180984366,
785  -5.634052164349972147249920483307154,
786  -5.318325224633270857323649515199378,
787  -5.007779602198768196443702627184136,
788  -4.701815647407499816097538015812822,
789  -4.399917168228137647767932535438923,
790  -4.101634474566656714970981238455522,
791  -3.806571513945360461165972000460225,
792  -3.514375935740906211539950586474333,
793  -3.224731291992035725848171110188419,
794  -2.93735082300462180968533902619139,
795  -2.651972435430635011005457785998431,
796  -2.368354588632401404111511265341516,
797  -2.086272879881762020832563302363221,
798  -1.805517171465544918903773574186889,
799  -1.525889140209863662948970133151528,
800  -1.24720015694311794069356453069359,
801  -0.9692694230711780167435414890191023,
802  -0.6919223058100445772682192875955947,
803  -0.4149888241210786845769291291996859,
804  -0.1383022449870097241150497679666744,
805  0.1383022449870097241150497679666744,
806  0.4149888241210786845769291291996859,
807  0.6919223058100445772682192875955947,
808  0.9692694230711780167435414890191023,
809  1.24720015694311794069356453069359,
810  1.525889140209863662948970133151528,
811  1.805517171465544918903773574186889,
812  2.086272879881762020832563302363221,
813  2.368354588632401404111511265341516,
814  2.651972435430635011005457785998431,
815  2.93735082300462180968533902619139,
816  3.224731291992035725848171110188419,
817  3.514375935740906211539950586474333,
818  3.806571513945360461165972000460225,
819  4.101634474566656714970981238455522,
820  4.399917168228137647767932535438923,
821  4.701815647407499816097538015812822,
822  5.007779602198768196443702627184136,
823  5.318325224633270857323649515199378,
824  5.634052164349972147249920483307154,
825  5.955666326799486045344567180984366,
826  6.284011228774828235418093195070243,
827  6.620112262636027379036660108937914,
828  6.965241120551107529242642193492688,
829  7.321013032780949201189569363719477,
830  7.68954016404049682844780422986949,
831  8.073687285010225225858791140758144,
832  8.477529083379863090564166344821916,
833  8.907249099964769757295972885642943,
834  9.373159549646721162545652439723862,
835  9.895287586829539021204461477159608,
836  10.52612316796054588332682628381528
837  };
838 
839  // Gauss-Hermite weights
840  static float64_t wgh[n]=
841  {
842  5.535706535856942820575463300987E-49,
843  1.6797479901081592186662883306299E-43,
844  3.4211380112557405043272218281457E-39,
845  1.557390624629763802309335380265E-35,
846  2.549660899112999256604766580441E-32,
847  1.92910359546496685030196877906707E-29,
848  7.8617977889259103690999914962788E-27,
849  1.911706883300642829958456965534449E-24,
850  2.982862784279851154478700702016E-22,
851  3.15225456650378141612134668341E-20,
852  2.35188471067581911695767591555844E-18,
853  1.28009339132243804163956329526337E-16,
854  5.218623726590847522957808513052588E-15,
855  1.628340730709720362084307081240893E-13,
856  3.95917776694772392723644586425458E-12,
857  7.61521725014545135331529567531937E-11,
858  1.1736167423215493435425064670822E-9,
859  1.465125316476109354926622003804004E-8,
860  1.495532936727247061102461692934817E-7,
861  1.258340251031184576157842180019028E-6,
862  8.7884992308503591814440474067043E-6,
863  5.125929135786274660821911412739621E-5,
864  2.509836985130624860823620179819094E-4,
865  0.001036329099507577663456741746283101,
866  0.00362258697853445876066812537162265,
867  0.01075604050987913704946517278667313,
868  0.0272031289536889184538348212614932,
869  0.0587399819640994345496889462518317,
870  0.1084983493061868406330258455060973,
871  0.1716858423490837020007279701237768,
872  0.2329947860626780466505660293325675,
873  0.2713774249413039779456065084184279,
874  0.2713774249413039779456065084184279,
875  0.2329947860626780466505660293325675,
876  0.1716858423490837020007279701237768,
877  0.1084983493061868406330258455060973,
878  0.0587399819640994345496889462518317,
879  0.0272031289536889184538348212614932,
880  0.01075604050987913704946517278667313,
881  0.00362258697853445876066812537162265,
882  0.001036329099507577663456741746283101,
883  2.509836985130624860823620179819094E-4,
884  5.125929135786274660821911412739621E-5,
885  8.7884992308503591814440474067043E-6,
886  1.258340251031184576157842180019028E-6,
887  1.495532936727247061102461692934817E-7,
888  1.465125316476109354926622003804004E-8,
889  1.1736167423215493435425064670822E-9,
890  7.61521725014545135331529567531937E-11,
891  3.95917776694772392723644586425458E-12,
892  1.628340730709720362084307081240893E-13,
893  5.218623726590847522957808513052588E-15,
894  1.28009339132243804163956329526337E-16,
895  2.35188471067581911695767591555844E-18,
896  3.15225456650378141612134668341E-20,
897  2.982862784279851154478700702016E-22,
898  1.911706883300642829958456965534449E-24,
899  7.8617977889259103690999914962788E-27,
900  1.92910359546496685030196877906707E-29,
901  2.549660899112999256604766580441E-32,
902  1.557390624629763802309335380265E-35,
903  3.4211380112557405043272218281457E-39,
904  1.6797479901081592186662883306299E-43,
905  5.535706535856942820575463300987E-49
906  };
907 
908  return evaluate_quadgh(f, n, xgh, wgh);
909 }
910 }
911 
912 #endif /* HAVE_EIGEN3 */

SHOGUN Machine Learning Toolbox - Documentation