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

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