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

SHOGUN Machine Learning Toolbox - Documentation