31 using namespace Eigen;
36 #ifndef DOXYGEN_SHOULD_SKIP_THIS
48 class CITransformFunction :
public CFunction
61 virtual ~CITransformFunction() {
SG_UNREF(m_f); }
76 return (*m_f)(gx)*dgx;
99 class CILTransformFunction :
public CFunction
114 virtual ~CILTransformFunction() {
SG_UNREF(m_f); }
129 return -(*m_f)(m_b-
CMath::sq(gx))*2*gx*dgx;
155 class CIUTransformFunction :
public CFunction
170 virtual ~CIUTransformFunction() {
SG_UNREF(m_f); }
185 return (*m_f)(m_a+
CMath::sq(gx))*2*gx*dgx;
206 class CTransformFunction :
public CFunction
223 virtual ~CTransformFunction() {
SG_UNREF(m_f); }
239 return (*m_f)(gx)*dgx;
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",
263 REQUIRE(rel_tol>0.0,
"Relative tolerance must be positive, but is %f\n",
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, "
271 typedef void TQuadGKEvaluationFunction(
CFunction* f,
275 TQuadGKEvaluationFunction* evaluate_quadgk;
299 tf=
new CITransformFunction(f);
300 evaluate_quadgk=&evaluate_quadgk15;
306 tf=
new CILTransformFunction(f, b);
307 evaluate_quadgk=&evaluate_quadgk15;
313 tf=
new CIUTransformFunction(f, a);
314 evaluate_quadgk=&evaluate_quadgk15;
320 tf=
new CTransformFunction(f, a, b);
321 evaluate_quadgk=&evaluate_quadgk21;
343 evaluate_quadgk(tf, subs, q_subs, err_subs);
363 while (err>tol && iter<max_iter)
370 (*subs)[2*i])/(tb-ta))
373 float64_t mid=((*subs)[2*i]+(*subs)[2*i+1])/2.0;
397 evaluate_quadgk(tf, subs, q_subs, err_subs);
415 SG_SWARNING(
"Error tolerance not met. Estimated error is equal to %g "
416 "after %d iterations\n", err, iter)
444 "The length of node array (%d) and weight array (%d) should be the same\n",
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 "
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")
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);
479 VectorXd eigen_hw=(eigen_subs.row(1)-eigen_subs.row(0))/2.0;
480 VectorXd eigen_center=eigen_subs.colwise().sum()/2.0;
483 MatrixXd x=eigen_hw*eigen_xgk.adjoint()+eigen_center*
484 (VectorXd::Ones(n)).adjoint();
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));
494 VectorXd eigen_q=((ygk*eigen_wgk.asDiagonal()).rowwise().sum()).cwiseProduct(
496 q->set_array(eigen_q.data(), eigen_q.size());
499 MatrixXd yg(ygk.rows(), ygk.cols()/2);
501 for (
index_t i=1, j=0; i<ygk.cols(); i+=2, j++)
502 yg.col(j)=ygk.col(i);
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());
513 "The length of node array (%d) and weight array (%d) should be the same\n",
520 generate_gauher20(xgh, wgh);
527 eigen_xgh = MatrixXd::Zero(n,1);
528 eigen_wgh = MatrixXd::Ones(n,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);
541 eigen_wgh = eig.eigenvectors().row(0).transpose().real().array().pow(2);
544 eigen_xgh = eig.eigenvalues().real()*sqrt(2.0);
552 "The length of node array (%d) and weight array (%d) should be the same\n",
554 REQUIRE(xgh.
vlen == 20,
"The length of xgh and wgh should be 20\n");
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
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
606 for (
index_t idx = 0; idx < n; idx++)
608 wgh[idx] = wgh_pre[idx];
609 xgh[idx] = xgh_pre[idx];
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
642 0.129484966168869693270611432679082,
643 0.279705391489276667901467771423780,
644 0.381830050505118944950369775488975,
645 0.417959183673469387755102040816327,
646 0.381830050505118944950369775488975,
647 0.279705391489276667901467771423780,
648 0.129484966168869693270611432679082
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
672 evaluate_quadgk(f, subs, q, err, n, xgk, wg, wgk);
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
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
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
747 evaluate_quadgk(f, subs, q, err, n, xgk, wg, wgk);
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");
761 q+=wgh[i]*(*f)(xgh[i]);
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
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
908 return evaluate_quadgh(f, n, xgh, wgh);
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
Class that contains certain methods related to numerical integration.
Class of a function of one variable.
static float64_t integrate_quadgh(CFunction *f)
static void generate_gauher20(SGVector< float64_t > xgh, SGVector< float64_t > wgh)
void set_array(T *p_array, int32_t p_num_elements, int32_t array_size)
all of classes and functions are contained in the shogun namespace
int32_t get_num_elements() const
static float64_t integrate_quadgh_customized(CFunction *f, SGVector< float64_t > xgh, SGVector< float64_t > wgh)