cubature.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2015-2020, Andrew W. Steiner
5 
6  This file is part of O2scl. It has been adapted from cubature
7  written by Steven G. Johnson.
8 
9  O2scl is free software; you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation; either version 3 of the License, or
12  (at your option) any later version.
13 
14  O2scl is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
21 
22  -------------------------------------------------------------------
23 */
24 /* Adaptive multidimensional integration of a vector of integrands.
25  *
26  * Copyright (c) 2005-2013 Steven G. Johnson
27  *
28  * Portions (see comments) based on HIntLib (also distributed under
29  * the GNU GPL, v2 or later), copyright (c) 2002-2005 Rudolf Schuerer.
30  * (http://www.cosy.sbg.ac.at/~rschuer/hintlib/)
31  *
32  * Portions (see comments) based on GNU GSL (also distributed under
33  * the GNU GPL, v2 or later), copyright (c) 1996-2000 Brian Gough.
34  * (http://www.gnu.org/software/gsl/)
35  *
36  * This program is free software; you can redistribute it and/or modify
37  * it under the terms of the GNU General Public License as published by
38  * the Free Software Foundation; either version 2 of the License, or
39  * (at your option) any later version.
40  *
41  * This program is distributed in the hope that it will be useful,
42  * but WITHOUT ANY WARRANTY; without even the implied warranty of
43  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
44  * GNU General Public License for more details.
45  *
46  * You should have received a copy of the GNU General Public License
47  * along with this program; if not, write to the Free Software
48  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
49  *
50  */
51 /** \file cubature.h
52  \brief File for definitions of \ref o2scl::inte_hcubature and
53  \ref o2scl::inte_pcubature
54 */
55 #ifndef O2SCL_CUBATURE_H
56 #define O2SCL_CUBATURE_H
57 
58 #include <cmath>
59 #include <functional>
60 #include <boost/numeric/ublas/vector.hpp>
61 
62 #ifndef O2SCL_CLENCURT_H
63 #define O2SCL_CLENCURT_H
64 #include <o2scl/clencurt.h>
65 #endif
66 #include <o2scl/err_hnd.h>
67 #include <o2scl/vector.h>
68 
69 namespace o2scl {
70 
71  /** \brief Base class for integration routines from the
72  Cubature library
73  */
75 
76  public:
77 
78  /** \brief Different ways of measuring the absolute and
79  relative error
80 
81  Error estimates given a vector e of error estimates in the
82  individual components of a vector v of integrands. These are
83  all equivalent when there is only a single integrand.
84  */
85  typedef enum {
86  /** \brief individual relerr criteria in each component */
88  /** \brief paired L2 norms of errors in each component,
89  mainly for integrating vectors of complex numbers */
91  /** \brief abserr is L_2 norm |e|, and relerr is |e|/|v| */
93  /** \brief abserr is L_1 norm |e|, and relerr is |e|/|v| */
95  /** \brief abserr is \f$ L_{\infty} \f$ norm |e|, and
96  relerr is |e|/|v| */
98  } error_norm;
99 
100  };
101 
102  /** \brief Adaptive multidimensional integration on hyper-rectangles
103  using cubature rules from the Cubature library
104 
105  This class is experimental.
106 
107  \hline
108  \b Documentation \b adapted \b from \b Cubature
109 
110  A cubature rule takes a function and a hypercube and evaluates
111  the function at a small number of points, returning an estimate
112  of the integral as well as an estimate of the error, and also
113  a suggested dimension of the hypercube to subdivide.
114 
115  Given such a rule, the adaptive integration is simple:
116 
117  1) Evaluate the cubature rule on the hypercube(s).
118  Stop if converged.
119 
120  2) Pick the hypercube with the largest estimated error,
121  and divide it in two along the suggested dimension.
122 
123  3) Goto (1).
124 
125  The basic algorithm is based on the adaptive cubature described
126  in \ref Genz80 and subsequently extended to integrating a vector
127  of integrands in \ref Berntsen91 .
128 
129  \hline
130 
131  */
132  template<class func_t>
134 
135  public:
136 
138  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
139 
140  protected:
141 
142  /** \brief A value and error
143  */
144  class esterr {
145  public:
146  esterr() {
147  val=0.0;
148  err=0.0;
149  }
150  esterr(const esterr &e) {
151  val=e.val;
152  err=e.err;
153  }
154  esterr &operator=(const esterr &e) {
155  if (this!=&e) {
156  val=e.val;
157  err=e.err;
158  }
159  return *this;
160  }
161  /** \brief Desc */
162  double val;
163  /** \brief Desc */
164  double err;
165  };
166 
167  /** \brief Return the maximum error from \c ee
168  */
169  double errMax(const std::vector<esterr> &ee) {
170  double errmax = 0.0;
171  for (size_t k = 0; k < ee.size(); ++k) {
172  if (ee[k].err > errmax) errmax = ee[k].err;
173  }
174  return errmax;
175  }
176 
177  /** \brief Desc
178  */
179  class hypercube {
180  public:
181  hypercube() {
182  dim=0;
183  vol=0.0;
184  }
185  hypercube(const hypercube &e) {
186  dim=e.dim;
187  data=e.data;
188  vol=e.vol;
189  }
190  hypercube &operator=(const hypercube &e) {
191  if (this!=&e) {
192  dim=e.dim;
193  data=e.data;
194  vol=e.vol;
195  }
196  return *this;
197  }
198  /** \brief Desc */
199  size_t dim;
200  /** \brief length 2*dim = center followed by half-widths */
201  std::vector<double> data;
202  /** \brief cache volume = product of widths */
203  double vol;
204  };
205 
206  /** \brief Desc
207  */
208  double compute_vol(const hypercube &h) {
209  double vol = 1;
210  for (size_t i = 0; i < h.dim; ++i) {
211  vol *= 2 * h.data[i + h.dim];
212  }
213  return vol;
214  }
215 
216  /** \brief Desc
217  */
218  template<class vec_t>
219  void make_hypercube(size_t dim, const vec_t &center,
220  const vec_t &halfwidth, hypercube &h) {
221 
222  h.dim = dim;
223  h.data.resize(dim*2);
224  h.vol = 0;
225  for (size_t i = 0; i < dim; ++i) {
226  h.data[i] = center[i];
227  h.data[i + dim] = halfwidth[i];
228  }
229  h.vol = compute_vol(h);
230  return;
231  }
232 
233  /** \brief Desc
234  */
235  void make_hypercube2(size_t dim, const std::vector<double> &dat,
236  hypercube &h) {
237 
238  h.dim = dim;
239  h.data.resize(dim*2);
240  h.vol = 0;
241  for (size_t i = 0; i < dim; ++i) {
242  h.data[i] = dat[i];
243  h.data[i + dim] = dat[i+dim];
244  }
245  h.vol = compute_vol(h);
246  return;
247  }
248 
249  /** \brief Desc
250  */
251  template<class vec_t>
253  (size_t dim, const vec_t &xmin, const vec_t &xmax, hypercube &h) {
254 
255  make_hypercube(dim,xmin,xmax,h);
256  for (size_t i = 0; i < dim; ++i) {
257  h.data[i] = 0.5 * (xmin[i] + xmax[i]);
258  h.data[i + dim] = 0.5 * (xmax[i] - xmin[i]);
259  }
260  h.vol = compute_vol(h);
261  return;
262  }
263 
264  /** \brief Desc
265  */
267  h.data.clear();
268  h.dim = 0;
269  return;
270  }
271 
272  /** \brief Desc
273  */
274  class region {
275  public:
276  region() {
277  splitDim=0;
278  fdim=0;
279  errmax=0.0;
280  }
281  region(const region &e) {
282  h=e.h;
283  splitDim=e.splitDim;
284  fdim=e.fdim;
285  ee=e.ee;
286  errmax=e.errmax;
287  }
288  region &operator=(const region &e) {
289  if (this!=&e) {
290  h=e.h;
291  splitDim=e.splitDim;
292  fdim=e.fdim;
293  ee=e.ee;
294  errmax=e.errmax;
295  }
296  return *this;
297  }
298  /** \brief Desc */
300  /** \brief Desc */
301  size_t splitDim;
302  /** \brief dimensionality of vector integrand */
303  size_t fdim;
304  /** \brief array of length fdim */
305  std::vector<esterr> ee;
306  /** \brief max ee[k].err */
307  double errmax;
308  };
309 
310  /** \brief Desc
311  */
312  void make_region(const hypercube &h, size_t fdim, region &R) {
313 
314  make_hypercube2(h.dim, h.data,R.h);
315  R.splitDim = 0;
316  R.fdim = fdim;
317  R.ee.resize(fdim);
318  R.errmax = HUGE_VAL;
319 
320  return;
321  }
322 
323  /** \brief Desc
324  */
325  int cut_region(region &R, region &R2) {
326 
327  size_t d = R.splitDim, dim = R.h.dim;
328  R2=R;
329  R.h.data[d + dim] *= 0.5;
330  R.h.vol *= 0.5;
331  make_hypercube2(dim, R.h.data,R2.h);
332  R.h.data[d] -= R.h.data[d + dim];
333  R2.h.data[d] += R.h.data[d + dim];
334  R2.ee.resize(R2.fdim);
335  return 0;
336  }
337 
338  /** \brief Desc
339  */
340  class rule {
341 
342  public:
343  rule() {
344  dim=0;
345  fdim=0;
346  num_points=0;
347  num_regions=0;
348  vals_ix=0;
349  }
350  rule(const rule &e) {
351  dim=e.dim;
352  fdim=e.fdim;
355  pts=e.pts;
356  vals_ix=e.vals_ix;
357  }
358  rule &operator=(const rule &e) {
359  if (this!=&e) {
360  dim=e.dim;
361  fdim=e.fdim;
364  pts=e.pts;
365  vals_ix=e.vals_ix;
366  }
367  }
368 
369  /** \brief The dimensionality */
370  size_t dim;
371  /** \brief The number of functions */
372  size_t fdim;
373  /** \brief The number of evaluation points */
374  size_t num_points;
375  /** \brief The max number of regions evaluated at once */
376  size_t num_regions;
377  /** \brief points to eval: num_regions * num_points * dim */
379  /** \brief num_regions * num_points * fdim */
380  size_t vals_ix;
381  };
382 
383  /** \brief Desc
384  */
385  void alloc_rule_pts(rule &r, size_t num_regions) {
386  if (num_regions > r.num_regions) {
387  r.num_regions = 0;
388 
389  /* Allocate extra so that repeatedly calling alloc_rule_pts
390  with growing num_regions only needs a logarithmic number of
391  allocations
392  */
393  num_regions *= 2;
394 
395  r.pts.resize((num_regions
396  * r.num_points * (r.dim + r.fdim)));
397  r.vals_ix=num_regions * r.num_points * r.dim;
398  r.num_regions = num_regions;
399  }
400  return;
401  }
402 
403  /** \brief Desc
404  */
405  void make_rule(size_t dim, size_t fdim, size_t num_points, rule &r) {
406 
407  r.vals_ix=0;
408  r.num_regions = 0;
409  r.dim = dim;
410  r.fdim = fdim;
411  r.num_points = num_points;
412  return;
413  }
414 
415  /** \brief Desc
416 
417  \note All regions must have same fdim
418  */
419  int eval_regions(size_t nR, std::vector<region> &R, func_t &f, rule &r) {
420 
421  size_t iR;
422  if (nR == 0) {
423  /* nothing to evaluate */
424  return o2scl::success;
425  }
426  if (r.dim==1) {
427  if (rule15gauss_evalError(r, R[0].fdim, f, nR, R)) {
428  return o2scl::gsl_failure;
429  }
430  } else {
431  if (rule75genzmalik_evalError(r, R[0].fdim, f, nR, R)) {
432  return o2scl::gsl_failure;
433  }
434  }
435  for (iR = 0; iR < nR; ++iR) {
436  R[iR].errmax = errMax(R[iR].ee);
437  }
438  return o2scl::success;
439  }
440 
441  /** \brief Functions to loop over points in a hypercube.
442 
443  Based on orbitrule.cpp in HIntLib-0.0.10
444 
445  ls0 returns the least-significant 0 bit of n (e.g. it returns
446  0 if the LSB is 0, it returns 1 if the 2 LSBs are 01, etcetera).
447  */
448  size_t ls0(size_t n) {
449 
450 #if defined(__GNUC__) && \
451  ((__GNUC__ == 3 && __GNUC_MINOR__ >= 4) || __GNUC__ > 3)
452  return __builtin_ctz(~n); /* gcc builtin for version >= 3.4 */
453 #else
454  const size_t bits[256] = {
455  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
456  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
457  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
458  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6,
459  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
460  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
461  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
462  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 7,
463  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
464  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
465  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
466  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6,
467  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
468  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
469  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
470  0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 8,
471  };
472  size_t bit = 0;
473  while ((n & 0xff) == 0xff) {
474  n >>= 8;
475  bit += 8;
476  }
477  return bit + bits[n & 0xff];
478 #endif
479  }
480 
481  /** \brief Evaluate the integration points for all \f$ 2^n \f$ points
482  (+/-r,...+/-r)
483 
484  A Gray-code ordering is used to minimize the number of
485  coordinate updates in p, although this doesn't matter as much
486  now that we are saving all pts.
487  */
488  void evalR_Rfs2(ubvector &pts, size_t pts_ix, size_t dim,
489  std::vector<double> &p, size_t p_ix,
490  const std::vector<double> &c, size_t c_ix,
491  const std::vector<double> &r, size_t r_ix) {
492 
493  size_t i;
494  /* 0/1 bit = +/- for corresponding element of r[] */
495  size_t signs = 0;
496 
497  /* We start with the point where r is ADDed in every coordinate
498  (this implies signs=0). */
499  for (i = 0; i < dim; ++i) {
500  p[p_ix+i] = c[c_ix+i] + r[r_ix+i];
501  }
502 
503  /* Loop through the points in Gray-code ordering */
504  for (i = 0;; ++i) {
505  size_t mask, d;
506 
507  for(size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
508  pts_ix+=dim;
509 
510  d = ls0(i); /* which coordinate to flip */
511  if (d >= dim) {
512  break;
513  }
514 
515  /* flip the d-th bit and add/subtract r[d] */
516  mask = 1U << d;
517  signs ^= mask;
518  p[p_ix+d] = (signs & mask) ? c[c_ix+d] - r[r_ix+d] :
519  c[c_ix+d] + r[r_ix+d];
520  }
521  return;
522  }
523 
524  /** \brief Desc
525  */
526  void evalRR0_0fs2(ubvector &pts, size_t pts_ix, size_t dim,
527  std::vector<double> &p, size_t p_ix,
528  const std::vector<double> &c, size_t c_ix,
529  const std::vector<double> &r, size_t r_ix) {
530 
531  for (size_t i = 0; i < dim - 1; ++i) {
532  p[p_ix+i] = c[c_ix+i] - r[r_ix+i];
533  for (size_t j = i + 1; j < dim; ++j) {
534  p[p_ix+j] = c[c_ix+j] - r[r_ix+j];
535  for(size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
536  pts_ix+=dim;
537  p[p_ix+i] = c[c_ix+i] + r[r_ix+i];
538  for(size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
539  pts_ix+=dim;
540  p[p_ix+j] = c[c_ix+j] + r[r_ix+j];
541  for(size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
542  pts_ix+=dim;
543  p[p_ix+i] = c[c_ix+i] - r[r_ix+i];
544  for(size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
545  pts_ix+=dim;
546 
547  // Done with j -> Restore p[j]
548  p[p_ix+j] = c[c_ix+j];
549  }
550  // Done with i -> Restore p[i]
551  p[p_ix+i] = c[c_ix+i];
552  }
553  return;
554  }
555 
556  /** \brief Desc
557  */
558  void evalR0_0fs4d2
559  (ubvector &pts, size_t pts_ix, size_t dim,
560  std::vector<double> &p, size_t p_ix,
561  const std::vector<double> &c, size_t c_ix,
562  const std::vector<double> &r1, size_t r1_ix,
563  const std::vector<double> &r2, size_t r2_ix) {
564 
565  for(size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
566  pts_ix+=dim;
567 
568  for (size_t i = 0; i < dim; i++) {
569  p[p_ix+i] = c[c_ix+i] - r1[r1_ix+i];
570  for(size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
571  pts_ix+=dim;
572 
573  p[p_ix+i] = c[c_ix+i] + r1[r1_ix+i];
574  for(size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
575  pts_ix+=dim;
576 
577  p[p_ix+i] = c[c_ix+i] - r2[r2_ix+i];
578  for(size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
579  pts_ix+=dim;
580 
581  p[p_ix+i] = c[c_ix+i] + r2[r2_ix+i];
582  for(size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
583  pts_ix+=dim;
584 
585  p[p_ix+i] = c[c_ix+i];
586  }
587  return;
588  }
589 
590  /** \brief Desc
591  */
592  size_t num0_0(size_t dim) { return 1; }
593  /** \brief Desc
594  */
595  size_t numR0_0fs(size_t dim) { return 2*dim; }
596  /** \brief Desc
597  */
598  size_t numRR0_0fs(size_t dim) { return 2*dim*(dim-1); }
599  /** \brief Desc
600  */
601  size_t numR_Rfs(size_t dim) { return 1U << dim; }
602 
603  /** \brief Desc
604 
605  Based on rule75genzmalik.cpp in HIntLib-0.0.10: An embedded
606  cubature rule of degree 7 (embedded rule degree 5)
607  from \ref Genz83.
608  */
609  class rule75genzmalik : public rule {
610 
611  public:
612 
613  rule75genzmalik() : rule() {
614  }
615  rule75genzmalik(const rule75genzmalik &e) {
616  this->dim=e.dim;
617  this->fdim=e.fdim;
618  this->num_points=e.num_points;
619  this->num_regions=e.num_regions;
620  this->pts=e.pts;
621  this->vals_ix=e.vals_ix;
622  p=e.p;
623  weight1=e.weight1;
624  weight3=e.weight3;
625  weight5=e.weight5;
626  weightE1=e.weightE1;
627  weightE3=e.weightE3;
628  }
629  rule75genzmalik &operator=(const rule75genzmalik &e) {
630  if (this!=&e) {
631  this->dim=e.dim;
632  this->fdim=e.fdim;
633  this->num_points=e.num_points;
634  this->num_regions=e.num_regions;
635  this->pts=e.pts;
636  this->vals_ix=e.vals_ix;
637  p=e.p;
638  weight1=e.weight1;
639  weight3=e.weight3;
640  weight5=e.weight5;
641  weightE1=e.weightE1;
642  weightE3=e.weightE3;
643  }
644  }
645 
646  /** \brief Desc */
647  std::vector<double> p;
648 
649  /** \brief dimension-dependent constants */
650  double weight1;
651  /** \brief Desc */
652  double weight3;
653  /** \brief Desc */
654  double weight5;
655  /** \brief Desc */
656  double weightE1;
657  /** \brief Desc */
658  double weightE3;
659  };
660 
661 #ifdef O2SCL_NEVER_DEFINED
662  }{
663 #endif
664 
665  /** \brief Desc
666  */
667  int rule75genzmalik_evalError
668  (rule &runder, size_t fdim, func_t &f, size_t nR,
669  std::vector<region> &R) {
670 
671  /* lambda2 = sqrt(9/70), lambda4 = sqrt(9/10), lambda5 = sqrt(9/19) */
672  const double lambda2 = 0.3585685828003180919906451539079374954541;
673  const double lambda4 = 0.9486832980505137995996680633298155601160;
674  const double lambda5 = 0.6882472016116852977216287342936235251269;
675  const double weight2 = 980. / 6561.;
676  const double weight4 = 200. / 19683.;
677  const double weightE2 = 245. / 486.;
678  const double weightE4 = 25. / 729.;
679  const double ratio = (lambda2 * lambda2) / (lambda4 * lambda4);
680 
681  rule75genzmalik *r = (rule75genzmalik *) (&runder);
682  size_t i, j, iR, dim = runder.dim;
683  size_t npts = 0;
684  double *vals;
685  ubvector &pts2=runder.pts;
686 
687  alloc_rule_pts(runder, nR);
688  vals = &(runder.pts[runder.vals_ix]);
689 
690  for (iR = 0; iR < nR; ++iR) {
691  std::vector<double> &center2=R[iR].h.data;
692 
693  for (i = 0; i < dim; ++i) {
694  r->p[i] = center2[i];
695  }
696 
697  for (i = 0; i < dim; ++i) {
698  r->p[i+2*dim] = center2[i+dim] * lambda2;
699  }
700  for (i = 0; i < dim; ++i) {
701  r->p[i+dim] = center2[i+dim] * lambda4;
702  }
703 
704  /* Evaluate points in the center, in (lambda2,0,...,0) and
705  (lambda3=lambda4, 0,...,0). */
706  evalR0_0fs4d2(runder.pts,npts*dim, dim, r->p,0,center2,0,
707  r->p,2*dim,r->p,dim);
708  npts += num0_0(dim) + 2 * numR0_0fs(dim);
709 
710  /* Calculate points for (lambda4, lambda4, 0, ...,0) */
711  evalRR0_0fs2(runder.pts,npts*dim,dim,r->p,0,center2,0,r->p,dim);
712  npts += numRR0_0fs(dim);
713 
714  /* Calculate points for (lambda5, lambda5, ..., lambda5) */
715  for (i = 0; i < dim; ++i) {
716  r->p[i+dim] = center2[i+dim] * lambda5;
717  }
718  evalR_Rfs2(runder.pts,npts*dim,dim,r->p,0,center2,0,r->p,dim);
719  npts += numR_Rfs(dim);
720  }
721 
722  /* Evaluate the integrand function(s) at all the points */
723  if (f(dim, npts, &(pts2[0]), fdim, vals)) {
724  return o2scl::gsl_failure;
725  }
726 
727  /* we are done with the points, and so we can re-use the pts
728  array to store the maximum difference diff[i] in each dimension
729  for each hypercube */
730  //diff = pts;
731  for (i = 0; i < dim * nR; ++i) pts2[i] = 0;
732 
733  for (j = 0; j < fdim; ++j) {
734 
735  const double *v = vals + j;
736 
737  for (iR = 0; iR < nR; ++iR) {
738  double result, res5th;
739  double val0, sum2=0, sum3=0, sum4=0, sum5=0;
740  size_t k, k0 = 0;
741  /* accumulate j-th function values into j-th integrals
742  NOTE: this relies on the ordering of the eval functions
743  above, as well as on the internal structure of
744  the evalR0_0fs4d function */
745 
746  val0 = v[0]; /* central point */
747  k0 += 1;
748 
749  for (k = 0; k < dim; ++k) {
750  double v0 = v[fdim*(k0 + 4*k)];
751  double v1 = v[fdim*((k0 + 4*k) + 1)];
752  double v2 = v[fdim*((k0 + 4*k) + 2)];
753  double v3 = v[fdim*((k0 + 4*k) + 3)];
754 
755  sum2 += v0 + v1;
756  sum3 += v2 + v3;
757 
758  pts2[iR * dim + k] +=
759  fabs(v0 + v1 - 2*val0 - ratio * (v2 + v3 - 2*val0));
760  }
761 #ifdef O2SCL_NEVER_DEFINED
762  }{
763 #endif
764  k0 += 4*k;
765 
766  for (k = 0; k < numRR0_0fs(dim); ++k) {
767  sum4 += v[fdim*(k0 + k)];
768  }
769  k0 += k;
770 
771  for (k = 0; k < numR_Rfs(dim); ++k) {
772  sum5 += v[fdim*(k0 + k)];
773  }
774 
775  /* Calculate fifth and seventh order results */
776  result = R[iR].h.vol * (r->weight1 * val0 + weight2 * sum2 +
777  r->weight3 * sum3 + weight4 * sum4 +
778  r->weight5 * sum5);
779  res5th = R[iR].h.vol * (r->weightE1 * val0 + weightE2 * sum2 +
780  r->weightE3 * sum3 + weightE4 * sum4);
781 
782  R[iR].ee[j].val = result;
783  R[iR].ee[j].err = fabs(res5th - result);
784 
785  v += runder.num_points * fdim;
786  }
787  }
788 
789  /* figure out dimension to split: */
790  for (iR = 0; iR < nR; ++iR) {
791  double maxdiff = 0;
792  size_t dimDiffMax = 0;
793 
794  for (i = 0; i < dim; ++i) {
795  if (pts2[iR*dim + i] > maxdiff) {
796  maxdiff = pts2[iR*dim + i];
797  dimDiffMax = i;
798  }
799  }
800  R[iR].splitDim = dimDiffMax;
801  }
802  return o2scl::success;
803  }
804 
805 #ifdef O2SCL_NEVER_DEFINED
806  }{
807 #endif
808 
809  /** \brief Desc
810  */
811  void make_rule75genzmalik(size_t dim, size_t fdim,
812  rule75genzmalik &r) {
813 
814  if (dim < 2) {
815  O2SCL_ERR("this rule does not support 1d integrals",
817  }
818 
819  /* Because of the use of a bit-field in evalR_Rfs, we are limited
820  to be < 32 dimensions (or however many bits are in size_t).
821  This is not a practical limitation...long before you reach
822  32 dimensions, the Genz-Malik cubature becomes excruciatingly
823  slow and is superseded by other methods (e.g. Monte-Carlo). */
824  if (dim >= sizeof(size_t) * 8) {
825  O2SCL_ERR("this rule does not support large dims",
827  return;
828  }
829 
830  make_rule(dim, fdim,num0_0(dim) + 2 * numR0_0fs(dim)
831  + numRR0_0fs(dim) + numR_Rfs(dim),r);
832 
833  r.weight1=(12824.0-9120.0*dim+400.0*dim*dim)/19683.0;
834  r.weight3=(1820.0-400.0*dim)/19683.0;
835  r.weight5=6859.0/19683.0/((double)(1U << dim));
836  r.weightE1=(729.0-950.0*dim+50.0*dim*dim)/729.0;
837  r.weightE3=(265.0-100.0*dim)/1458.0;
838 
839  r.p.resize(dim*3);
840 
841  return;
842  }
843 
844  /** \brief 1d 15-point Gaussian quadrature rule, based on qk15.c
845  and qk.c in GNU GSL (which in turn is based on QUADPACK).
846  */
847  int rule15gauss_evalError
848  (rule &r, size_t fdim, func_t &f, size_t nR, std::vector<region> &R) {
849 
850  static const double cub_dbl_min=std::numeric_limits<double>::min();
851  static const double cub_dbl_eps=std::numeric_limits<double>::epsilon();
852 
853  /* Gauss quadrature weights and kronrod quadrature abscissae and
854  weights as evaluated with 80 decimal digit arithmetic by
855  L. W. Fullerton, Bell Labs, Nov. 1981. */
856  const size_t n = 8;
857  const double xgk[8] = { /* abscissae of the 15-point kronrod rule */
858  0.991455371120812639206854697526329,
859  0.949107912342758524526189684047851,
860  0.864864423359769072789712788640926,
861  0.741531185599394439863864773280788,
862  0.586087235467691130294144838258730,
863  0.405845151377397166906606412076961,
864  0.207784955007898467600689403773245,
865  0.000000000000000000000000000000000
866  /* xgk[1], xgk[3], ... abscissae of the 7-point gauss rule.
867  xgk[0], xgk[2], ... to optimally extend the 7-point gauss rule */
868  };
869  static const double wg[4] = { /* weights of the 7-point gauss rule */
870  0.129484966168869693270611432679082,
871  0.279705391489276667901467771423780,
872  0.381830050505118944950369775488975,
873  0.417959183673469387755102040816327
874  };
875  static const double wgk[8] = { /* weights of the 15-point kronrod rule */
876  0.022935322010529224963732008058970,
877  0.063092092629978553290700663189204,
878  0.104790010322250183839876322541518,
879  0.140653259715525918745189590510238,
880  0.169004726639267902826583426598550,
881  0.190350578064785409913256402421014,
882  0.204432940075298892414161999234649,
883  0.209482141084727828012999174891714
884  };
885  size_t j, k, iR;
886  size_t npts = 0;
887  double *pts, *vals;
888 
889  alloc_rule_pts(r, nR);
890  pts = &((r.pts)[0]);
891  vals = &(r.pts[r.vals_ix]);
892 
893  for (iR = 0; iR < nR; ++iR) {
894  const double center = R[iR].h.data[0];
895  const double halfwidth = R[iR].h.data[1];
896 
897  pts[npts++] = center;
898 
899  for (j = 0; j < (n - 1) / 2; ++j) {
900  int j2 = 2*j + 1;
901  double w = halfwidth * xgk[j2];
902  pts[npts++] = center - w;
903  pts[npts++] = center + w;
904  }
905  for (j = 0; j < n/2; ++j) {
906  int j2 = 2*j;
907  double w = halfwidth * xgk[j2];
908  pts[npts++] = center - w;
909  pts[npts++] = center + w;
910  }
911 
912  R[iR].splitDim = 0; /* no choice but to divide 0th dimension */
913  }
914 
915  if (f(1, npts, pts, fdim, vals)) {
916  return o2scl::gsl_failure;
917  }
918 
919  for (k = 0; k < fdim; ++k) {
920  const double *vk = vals + k;
921  for (iR = 0; iR < nR; ++iR) {
922  const double halfwidth = R[iR].h.data[1];
923  double result_gauss = vk[0] * wg[n/2 - 1];
924  double result_kronrod = vk[0] * wgk[n - 1];
925  double result_abs = fabs(result_kronrod);
926  double result_asc, mean, err;
927 
928  /* accumulate integrals */
929  npts = 1;
930  for (j = 0; j < (n - 1) / 2; ++j) {
931  int j2 = 2*j + 1;
932  double v = vk[fdim*npts] + vk[fdim*npts+fdim];
933  result_gauss += wg[j] * v;
934  result_kronrod += wgk[j2] * v;
935  result_abs += wgk[j2] * (fabs(vk[fdim*npts])
936  + fabs(vk[fdim*npts+fdim]));
937  npts += 2;
938  }
939  for (j = 0; j < n/2; ++j) {
940  int j2 = 2*j;
941  result_kronrod += wgk[j2] * (vk[fdim*npts]
942  + vk[fdim*npts+fdim]);
943  result_abs += wgk[j2] * (fabs(vk[fdim*npts])
944  + fabs(vk[fdim*npts+fdim]));
945  npts += 2;
946  }
947 
948  /* integration result */
949  R[iR].ee[k].val = result_kronrod * halfwidth;
950 
951  /* error estimate
952  (from GSL, probably dates back to QUADPACK
953  ... not completely clear to me why we don't just use
954  fabs(result_kronrod - result_gauss) * halfwidth */
955  mean = result_kronrod * 0.5;
956  result_asc = wgk[n - 1] * fabs(vk[0] - mean);
957  npts = 1;
958  for (j = 0; j < (n - 1) / 2; ++j) {
959  int j2 = 2*j + 1;
960  result_asc += wgk[j2] * (fabs(vk[fdim*npts]-mean)
961  + fabs(vk[fdim*npts+fdim]-mean));
962  npts += 2;
963  }
964  for (j = 0; j < n/2; ++j) {
965  int j2 = 2*j;
966  result_asc += wgk[j2] * (fabs(vk[fdim*npts]-mean)
967  + fabs(vk[fdim*npts+fdim]-mean));
968  npts += 2;
969  }
970  err = fabs(result_kronrod - result_gauss) * halfwidth;
971  result_abs *= halfwidth;
972  result_asc *= halfwidth;
973  if (result_asc != 0 && err != 0) {
974  double scale = pow((200 * err / result_asc), 1.5);
975  err = (scale < 1) ? result_asc * scale : result_asc;
976  }
977  if (result_abs > cub_dbl_min / (50 * cub_dbl_eps)) {
978  double min_err = 50 * cub_dbl_eps * result_abs;
979  if (min_err > err) err = min_err;
980  }
981  R[iR].ee[k].err = err;
982 
983  /* increment vk to point to next batch of results */
984  vk += 15*fdim;
985  }
986  }
987  return o2scl::success;
988  }
989 
990  /** \brief Desc
991  */
992  void make_rule15gauss(size_t dim, size_t fdim, rule &r) {
993 
994  if (dim != 1) {
995  O2SCL_ERR("this rule is only for 1d integrals.",o2scl::exc_esanity);
996  }
997 
998  return make_rule(dim,fdim,15,r);
999  }
1000 
1001  /** \name Binary heap implementation
1002 
1003  Based on \ref Cormen09 and used as a priority queue of
1004  regions to integrate.
1005  */
1006  //@{
1007  /** \brief Desc
1008  */
1010 
1011  /** \brief Desc
1012  */
1013  class heap {
1014  public:
1015  heap() {
1016  n=0;
1017  nalloc=0;
1018  fdim=0;
1019  }
1020  heap(const heap &e) {
1021  n=e.n;
1022  nalloc=e.nalloc;
1023  items=e.items;
1024  fdim=e.fdim;
1025  ee=e.ee;
1026  }
1027  heap &operator=(const heap &e) {
1028  if (this!=&e) {
1029  n=e.n;
1030  nalloc=e.nalloc;
1031  items=e.items;
1032  fdim=e.fdim;
1033  ee=e.ee;
1034  }
1035  return *this;
1036  }
1037  /** \brief Desc */
1038  size_t n;
1039  /** \brief Desc */
1040  size_t nalloc;
1041  /** \brief Desc */
1042  std::vector<heap_item> items;
1043  /** \brief Desc */
1044  size_t fdim;
1045  /** array of length fdim of the total integrand & error */
1046  std::vector<esterr> ee;
1047  };
1048 
1049  /** \brief Desc
1050  */
1051  void heap_resize(heap &h, size_t nalloc) {
1052 
1053  h.nalloc = nalloc;
1054  if (nalloc) {
1055  h.items.resize(nalloc);
1056  }
1057  return;
1058  }
1059 
1060  /** \brief Desc
1061  */
1062  heap heap_alloc(size_t nalloc, size_t fdim) {
1063 
1064  heap h;
1065  h.n = 0;
1066  h.nalloc = 0;
1067  h.fdim = fdim;
1068  h.ee.resize(fdim);
1069  for (size_t i = 0; i < fdim; ++i) {
1070  h.ee[i].val = 0.0;
1071  h.ee[i].err = 0.0;
1072  }
1073  heap_resize(h, nalloc);
1074  return h;
1075  }
1076 
1077  /** \brief Note that heap_free does not deallocate anything referenced by
1078  the items */
1079  void heap_free(heap &h) {
1080 
1081  h.n = 0;
1082  heap_resize(h, 0);
1083  h.fdim = 0;
1084  h.ee.clear();
1085  return;
1086  }
1087 
1088  /** \brief Desc
1089  */
1090  int heap_push(heap &h, heap_item &hi) {
1091 
1092  int insert;
1093  size_t fdim = h.fdim;
1094 
1095  for (size_t i = 0; i < fdim; ++i) {
1096  h.ee[i].val += hi.ee[i].val;
1097  h.ee[i].err += hi.ee[i].err;
1098  }
1099  insert = h.n;
1100  if (++(h.n) > h.nalloc) {
1101  heap_resize(h, h.n * 2);
1102  }
1103 
1104  while (insert) {
1105  int parent = (insert - 1) / 2;
1106  if (hi.errmax <= h.items[parent].errmax) {
1107  break;
1108  }
1109  h.items[insert] = h.items[parent];
1110  insert = parent;
1111  }
1112  h.items[insert] = hi;
1113  return o2scl::success;
1114  }
1115 
1116  /** \brief Desc
1117  */
1118  int heap_push_many(heap &h, size_t ni, heap_item *hi) {
1119  for (size_t i = 0; i < ni; ++i) {
1120  if (heap_push(h, hi[i])) return o2scl::gsl_failure;
1121  }
1122  return o2scl::success;
1123  }
1124 
1125  /** \brief Desc
1126  */
1128 
1129  heap_item ret;
1130  int i, n, child;
1131 
1132  if (!(h.n)) {
1133  O2SCL_ERR("Attempted to pop an empty heap in cubature.",
1135  }
1136 
1137  ret = h.items[0];
1138  h.items[i = 0] = h.items[n = --(h.n)];
1139 
1140  while ((child = i * 2 + 1) < n) {
1141 
1142  int largest;
1143  heap_item swap;
1144 
1145  if (h.items[child].errmax <= h.items[i].errmax) {
1146  largest = i;
1147  } else {
1148  largest = child;
1149  }
1150 
1151  if (++child < n && h.items[largest].errmax <
1152  h.items[child].errmax) {
1153  largest = child;
1154  }
1155  if (largest == i) {
1156  break;
1157  }
1158  swap = h.items[i];
1159  h.items[i] = h.items[largest];
1160  h.items[i = largest] = swap;
1161  }
1162 
1163  {
1164  size_t ii, fdim = h.fdim;
1165  for (ii = 0; ii < fdim; ++ii) {
1166  h.ee[ii].val -= ret.ee[ii].val;
1167  h.ee[ii].err -= ret.ee[ii].err;
1168  }
1169  }
1170  return ret;
1171  }
1172  //@}
1173 
1174  /** \brief Desc
1175  */
1176  int converged(size_t fdim, const std::vector<esterr> &ee,
1177  double reqAbsError, double reqRelError,
1178  error_norm norm) {
1179 
1180  size_t j;
1181 
1182  switch (norm) {
1183 
1184  case ERROR_INDIVIDUAL:
1185 
1186  for (j = 0; j < fdim; ++j) {
1187  if (ee[j].err > reqAbsError && ee[j].err >
1188  fabs(ee[j].val)*reqRelError) {
1189  return 0;
1190  }
1191  }
1192  return 1;
1193 
1194  case ERROR_PAIRED:
1195 
1196  for (j = 0; j+1 < fdim; j += 2) {
1197  double maxerr, serr, err, maxval, sval, val;
1198  /* scale to avoid overflow/underflow */
1199  maxerr = ee[j].err > ee[j+1].err ? ee[j].err : ee[j+1].err;
1200  maxval = ee[j].val > ee[j+1].val ? ee[j].val : ee[j+1].val;
1201  serr = maxerr > 0 ? 1/maxerr : 1;
1202  sval = maxval > 0 ? 1/maxval : 1;
1203  err=std::hypot(ee[j].err*serr,ee[j+1].err*serr)*maxerr;
1204  val=std::hypot(ee[j].val*sval,ee[j+1].val*sval)*maxval;
1205  if (err > reqAbsError && err > val*reqRelError) {
1206  return 0;
1207  }
1208  }
1209 
1210  /* fdim is odd, do last dimension individually */
1211  if (j < fdim) {
1212  if (ee[j].err > reqAbsError && ee[j].err >
1213  fabs(ee[j].val)*reqRelError) {
1214  return 0;
1215  }
1216  }
1217 
1218  return 1;
1219 
1220  case ERROR_L1:
1221  {
1222  double err = 0, val = 0;
1223  for (j = 0; j < fdim; ++j) {
1224  err += ee[j].err;
1225  val += fabs(ee[j].val);
1226  }
1227  return err <= reqAbsError || err <= val*reqRelError;
1228  }
1229 
1230  case ERROR_LINF:
1231  {
1232  double err = 0, val = 0;
1233  for (j = 0; j < fdim; ++j) {
1234  double absval = fabs(ee[j].val);
1235  if (ee[j].err > err) err = ee[j].err;
1236  if (absval > val) val = absval;
1237  }
1238  return err <= reqAbsError || err <= val*reqRelError;
1239  }
1240 
1241  case ERROR_L2:
1242  {
1243  double maxerr = 0, maxval = 0, serr, sval, err = 0, val = 0;
1244  /* scale values by 1/max to avoid overflow/underflow */
1245  for (j = 0; j < fdim; ++j) {
1246  double absval = fabs(ee[j].val);
1247  if (ee[j].err > maxerr) maxerr = ee[j].err;
1248  if (absval > maxval) maxval = absval;
1249  }
1250  serr = maxerr > 0 ? 1/maxerr : 1;
1251  sval = maxval > 0 ? 1/maxval : 1;
1252  for (j = 0; j < fdim; ++j) {
1253  err += (ee[j].err * serr)*(ee[j].err * serr);
1254  val += fabs(ee[j].val) * sval*fabs(ee[j].val) * sval;
1255  }
1256  err = sqrt(err) * maxerr;
1257  val = sqrt(val) * maxval;
1258  return err <= reqAbsError || err <= val*reqRelError;
1259  }
1260 
1261  }
1262  return 1; /* unreachable */
1263  }
1264 
1265  /** \brief Desc
1266  */
1267  template<class vec_t>
1268  int rulecubature(rule &r, size_t fdim, func_t &f,
1269  const hypercube &h, size_t maxEval,
1270  double reqAbsError, double reqRelError,
1271  error_norm norm, vec_t &val,
1272  vec_t &err, int parallel) {
1273 
1274  size_t numEval = 0;
1275  heap regions;
1276  size_t i, j;
1277  /* array of regions to evaluate */
1278  std::vector<region> R;
1279  size_t nR_alloc = 0;
1280  std::vector<esterr> ee(fdim);
1281 
1282  /* norm is irrelevant */
1283  if (fdim <= 1) norm = ERROR_INDIVIDUAL;
1284  /* invalid norm */
1285  if (norm < 0 || norm > ERROR_LINF) return o2scl::gsl_failure;
1286 
1287  regions = heap_alloc(1, fdim);
1288 
1289  nR_alloc = 2;
1290  R.resize(nR_alloc);
1291  make_region(h, fdim, R[0]);
1292  if (eval_regions(1, R, f, r) || heap_push(regions, R[0])) {
1293  heap_free(regions);
1294  //delete R;
1295  return o2scl::gsl_failure;
1296  }
1297  numEval += r.num_points;
1298 
1299  while (numEval < maxEval || !maxEval) {
1300 
1301  if (converged(fdim, regions.ee, reqAbsError, reqRelError, norm)) {
1302  break;
1303  }
1304 
1305  /* maximize potential parallelism */
1306  if (parallel) {
1307 
1308  /* Adapted from I. Gladwell, "Vectorization of one
1309  dimensional quadrature codes," pp. 230--238 in _Numerical
1310  Integration. Recent Developments, Software and
1311  Applications_, G. Fairweather and P. M. Keast, eds., NATO
1312  ASI Series C203, Dordrecht (1987), as described in J. M.
1313  Bull and T. L. Freeman, "Parallel Globally Adaptive
1314  Algorithms for Multi-dimensional Integration,"
1315  http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.6638
1316  (1994).
1317 
1318  Basically, this evaluates in one shot all regions that
1319  *must* be evaluated in order to reduce the error to the
1320  requested bound: the minimum set of largest-error regions
1321  whose errors push the total error over the bound.
1322 
1323  [Note: Bull and Freeman claim that the Gladwell approach
1324  is intrinsically inefficent because it "requires
1325  sorting", and propose an alternative algorithm that
1326  "only" requires three passes over the entire set of
1327  regions. Apparently, they didn't realize that one could
1328  use a heap data structure, in which case the time to pop
1329  K biggest-error regions out of N is only O(K log N), much
1330  better than the O(N) cost of the Bull and Freeman
1331  algorithm if K << N, and it is also much simpler.]
1332  */
1333  size_t nR = 0;
1334  for (j = 0; j < fdim; ++j) ee[j] = regions.ee[j];
1335  do {
1336 
1337  if (nR + 2 > nR_alloc) {
1338  nR_alloc = (nR + 2) * 2;
1339  R.resize(nR_alloc);
1340  }
1341  R[nR] = heap_pop(regions);
1342  for (j = 0; j < fdim; ++j) ee[j].err -= R[nR].ee[j].err;
1343  numEval += r.num_points * 2;
1344  nR += 2;
1345  if (converged(fdim, ee, reqAbsError, reqRelError, norm)) {
1346  /* other regions have small errs */
1347  break;
1348  }
1349  } while (regions.n > 0 && (numEval < maxEval || !maxEval));
1350 
1351  if (eval_regions(nR, R, f, r)
1352  || heap_push_many(regions, nR, &(R[0]))) {
1353  heap_free(regions);
1354  //delete R;
1355  return o2scl::gsl_failure;
1356  }
1357 
1358  // End of 'if (parallel)'
1359  } else {
1360 
1361  /* minimize number of function evaluations */
1362 
1363  /* get worst region */
1364  R[0] = heap_pop(regions);
1365  if (cut_region(R[0], R[1]) || eval_regions(2, R, f, r)
1366  || heap_push_many(regions, 2, &(R[0]))) {
1367  heap_free(regions);
1368  //delete R;
1369  return o2scl::gsl_failure;
1370  }
1371  numEval += r.num_points * 2;
1372  }
1373  }
1374 
1375  /* re-sum integral and errors */
1376  for (j = 0; j < fdim; ++j) {
1377  val[j] = 0.0;
1378  err[j] = 0.0;
1379  }
1380  for (i = 0; i < regions.n; ++i) {
1381  for (j = 0; j < fdim; ++j) {
1382  val[j] += regions.items[i].ee[j].val;
1383  err[j] += regions.items[i].ee[j].err;
1384  }
1385  {
1386  destroy_hypercube(regions.items[i].h);
1387  regions.items[i].ee.clear();
1388  }
1389  }
1390 
1391  heap_free(regions);
1392  //delete R;
1393 
1394  return o2scl::success;
1395  }
1396 
1397  /** \brief Desc
1398  */
1399  template<class vec_t>
1400  int cubature(size_t fdim, func_t &f, size_t dim, const vec_t &xmin,
1401  const vec_t &xmax, size_t maxEval, double reqAbsError,
1402  double reqRelError, error_norm norm, vec_t &val,
1403  vec_t &err, int parallel) {
1404 
1405  hypercube h;
1406 
1407  if (fdim == 0) {
1408  /* nothing to do */
1409  return o2scl::success;
1410  }
1411  if (dim == 0) {
1412  /* trivial integration */
1413  if (f(0, 1, &(xmin[0]), fdim, &(val[0]))) {
1414  return o2scl::gsl_failure;
1415  }
1416  for (size_t i = 0; i < fdim; ++i) err[i] = 0;
1417  return o2scl::success;
1418  }
1419  int status;
1420  if (dim==1) {
1421  rule r;
1422  make_rule15gauss(dim,fdim,r);
1423  make_hypercube_range(dim,xmin,xmax,h);
1424  status = rulecubature(r,fdim,f,h,maxEval,reqAbsError,
1425  reqRelError,norm,val,err,parallel);
1426  } else {
1427  rule75genzmalik r;
1428  make_rule75genzmalik(dim,fdim,r);
1429  make_hypercube_range(dim,xmin,xmax,h);
1430  status = rulecubature(r,fdim,f,h,maxEval,reqAbsError,
1431  reqRelError,norm,val,err,parallel);
1432  }
1433  destroy_hypercube(h);
1434 
1435  return status;
1436  }
1437 
1438  public:
1439 
1440  /// Desc
1442 
1443  inte_hcubature() {
1444  use_parallel=0;
1445  }
1446 
1447  /** \brief Desc
1448  */
1449  template<class vec_t>
1450  int integ(size_t fdim, func_t &f, size_t dim,
1451  const vec_t &xmin, const vec_t &xmax, size_t maxEval,
1452  double reqAbsError, double reqRelError, error_norm norm,
1453  vec_t &val, vec_t &err) {
1454 
1455  if (fdim == 0) {
1456  /* nothing to do */
1457  return o2scl::success;
1458  }
1459  return cubature(fdim,f,dim,xmin,xmax,maxEval,reqAbsError,
1460  reqRelError,norm,val,err,use_parallel);
1461 
1462  }
1463 
1464  };
1465 
1466 #ifdef O2SCL_NEVER_DEFINED
1467 }{
1468 #endif
1469 
1470  /** \brief Integration by p-adaptive cubature from the Cubature library
1471 
1472  This class is experimental.
1473 
1474  \hline
1475  \b Documentation \b adapted \b from \b Cubature
1476 
1477  This class performs adaptive integration by increasing the
1478  degree of the cubature rule rather than subdividing the domain,
1479  using products of Clenshaw-Curtis rules. This algorithm may be
1480  superior to Genz-Malik for smooth integrands lacking
1481  strongly-localized features, in moderate dimensions.
1482 
1483  \hline
1484  */
1485  template<class func_t, class vec_t, class vec_crange_t, class vec_range_t>
1487 
1488  protected:
1489 
1490  /** \brief Maximum integral dimension
1491  */
1492  static const size_t MAXDIM=20;
1493 
1494  /** \brief Cache of the values for the m[dim] grid.
1495 
1496  For adaptive cubature, thanks to the nesting of the C-C rules, we
1497  can re-use the values from coarser grids for finer grids, and the
1498  coarser grids are also used for error estimation.
1499 
1500  A grid is determined by an m[dim] array, where m[i] denotes
1501  2^(m[i]+1)+1 points in the i-th dimension.
1502 
1503  If mi < dim, then we only store the values corresponding to
1504  the difference between the m grid and the grid with m[mi] ->
1505  m[mi]-1. (m[mi]-1 == -1 corresponds to the trivial grid of one
1506  point in the center.)
1507  */
1508  class cache {
1509  public:
1510  cache() {
1511  m.resize(MAXDIM);
1512  mi=0;
1513  }
1514  cache(const cache &e) {
1515  m=e.m;
1516  mi=e.mi;
1517  val=e.val;
1518  }
1519  cache &operator=(const cache &e) {
1520  if (this!=&e) {
1521  m=e.m;
1522  mi=e.mi;
1523  val=e.val;
1524  }
1525  }
1526  /** \brief Desc */
1527  std::vector<size_t> m;
1528  /** \brief Desc */
1529  size_t mi;
1530  /** \brief Desc */
1531  vec_t val;
1532  };
1533 
1534  /** \brief Desc
1535 
1536  recursive loop over all cubature points for the given (m,mi)
1537  cache entry: add each point to the buffer buf, evaluating all
1538  at once whenever the buffer is full or when we are done
1539  */
1540  int compute_cacheval(const std::vector<size_t> &m, size_t mi, vec_t &val,
1541  size_t &vali, size_t fdim, func_t &f, size_t dim,
1542  size_t id, std::vector<double> &p, const vec_t &xmin,
1543  const vec_t &xmax, vec_t &buf, size_t nbuf,
1544  size_t &ibuf) {
1545 
1546  if (id == dim) {
1547  /* add point to buffer of points */
1548  for(size_t k=0;k<dim;k++) {
1549  buf[k+ibuf*dim]=p[k];
1550  }
1551  ibuf++;
1552  if (ibuf == nbuf) {
1553  const vec_crange_t buf2=o2scl::const_vector_range
1554  (((const vec_t &)buf),0,buf.size());
1555  vec_range_t val2=o2scl::vector_range(val,vali,val.size());
1556  /* flush buffer */
1557  if (f(dim, nbuf, &(buf[0]), fdim, &(val[0]) + vali)) {
1558  return o2scl::gsl_failure;
1559  }
1560  vali += ibuf * fdim;
1561  ibuf = 0;
1562  }
1563 
1564  } else {
1565 
1566  double c = (xmin[id] + xmax[id]) * 0.5;
1567  double r = (xmax[id] - xmin[id]) * 0.5;
1568  const double *x = clencurt_x
1569  + ((id == mi) ? (m[id] ? (1 << (m[id] - 1)) : 0) : 0);
1570  size_t i;
1571  size_t nx = (id == mi ? (m[id] ? (1 << (m[id] - 1)) : 1)
1572  : (1 << (m[id])));
1573  if (id != mi) {
1574  p[id] = c;
1575  if (compute_cacheval(m, mi, val, vali, fdim, f,
1576  dim, id + 1, p,
1577  xmin, xmax, buf, nbuf, ibuf)) {
1578  return o2scl::gsl_failure;
1579  }
1580  }
1581  for (i = 0; i < nx; ++i) {
1582  p[id] = c + r * x[i];
1583  if (compute_cacheval(m, mi, val, vali, fdim, f,
1584  dim, id + 1, p,
1585  xmin, xmax, buf, nbuf, ibuf)) {
1586  return o2scl::gsl_failure;
1587  }
1588  p[id] = c - r * x[i];
1589  if (compute_cacheval(m, mi, val, vali, fdim, f,
1590  dim, id + 1, p,
1591  xmin, xmax, buf, nbuf, ibuf)) {
1592  return o2scl::gsl_failure;
1593  }
1594  }
1595  }
1596  return o2scl::success;
1597  }
1598 
1599  /** \brief Desc
1600  */
1601  size_t num_cacheval(const std::vector<size_t> &m, size_t mi, size_t dim,
1602  size_t i_shift) {
1603 
1604  size_t nval = 1;
1605  for (size_t i = 0; i < dim; ++i) {
1606  if (i == mi) {
1607  if (m[i+i_shift]==0) {
1608  nval*=2;
1609  } else {
1610  nval*= (1 << (m[i+i_shift]));
1611  }
1612  } else {
1613  nval *= (1 << (m[i+i_shift] + 1)) + 1;
1614  }
1615  }
1616  return nval;
1617  }
1618 
1619  /** \brief Desc
1620  */
1621  int add_cacheval(std::vector<cache> &vc,
1622  const std::vector<size_t> &m,
1623  size_t mi, size_t fdim, func_t &f, size_t dim,
1624  const vec_t &xmin, const vec_t &xmax, vec_t &buf,
1625  size_t nbuf) {
1626 
1627  size_t ic = vc.size();
1628  size_t nval, vali = 0, ibuf = 0;
1629  std::vector<double> p(MAXDIM);
1630 
1631  cache c;
1632  vc.push_back(c);
1633 
1634  vc[ic].mi = mi;
1635  for(size_t j=0;j<dim;j++) {
1636  vc[ic].m[j]=m[j];
1637  }
1638  nval = fdim * num_cacheval(m,mi,dim,0);
1639  vc[ic].val.resize(nval);
1640 
1641  if (compute_cacheval(m,mi,vc[ic].val,vali,fdim,f,
1642  dim,0,p,xmin,xmax,buf,nbuf,ibuf)) {
1643  return o2scl::gsl_failure;
1644  }
1645 
1646  if (ibuf > 0) {
1647  /* flush remaining buffer */
1648  const vec_crange_t buf2=o2scl::const_vector_range
1649  (((const vec_t &)buf),0,buf.size());
1650  vec_range_t val2=o2scl::vector_range(vc[ic].val,vali,vc[ic].val.size());
1651  return f(dim, ibuf, &(buf[0]), fdim, &((vc[ic].val)[vali]));
1652  }
1653 
1654  return o2scl::success;
1655  }
1656 
1657  /** \brief Desc
1658 
1659  Recursive loop to evaluate the integral contribution from the
1660  cache entry c, accumulating in val, for the given m[] except
1661  with m[md] -> m[md] - 1 if md < dim, using the cached values
1662  (cm,cmi,cval). id is the current loop dimension (from 0 to
1663  dim-1).
1664  */
1665  size_t eval(const std::vector<size_t> &cm, size_t cmi,
1666  vec_t &cval, const std::vector<size_t> &m, size_t md,
1667  size_t fdim, size_t dim, size_t id,
1668  double weight, vec_t &val, size_t voff2) {
1669 
1670  size_t voff = 0; /* amount caller should offset cval array afterwards */
1671  if (id == dim) {
1672  size_t i;
1673  for (i = 0; i < fdim; ++i) {
1674  val[i] += cval[i+voff2] * weight;
1675  }
1676  voff = fdim;
1677 
1678  } else if (m[id] == 0 && id == md) {
1679 
1680  /* using trivial rule for this dim */
1681  voff = eval(cm, cmi, cval, m, md, fdim, dim, id+1, weight*2, val,voff2);
1682  voff += fdim * (1 << cm[id]) * 2
1683  * num_cacheval(cm, cmi - (id+1), dim - (id+1),id+1);
1684 
1685  } else {
1686 
1687  size_t i;
1688  /* order of C-C rule */
1689  size_t mid = m[id] - (id == md);
1690  const double *w = clencurt_w + mid + (1 << mid) - 1
1691  + (id == cmi ? (cm[id] ? 1 + (1 << (cm[id]-1)) : 1) : 0);
1692  size_t cnx = (id == cmi ? (cm[id] ? (1 << (cm[id]-1)) : 1)
1693  : (1 << (cm[id])));
1694  size_t nx = cm[id] <= mid ? cnx : (1 << mid);
1695 
1696  if (id != cmi) {
1697  voff = eval(cm, cmi, cval, m, md, fdim, dim, id + 1,
1698  weight * w[0], val,voff2);
1699  ++w;
1700  }
1701  for (i = 0; i < nx; ++i) {
1702  voff += eval(cm, cmi, cval, m, md, fdim, dim, id + 1,
1703  weight * w[i], val,voff+voff2);
1704  voff += eval(cm, cmi, cval, m, md, fdim, dim, id + 1,
1705  weight * w[i], val,voff+voff2);
1706  }
1707 
1708  voff += (cnx - nx) * fdim * 2
1709  * num_cacheval(cm, cmi - (id+1), dim - (id+1),id+1);
1710  }
1711  return voff;
1712  }
1713 
1714  /** \brief Desc
1715 
1716  Loop over all cache entries that contribute to the integral,
1717  (with m[md] decremented by 1)
1718  */
1719  void evals(std::vector<cache> &vc, const std::vector<size_t> &m,
1720  size_t md, size_t fdim, size_t dim, double V, vec_t &val) {
1721 
1722  for(size_t k=0;k<fdim;k++) {
1723  val[k]=0.0;
1724  }
1725  for (size_t i = 0; i < vc.size(); ++i) {
1726  if (vc[i].mi >= dim ||
1727  vc[i].m[vc[i].mi] + (vc[i].mi == md) <= m[vc[i].mi]) {
1728  eval(vc[i].m, vc[i].mi, vc[i].val,
1729  m, md, fdim, dim, 0, V, val,0);
1730  }
1731  }
1732  return;
1733  }
1734 
1735  /** \brief Desc
1736 
1737  Evaluate the integrals for the given m[] using the cached
1738  values in vc, storing the integrals in val[], the error
1739  estimate in err[], and the dimension to subdivide next (the
1740  largest error contribution) in *mi
1741  */
1742  void eval_integral(std::vector<cache> &vc,
1743  const std::vector<size_t> &m,
1744  size_t fdim, size_t dim, double V,
1745  size_t &mi, vec_t &val, vec_t &err, vec_t &val1) {
1746 
1747  double maxerr = 0;
1748  size_t i, j;
1749 
1750  evals(vc, m, dim, fdim, dim, V, val);
1751 
1752  /* error estimates along each dimension by comparing val with
1753  lower-order rule in that dimension; overall (conservative)
1754  error estimate from maximum error of lower-order rules. */
1755  for(j=0;j<fdim;j++) {
1756  err[j]=0.0;
1757  }
1758  mi = 0;
1759  for (i = 0; i < dim; ++i) {
1760  double emax = 0;
1761  evals(vc, m, i, fdim, dim, V, val1);
1762  for (j = 0; j < fdim; ++j) {
1763  double e = fabs(val[j] - val1[j]);
1764  if (e > emax) emax = e;
1765  if (e > err[j]) err[j] = e;
1766  }
1767  if (emax > maxerr) {
1768  maxerr = emax;
1769  mi = i;
1770  }
1771  }
1772  return;
1773  }
1774 
1775  /** \brief Desc
1776  */
1777  template<class vec2_t>
1778  int converged(size_t fdim, const vec2_t &vals, const vec2_t &errs,
1779  double reqAbsError, double reqRelError, error_norm norm) {
1780 
1781  switch (norm) {
1782 
1783  case ERROR_INDIVIDUAL:
1784 
1785  for (size_t j = 0; j < fdim; ++j) {
1786  if (errs[j] > reqAbsError && errs[j] > fabs(vals[j])*reqRelError) {
1787  return 0;
1788  }
1789  }
1790  return 1;
1791 
1792  case ERROR_PAIRED:
1793 
1794  {
1795  size_t j;
1796 
1797  for (j = 0; j+1 < fdim; j += 2) {
1798  double maxerr, serr, err, maxval, sval, val;
1799  /* scale to avoid overflow/underflow */
1800  maxerr = errs[j] > errs[j+1] ? errs[j] : errs[j+1];
1801  maxval = vals[j] > vals[j+1] ? vals[j] : vals[j+1];
1802  serr = maxerr > 0 ? 1/maxerr : 1;
1803  sval = maxval > 0 ? 1/maxval : 1;
1804  err=std::hypot(errs[j]*serr,errs[j+1]*serr)*maxerr;
1805  val=std::hypot(vals[j]*sval,vals[j+1]*sval)*maxval;
1806  if (err > reqAbsError && err > val*reqRelError) {
1807  return 0;
1808  }
1809  }
1810  /* fdim is odd, do last dimension individually */
1811  if (j < fdim) {
1812  if (errs[j] > reqAbsError && errs[j] > fabs(vals[j])*reqRelError) {
1813  return 0;
1814  }
1815  }
1816  }
1817  return 1;
1818 
1819  case ERROR_L1:
1820 
1821  {
1822  double err = 0, val = 0;
1823  for (size_t j = 0; j < fdim; ++j) {
1824  err += errs[j];
1825  val += fabs(vals[j]);
1826  }
1827  return err <= reqAbsError || err <= val*reqRelError;
1828  }
1829 
1830  case ERROR_LINF:
1831 
1832  {
1833  double err = 0, val = 0;
1834  for (size_t j = 0; j < fdim; ++j) {
1835  double absval = fabs(vals[j]);
1836  if (errs[j] > err) err = errs[j];
1837  if (absval > val) val = absval;
1838  }
1839  return err <= reqAbsError || err <= val*reqRelError;
1840  }
1841 
1842  case ERROR_L2:
1843 
1844  {
1845  double maxerr = 0, maxval = 0, serr, sval, err = 0, val = 0;
1846  /* scale values by 1/max to avoid overflow/underflow */
1847  for (size_t j = 0; j < fdim; ++j) {
1848  double absval = fabs(vals[j]);
1849  if (errs[j] > maxerr) maxerr = errs[j];
1850  if (absval > maxval) maxval = absval;
1851  }
1852  serr = maxerr > 0 ? 1/maxerr : 1;
1853  sval = maxval > 0 ? 1/maxval : 1;
1854  for (size_t j = 0; j < fdim; ++j) {
1855  err += (errs[j] * serr)*(errs[j] * serr);
1856  val += (fabs(vals[j]) * sval)*(fabs(vals[j]) * sval);
1857  }
1858  err = sqrt(err) * maxerr;
1859  val = sqrt(val) * maxval;
1860  return err <= reqAbsError || err <= val*reqRelError;
1861  }
1862  }
1863 
1864  O2SCL_ERR("Improper value of 'norm' in cubature::converged().",
1866  return 1;
1867  }
1868 
1869  public:
1870 
1871  /** \brief Desc
1872 
1873  Vectorized version with user-supplied buffer to store points
1874  and values. The buffer *buf should be of length *nbuf * dim on
1875  entry (these parameters are changed upon return to the final
1876  buffer and length that was used). The buffer length will be
1877  kept <= max(max_nbuf, 1) * dim.
1878 
1879  Also allows the caller to specify an array m[dim] of starting
1880  degrees for the rule, which upon return will hold the final
1881  degrees. The number of points in each dimension i is 2^(m[i]+1)
1882  + 1.
1883  */
1884  int integ_v_buf(size_t fdim, func_t &f, size_t dim, const vec_t &xmin,
1885  const vec_t &xmax, size_t maxEval, double reqAbsError,
1886  double reqRelError, error_norm norm, std::vector<size_t> &m,
1887  vec_t &buf, size_t &nbuf, size_t max_nbuf, vec_t &val,
1888  vec_t &err) {
1889 
1890  int ret = o2scl::gsl_failure;
1891  double V = 1;
1892  size_t numEval = 0, new_nbuf;
1893  size_t i;
1894  std::vector<cache> vc;
1895 
1896  vec_t val1(fdim);
1897 
1898  /* norm is irrelevant */
1899  if (fdim <= 1) norm = ERROR_INDIVIDUAL;
1900  /* invalid norm */
1901  if (norm < 0 || norm > ERROR_LINF) return o2scl::gsl_failure;
1902  /* nothing to do */
1903  if (fdim == 0) return o2scl::success;
1904  /* unsupported */
1905  if (dim > MAXDIM) return o2scl::gsl_failure;
1906  /* trivial case */
1907  if (dim == 0) {
1908  // AWS: this is one location where vector types need sync'ing
1909  const vec_crange_t xmin2=o2scl::const_vector_range(xmin,0,xmin.size());
1910  vec_range_t val2=o2scl::vector_range(val,0,val.size());
1911  if (f(0, 1, &xmin[0], fdim, &(val[0]))) return o2scl::gsl_failure;
1912  for (i = 0; i < fdim; ++i) err[i] = 0;
1913  return o2scl::success;
1914  }
1915 
1916  for (i = 0; i < fdim; ++i) {
1917  val[i] = 0;
1918  err[i] = HUGE_VAL;
1919  }
1920 
1921  for (i = 0; i < dim; ++i) {
1922  /* scale factor for C-C volume */
1923  V *= (xmax[i] - xmin[i]) * 0.5;
1924  }
1925 
1926  new_nbuf = num_cacheval(m,dim,dim,0);
1927 
1928  if (max_nbuf < 1) max_nbuf = 1;
1929  if (new_nbuf > max_nbuf) new_nbuf = max_nbuf;
1930  if (nbuf < new_nbuf) {
1931  buf.resize((nbuf=new_nbuf)*dim);
1932  }
1933 
1934  /* start by evaluating the m=0 cubature rule */
1935  if (add_cacheval(vc,m, dim, fdim, f, dim, xmin, xmax,
1936  buf, nbuf) != o2scl::success) {
1937  return ret;
1938  }
1939  while (1) {
1940  size_t mi;
1941 
1942  eval_integral(vc,m,fdim,dim,V,mi,val,err,val1);
1943 
1944  if (converged(fdim,val,err,reqAbsError, reqRelError, norm) ||
1945  (numEval > maxEval && maxEval)) {
1946  ret = o2scl::success;
1947  return ret;
1948  }
1949  m[mi] += 1;
1950  if (m[mi] > clencurt_M) {
1951  /* FAILURE */
1952  return ret;
1953  }
1954 
1955  new_nbuf = num_cacheval(m, mi, dim,0);
1956  if (new_nbuf > nbuf && nbuf < max_nbuf) {
1957  nbuf = new_nbuf;
1958  if (nbuf > max_nbuf) nbuf = max_nbuf;
1959  buf.resize(nbuf*dim);
1960  }
1961 
1962  if (add_cacheval(vc,m, mi, fdim, f,
1963  dim, xmin, xmax, buf, nbuf) != o2scl::success) {
1964  /* FAILURE */
1965  return ret;
1966  }
1967  numEval += new_nbuf;
1968  }
1969 
1970  return ret;
1971  }
1972 
1973  /** \brief Desc
1974  */
1975  static const size_t DEFAULT_MAX_NBUF=(1U << 20);
1976 
1977  /** \brief Desc
1978  */
1979  int integ(size_t fdim, func_t &f, size_t dim,
1980  const vec_t &xmin, const vec_t &xmax, size_t maxEval,
1981  double reqAbsError, double reqRelError, error_norm norm,
1982  vec_t &val, vec_t &err) {
1983 
1984  int ret;
1985  size_t nbuf = 0;
1986  std::vector<size_t> m(dim);
1987  vec_t buf;
1988 
1989  /* max_nbuf > 0 to amortize function overhead */
1990  ret = integ_v_buf(fdim,f,dim,xmin,xmax,
1991  maxEval,reqAbsError,reqRelError,norm,
1992  m,buf,nbuf,16,val,err);
1993 
1994  return ret;
1995  }
1996 
1997  };
1998 
1999 }
2000 
2001 #endif
o2scl::inte_hcubature::compute_vol
double compute_vol(const hypercube &h)
Desc.
Definition: cubature.h:208
o2scl::inte_hcubature::ls0
size_t ls0(size_t n)
Functions to loop over points in a hypercube.
Definition: cubature.h:448
o2scl::inte_hcubature::region::fdim
size_t fdim
dimensionality of vector integrand
Definition: cubature.h:303
o2scl::inte_hcubature::errMax
double errMax(const std::vector< esterr > &ee)
Return the maximum error from ee.
Definition: cubature.h:169
o2scl::inte_hcubature::rule::vals_ix
size_t vals_ix
num_regions * num_points * fdim
Definition: cubature.h:380
o2scl::inte_hcubature::rule::num_points
size_t num_points
The number of evaluation points.
Definition: cubature.h:374
o2scl::inte_pcubature::num_cacheval
size_t num_cacheval(const std::vector< size_t > &m, size_t mi, size_t dim, size_t i_shift)
Desc.
Definition: cubature.h:1601
o2scl::inte_hcubature::region::splitDim
size_t splitDim
Desc.
Definition: cubature.h:301
o2scl::inte_hcubature::hypercube::vol
double vol
cache volume = product of widths
Definition: cubature.h:203
boost::numeric::ublas::vector< double >
o2scl::inte_hcubature::rule::fdim
size_t fdim
The number of functions.
Definition: cubature.h:372
o2scl::inte_hcubature::numRR0_0fs
size_t numRR0_0fs(size_t dim)
Desc.
Definition: cubature.h:598
o2scl::inte_hcubature::alloc_rule_pts
void alloc_rule_pts(rule &r, size_t num_regions)
Desc.
Definition: cubature.h:385
o2scl::inte_hcubature::hypercube
Desc.
Definition: cubature.h:179
o2scl::inte_hcubature::rule75genzmalik::weight3
double weight3
Desc.
Definition: cubature.h:652
o2scl::inte_hcubature::heap_item
region heap_item
Desc.
Definition: cubature.h:1009
o2scl::inte_hcubature::esterr
A value and error.
Definition: cubature.h:144
o2scl::inte_hcubature
Adaptive multidimensional integration on hyper-rectangles using cubature rules from the Cubature libr...
Definition: cubature.h:133
o2scl::inte_hcubature::rule::pts
ubvector pts
points to eval: num_regions * num_points * dim
Definition: cubature.h:378
o2scl::inte_hcubature::rule75genzmalik::weightE3
double weightE3
Desc.
Definition: cubature.h:658
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::inte_hcubature::esterr::err
double err
Desc.
Definition: cubature.h:164
o2scl::inte_hcubature::rule75genzmalik::p
std::vector< double > p
Desc.
Definition: cubature.h:647
o2scl::inte_pcubature::add_cacheval
int add_cacheval(std::vector< cache > &vc, const std::vector< size_t > &m, size_t mi, size_t fdim, func_t &f, size_t dim, const vec_t &xmin, const vec_t &xmax, vec_t &buf, size_t nbuf)
Desc.
Definition: cubature.h:1621
o2scl::inte_pcubature
Integration by p-adaptive cubature from the Cubature library.
Definition: cubature.h:1486
o2scl::inte_cubature_base::ERROR_PAIRED
@ ERROR_PAIRED
paired L2 norms of errors in each component, mainly for integrating vectors of complex numbers
Definition: cubature.h:90
o2scl::inte_hcubature::heap_alloc
heap heap_alloc(size_t nalloc, size_t fdim)
Desc.
Definition: cubature.h:1062
o2scl::inte_hcubature::rule75genzmalik::weight5
double weight5
Desc.
Definition: cubature.h:654
o2scl::inte_hcubature::make_hypercube
void make_hypercube(size_t dim, const vec_t &center, const vec_t &halfwidth, hypercube &h)
Desc.
Definition: cubature.h:219
o2scl::inte_hcubature::heap::n
size_t n
Desc.
Definition: cubature.h:1038
o2scl::inte_hcubature::heap::items
std::vector< heap_item > items
Desc.
Definition: cubature.h:1042
o2scl::inte_pcubature::cache
Cache of the values for the m[dim] grid.
Definition: cubature.h:1508
o2scl::inte_hcubature::make_hypercube_range
void make_hypercube_range(size_t dim, const vec_t &xmin, const vec_t &xmax, hypercube &h)
Desc.
Definition: cubature.h:253
o2scl::inte_cubature_base::ERROR_LINF
@ ERROR_LINF
abserr is norm |e|, and relerr is |e|/|v|
Definition: cubature.h:97
o2scl::inte_pcubature::converged
int converged(size_t fdim, const vec2_t &vals, const vec2_t &errs, double reqAbsError, double reqRelError, error_norm norm)
Desc.
Definition: cubature.h:1778
o2scl::inte_hcubature::eval_regions
int eval_regions(size_t nR, std::vector< region > &R, func_t &f, rule &r)
Desc.
Definition: cubature.h:419
o2scl::inte_hcubature::region::h
hypercube h
Desc.
Definition: cubature.h:299
o2scl::inte_hcubature::make_rule
void make_rule(size_t dim, size_t fdim, size_t num_points, rule &r)
Desc.
Definition: cubature.h:405
o2scl::inte_hcubature::heap_pop
heap_item heap_pop(heap &h)
Desc.
Definition: cubature.h:1127
o2scl::inte_hcubature::converged
int converged(size_t fdim, const std::vector< esterr > &ee, double reqAbsError, double reqRelError, error_norm norm)
Desc.
Definition: cubature.h:1176
o2scl::inte_hcubature::heap::fdim
size_t fdim
Desc.
Definition: cubature.h:1044
o2scl::inte_cubature_base::error_norm
error_norm
Different ways of measuring the absolute and relative error.
Definition: cubature.h:85
o2scl::inte_hcubature::region::errmax
double errmax
max ee[k].err
Definition: cubature.h:307
o2scl::inte_hcubature::use_parallel
int use_parallel
Desc.
Definition: cubature.h:1441
o2scl::inte_hcubature::heap::nalloc
size_t nalloc
Desc.
Definition: cubature.h:1040
o2scl::success
@ success
Success.
Definition: err_hnd.h:47
o2scl::inte_hcubature::destroy_hypercube
void destroy_hypercube(hypercube &h)
Desc.
Definition: cubature.h:266
o2scl::inte_hcubature::heap_push
int heap_push(heap &h, heap_item &hi)
Desc.
Definition: cubature.h:1090
o2scl::inte_pcubature::eval
size_t eval(const std::vector< size_t > &cm, size_t cmi, vec_t &cval, const std::vector< size_t > &m, size_t md, size_t fdim, size_t dim, size_t id, double weight, vec_t &val, size_t voff2)
Desc.
Definition: cubature.h:1665
o2scl::inte_pcubature::evals
void evals(std::vector< cache > &vc, const std::vector< size_t > &m, size_t md, size_t fdim, size_t dim, double V, vec_t &val)
Desc.
Definition: cubature.h:1719
o2scl::inte_hcubature::evalRR0_0fs2
void evalRR0_0fs2(ubvector &pts, size_t pts_ix, size_t dim, std::vector< double > &p, size_t p_ix, const std::vector< double > &c, size_t c_ix, const std::vector< double > &r, size_t r_ix)
Desc.
Definition: cubature.h:526
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::inte_hcubature::hypercube::data
std::vector< double > data
length 2*dim = center followed by half-widths
Definition: cubature.h:201
o2scl::inte_pcubature::integ_v_buf
int integ_v_buf(size_t fdim, func_t &f, size_t dim, const vec_t &xmin, const vec_t &xmax, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, std::vector< size_t > &m, vec_t &buf, size_t &nbuf, size_t max_nbuf, vec_t &val, vec_t &err)
Desc.
Definition: cubature.h:1884
o2scl::inte_cubature_base::ERROR_L1
@ ERROR_L1
abserr is L_1 norm |e|, and relerr is |e|/|v|
Definition: cubature.h:94
o2scl::inte_hcubature::rule75genzmalik::weight1
double weight1
dimension-dependent constants
Definition: cubature.h:650
o2scl::inte_hcubature::esterr::val
double val
Desc.
Definition: cubature.h:162
o2scl::inte_hcubature::heap_push_many
int heap_push_many(heap &h, size_t ni, heap_item *hi)
Desc.
Definition: cubature.h:1118
o2scl::inte_hcubature::heap::ee
std::vector< esterr > ee
Definition: cubature.h:1046
o2scl::inte_hcubature::rule15gauss_evalError
int rule15gauss_evalError(rule &r, size_t fdim, func_t &f, size_t nR, std::vector< region > &R)
1d 15-point Gaussian quadrature rule, based on qk15.c and qk.c in GNU GSL (which in turn is based on ...
Definition: cubature.h:848
o2scl::inte_cubature_base
Base class for integration routines from the Cubature library.
Definition: cubature.h:74
o2scl::inte_hcubature::rule75genzmalik_evalError
int rule75genzmalik_evalError(rule &runder, size_t fdim, func_t &f, size_t nR, std::vector< region > &R)
Desc.
Definition: cubature.h:668
o2scl::const_vector_range
const dat_t * const_vector_range(const dat_t *v, size_t start, size_t last)
Vector range function for const pointers.
Definition: vector.h:3231
o2scl::inte_hcubature::rule
Desc.
Definition: cubature.h:340
o2scl::inte_hcubature::make_rule75genzmalik
void make_rule75genzmalik(size_t dim, size_t fdim, rule75genzmalik &r)
Desc.
Definition: cubature.h:811
o2scl::inte_hcubature::integ
int integ(size_t fdim, func_t &f, size_t dim, const vec_t &xmin, const vec_t &xmax, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, vec_t &val, vec_t &err)
Desc.
Definition: cubature.h:1450
o2scl::inte_cubature_base::ERROR_L2
@ ERROR_L2
abserr is L_2 norm |e|, and relerr is |e|/|v|
Definition: cubature.h:92
o2scl::inte_hcubature::evalR0_0fs4d2
void evalR0_0fs4d2(ubvector &pts, size_t pts_ix, size_t dim, std::vector< double > &p, size_t p_ix, const std::vector< double > &c, size_t c_ix, const std::vector< double > &r1, size_t r1_ix, const std::vector< double > &r2, size_t r2_ix)
Desc.
Definition: cubature.h:559
o2scl::inte_hcubature::cubature
int cubature(size_t fdim, func_t &f, size_t dim, const vec_t &xmin, const vec_t &xmax, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, vec_t &val, vec_t &err, int parallel)
Desc.
Definition: cubature.h:1400
o2scl::inte_hcubature::num0_0
size_t num0_0(size_t dim)
Desc.
Definition: cubature.h:592
o2scl::inte_pcubature::cache::val
vec_t val
Desc.
Definition: cubature.h:1531
o2scl::inte_hcubature::rule::dim
size_t dim
The dimensionality.
Definition: cubature.h:370
o2scl::inte_cubature_base::ERROR_INDIVIDUAL
@ ERROR_INDIVIDUAL
individual relerr criteria in each component
Definition: cubature.h:87
o2scl::inte_hcubature::heap_free
void heap_free(heap &h)
Note that heap_free does not deallocate anything referenced by the items.
Definition: cubature.h:1079
o2scl::inte_hcubature::make_region
void make_region(const hypercube &h, size_t fdim, region &R)
Desc.
Definition: cubature.h:312
o2scl::inte_pcubature::integ
int integ(size_t fdim, func_t &f, size_t dim, const vec_t &xmin, const vec_t &xmax, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, vec_t &val, vec_t &err)
Desc.
Definition: cubature.h:1979
o2scl::inte_pcubature::compute_cacheval
int compute_cacheval(const std::vector< size_t > &m, size_t mi, vec_t &val, size_t &vali, size_t fdim, func_t &f, size_t dim, size_t id, std::vector< double > &p, const vec_t &xmin, const vec_t &xmax, vec_t &buf, size_t nbuf, size_t &ibuf)
Desc.
Definition: cubature.h:1540
o2scl::vector_range
dat_t * vector_range(dat_t *v, size_t start, size_t last)
Vector range function for pointers.
Definition: vector.h:3221
o2scl::inte_pcubature::cache::mi
size_t mi
Desc.
Definition: cubature.h:1529
o2scl::inte_hcubature::rulecubature
int rulecubature(rule &r, size_t fdim, func_t &f, const hypercube &h, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, vec_t &val, vec_t &err, int parallel)
Desc.
Definition: cubature.h:1268
o2scl::exc_esanity
@ exc_esanity
sanity check failed - shouldn't happen
Definition: err_hnd.h:65
o2scl::inte_hcubature::heap_resize
void heap_resize(heap &h, size_t nalloc)
Desc.
Definition: cubature.h:1051
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::inte_hcubature::rule75genzmalik::weightE1
double weightE1
Desc.
Definition: cubature.h:656
o2scl::inte_hcubature::make_hypercube2
void make_hypercube2(size_t dim, const std::vector< double > &dat, hypercube &h)
Desc.
Definition: cubature.h:235
o2scl::inte_hcubature::rule75genzmalik
Desc.
Definition: cubature.h:609
o2scl::inte_hcubature::cut_region
int cut_region(region &R, region &R2)
Desc.
Definition: cubature.h:325
o2scl::inte_pcubature::cache::m
std::vector< size_t > m
Desc.
Definition: cubature.h:1527
o2scl::inte_pcubature::eval_integral
void eval_integral(std::vector< cache > &vc, const std::vector< size_t > &m, size_t fdim, size_t dim, double V, size_t &mi, vec_t &val, vec_t &err, vec_t &val1)
Desc.
Definition: cubature.h:1742
o2scl::inte_hcubature::make_rule15gauss
void make_rule15gauss(size_t dim, size_t fdim, rule &r)
Desc.
Definition: cubature.h:992
o2scl::inte_hcubature::hypercube::dim
size_t dim
Desc.
Definition: cubature.h:199
o2scl::gsl_failure
@ gsl_failure
Failure.
Definition: err_hnd.h:49
o2scl::inte_hcubature::region::ee
std::vector< esterr > ee
array of length fdim
Definition: cubature.h:305
o2scl::inte_hcubature::region
Desc.
Definition: cubature.h:274
o2scl::inte_hcubature::rule::num_regions
size_t num_regions
The max number of regions evaluated at once.
Definition: cubature.h:376
o2scl::inte_hcubature::heap
Desc.
Definition: cubature.h:1013
o2scl::inte_hcubature::numR_Rfs
size_t numR_Rfs(size_t dim)
Desc.
Definition: cubature.h:601
o2scl::inte_hcubature::numR0_0fs
size_t numR0_0fs(size_t dim)
Desc.
Definition: cubature.h:595
o2scl::inte_hcubature::evalR_Rfs2
void evalR_Rfs2(ubvector &pts, size_t pts_ix, size_t dim, std::vector< double > &p, size_t p_ix, const std::vector< double > &c, size_t c_ix, const std::vector< double > &r, size_t r_ix)
Evaluate the integration points for all points (+/-r,...+/-r)
Definition: cubature.h:488

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).