Skip to content

GitLab

Commit 9ea9156e authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC) 7/10/15 add quadrature integration features

parent d620dee7
Show whitespace changes
Inline Side-by-side
Showing with 1975 additions and 0 deletions
data/ClenshawCurtisRuleData-20.txt 0 → 100644
View file @ 9ea9156e
0.000000000000000 0.000346500346500 -0.001038541204746
0.001707753496665 0.003332307267386 0.003332307267386
0.006819348298639 0.006832484617804 -0.006453160953000
0.015299867030335 0.010133277051537 0.010133277051537
0.027091379149683 0.013430557749261 -0.013524623295350
0.042113336672471 0.016599075431932 0.016599075431932
0.060263124396756 0.019678233781133 -0.019616706437387
0.081416760868736 0.022605913597255 0.022605913597255
0.105429745301803 0.025392060329845 -0.025417036636148
0.132138044663434 0.027994547472752 0.027994547472752
0.161359214187130 0.030414189524314 -0.030392340702803
0.192893643655166 0.032618969763790 0.032618969763790
0.226525920938786 0.034607099581875 -0.034615949571773
0.262026303481463 0.036353347145552 0.036353347145552
0.299152287673515 0.037856229247877 -0.037846795338222
0.337650265397658 0.039095936840950 0.039095936840950
0.377257256429601 0.040072870771074 -0.040074749632227
0.417702704859633 0.040771981808914 0.040771981808914
0.458710327263834 0.041196523877068 -0.041193346401595
0.500000000000000 0.041335787586367 0.041335787586367
0.541289672736166 0.041196523877068 -0.041193346401595
0.582297295140367 0.040771981808914 0.040771981808914
0.622742743570400 0.040072870771074 -0.040074749632227
0.662349734602342 0.039095936840950 0.039095936840950
0.700847712326485 0.037856229247877 -0.037846795338222
0.737973696518537 0.036353347145552 0.036353347145552
0.773474079061213 0.034607099581875 -0.034615949571773
0.807106356344834 0.032618969763790 0.032618969763790
0.838640785812871 0.030414189524314 -0.030392340702803
0.867861955336566 0.027994547472752 0.027994547472752
0.894570254698197 0.025392060329845 -0.025417036636148
0.918583239131264 0.022605913597255 0.022605913597255
0.939736875603244 0.019678233781133 -0.019616706437387
0.957886663327529 0.016599075431932 0.016599075431932
0.972908620850317 0.013430557749261 -0.013524623295350
0.984700132969665 0.010133277051537 0.010133277051537
0.993180651701361 0.006832484617804 -0.006453160953000
0.998292246503335 0.003332307267386 0.003332307267386
1.000000000000000 0.000346500346500 -0.001038541204746
src/quadinteg.h 0 → 100644
View file @ 9ea9156e
#ifndef QUADINTEG_H_SEEN
#define QUADINTEG_H_SEEN
// include standard c/c++
#include <cmath>
#include <cfloat>
#include <stdlib.h>
#include <vector>
#include <list>
#include <utility>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <limits> //for precision features
#include "walltimer.h" //timing
#define NAMESPACETAG LagSHT
using namespace std;
namespace NAMESPACETAG {
/* 1D function Class (From Sophya)
*/
/*
\class GenericFunction
\brief Base class definition for objects representing functions
*/
class GenericFunction {
public:
GenericFunction() { }
virtual ~GenericFunction() { }
};
/*!
\class ClassFunc1D
\brief Abstract interface definition for 1 D real functions viewed as classes
This class represents real function of a single real argument : double f(double x)
Inheriting classes should define the double operator()(double x)
For convenience, ClassFunc has been defined using a typedef as an alias
*/
class ClassFunc1D : public GenericFunction {
public:
virtual double operator()(double x) const =0;
};
//! typedef definition of ClassFunc as ClassFunc1D for convenience
typedef ClassFunc1D ClassFunc ;
/*!
\class Function1D
\brief ClassFunc constructed from a a function:simple forward to function
For convenience, ClassFunc has been defined using a typedef as an alias
*/
//ClassFunc constructed from a a function:simple forward to function
class Function1D : public ClassFunc1D {
private:
double (*f)(double);
public:
Function1D(double (*g)(double)):f(g){}
virtual double operator()(double x) const { return f(x);}
};
typedef Function1D Function;
/*
\class Quadrature
\brief quadrature defined in [0,1]
*/
template <typename T, typename TDerived>
class Quadrature {
public:
typedef pair<T,T> values_t; //{intgral value, error}
typedef pair<T,T> bounds_t; //{a,b} bounds
typedef pair< bounds_t , values_t > result_t; // { {a,b},{int,err} }
typedef list< result_t > colresults_t; // collection of result_t
/*
Utility Sort funtion according to error (largest first)
*/
struct Comparator {
Quadrature* m;
Comparator(Quadrature* p): m(p) {}
bool operator()(const result_t& r1, const result_t& r2) {
return (r1.second).second > (r2.second).second;
}
};
public:
/*
Default Ctor
*/
Quadrature(string name="",size_t npts=0,string fName="", bool init=false):
name_(name),
npts_(npts),
func_(0){
if(init){
//compute abscissa, weights, error weights and save into fName
cout << "Init for " << name << " with " << npts << " pts" << endl;
TDerived& tDerivedObj = (TDerived&)*this;
try {
tDerivedObj.ComputeAbsWeights(fName);
if(fName.length()!=0) {//save
ofstream fout(fName.c_str(), ofstream::out);
int size = absc_.size();
for (int i=0; i<size;i++){
fout <<fixed << showpoint << setprecision(15) << absc_[i] << "\t" << absw_[i] << "\t" << errw_[i] << endl;
}
fout.close();
}//save
} catch (::exception& ex) {
string msg = "Init ";
msg += name + ": ERROR ";
msg += ex.what();
::cerr << msg << endl;
exit(EXIT_FAILURE);
} catch (...) {
::cerr << "Init " << name << " unexpected error " << endl;
exit(EXIT_FAILURE);
}
} else {
//read abscissa, weights, error weights from fName
ifstream inFile;
inFile.exceptions ( ifstream::eofbit | ifstream::failbit | ifstream::badbit );
absc_.resize(npts);
absw_.resize(npts);
errw_.resize(npts);
try {
inFile.open(fName.c_str());
for(size_t i=0;i<npts;i++){
inFile >> absc_[i]>> absw_[i] >> errw_[i];
}
inFile.close();
} catch (::ifstream::failure e) {
::cerr << name <<" Exception opening/reading/closing file <"<<fName<<">"<<::endl;
exit(EXIT_FAILURE);
} //try/catch read
}//initialisation
}//Ctor
/*
Dtor: care cannot be virtual because of template
*/
~Quadrature() {}
/*
Rescale abscissa and weights from [xInmin, xInmax] to [xOutmin, xOutmax]
Typical use-case the abs & weights are computed for (-1,1) range
and we need to translate them to (0,1) or vice-versa
*/
void RescaleAbsWeights(T xInmin = -1.0, T xInmax = 1.0,
T xOutmin=0.0, T xOutmax=1.0) {
T deltaXIn = xInmax-xInmin;
T deltaXOut= xOutmax-xOutmin;
T scale = deltaXOut/deltaXIn;
for(size_t i=0; i<npts_;i++){
absc_[i] = ((absc_[i]-xInmin)*xOutmax-(absc_[i]-xInmax)*xOutmin)/deltaXIn;
absw_[i] *= scale;
errw_[i] *= scale;
}
}//RescaleAbsWeights
/*
Set function
*/
void SetFunction(ClassFunc1D& aFunc) { func_ = &aFunc;}
/*
Set the integration bounds
*/
void SetBounds(T xmin, T xmax) { initBounds_ = make_pair(xmin,xmax); }
/*
Set both Function & bounds
*/
void SetFuncBounded(ClassFunc1D& aFunc, T xmin, T xmax) {
SetFunction(aFunc);
SetBounds(xmin,xmax);
}
/*
Get abcissa, abscissa weights, error weights
*/
vector<T> GetAbscissa() const {return absc_;}
vector<T> GetAbscissaW() const {return absw_;}
vector<T> GetErrorW() const {return errw_;}
/*
Get number of points used
*/
size_t GetNPoints() const {return npts_;}
size_t GetOrder() const {return 0;} //should be reimplemented
/*
Debug
*/
void ToOut(::ostream& out, const result_t& aResult) {
out << "{" << aResult.first.first << ",,"
<< aResult.first.second << "}, "
<< aResult.second.first << ", "
<< aResult.second.second << "}";
}//ToOut
void ToOut(::ostream& out, const colresults_t& aColRes) {
typename colresults_t::const_iterator iRes,iEnd;
iEnd = aColRes.end();
out << "{";
for(iRes = aColRes.begin();iRes!=iEnd;++iRes){
ToOut(out,*iRes); cout << ",";
}
out << "}" << endl;
}//ToOut
/*
Compute func integral between the current bounds
The abscissa/weights of the Quadratures are defined for (0,1) interval.
*/
values_t IntEval(const bounds_t& curBounds) {
T integral=0.;
T error = 0.;
T boundsDist = curBounds.second - curBounds.first;
for (size_t i=0;i<npts_;i++){
T xi = curBounds.first + absc_[i]*boundsDist;
T fi = (*func_)(xi);
integral += absw_[i] * fi;
error += errw_[i] * fi;
}
return make_pair(integral*boundsDist,fabs((double)error*boundsDist));
}//IntEval
/* Recursive computation of integral on a finite range
\input curBounds = {x0,x1}
\input cut: threshold to be compared to the integral error to stop or not
\output : {integral, error}
Todo: 1) improve the partitionning (eg. Mathematica uses the quadrature abscissa rescaled)
2) revisit the condition
*/
values_t LocalAdapStratRecur(const bounds_t& curBounds, T cut){
values_t curValues = IntEval(curBounds); //{int,err}
//Not used T integral = curValues.first;
T error = curValues.second;
if (error <= cut) return curValues;
//split the initial range in two partition and decrease by 2 the cut on each part
T xmiddle = (curBounds.first+curBounds.second)*0.5; //xmid
bounds_t lowerB = make_pair(curBounds.first,xmiddle); //{x0,xmid}
values_t valLow = LocalAdapStratRecur(lowerB,cut/2.0);
bounds_t upperB = make_pair(xmiddle,curBounds.second); //{xmid,x1}
values_t valUp = LocalAdapStratRecur(upperB,cut/2.0);
//total integral= sum of sub integral, idem for the error
return make_pair(valLow.first+valUp.first,valLow.second+valUp.second);
}//LocalAdapStratRecur
/*
Local Adaptative Strategy
*/
values_t LocalAdapStrat(T tol = (T)1.0e-6,
size_t MaxNumRecurssion = 1000, // not used
size_t MinNumRecurssion = 5 // not used
) {
bounds_t curBounds = initBounds_; // {a,b}
values_t curValues = IntEval(curBounds); //{int,err}
T integral = curValues.first;
T error = curValues.second;
result_t curResult = make_pair(curBounds,curValues); //{{a,b}, {int,err}}
T cut = tol*fabs((double)integral);
if (error <= cut) return curValues;
return LocalAdapStratRecur(curBounds,cut);
}//LocalAdapStrat
/*
Global Adaptative Strategy of integration to obtain error<tol*integral
add a maximum number of iterations. 1000 is large enough and so may
signal a numerical problem.
*/
values_t GlobalAdapStrat(T tol = (T)1.0e-6,
size_t MaxNumRecurssion = 1000,
size_t MinNumRecurssion = 5) {
bounds_t curBounds = initBounds_; // {a,b}
values_t curValues = IntEval(curBounds); //{int,err}
T integral = curValues.first;
T error = curValues.second;
result_t curResult = make_pair(curBounds,curValues); //{{a,b}, {int,err}}
colresults_t theResults;
typename colresults_t::const_iterator iRes, iEnd;
theResults.push_back(curResult);
size_t n=0;
while( ((error >= tol*fabs(integral)) && (n<MaxNumRecurssion)) ||
(n<MinNumRecurssion) ) {
n++;
curBounds = theResults.front().first; //{x0,x2}
T xmiddle = (curBounds.first+curBounds.second)*0.5; //x1
bounds_t lowerB = make_pair(curBounds.first,xmiddle); //{x0,x1}
values_t valLow = IntEval(lowerB);
result_t resLow = make_pair(lowerB,valLow);
bounds_t upperB = make_pair(xmiddle,curBounds.second); //{x1,x2}
values_t valUp = IntEval(upperB);
result_t resUp = make_pair(upperB,valUp);
//replace the first element (the one which has the highest error)
//of list of results by the two new results (resLow,resUp
theResults.pop_front();
theResults.push_back(resUp);
theResults.push_back(resLow);
//sort the list of results
Comparator compare(this);
theResults.sort(compare);
integral = 0.;
error = 0.;
iEnd=theResults.end();
for(iRes=theResults.begin();iRes!=iEnd;++iRes){
result_t iResCur=*iRes;
integral += (iResCur.second).first;
error += (iResCur.second).second;
}
}//eod while
//debug
/* cout << "GlobalAdapStrat end: {" << n << ", " << integral << ", " */
/* << error << "}" << endl; */
return make_pair(integral,error);
}//GlobalAdapStrat
protected:
string name_;
size_t npts_;
vector<T> absc_;
vector<T> absw_;
vector<T> errw_;
ClassFunc1D* func_;
bounds_t initBounds_;
private:
/*
Explicit prohibit Copy Ctor for the time beeing
*/
Quadrature(const Quadrature&);
/*
Explicit prohibit Assignment Operator for the time beeing
*/
Quadrature& operator=(const Quadrature&);
};//end Quadrature
//-------------------------------------------------------------------------------------
template <typename T>
class ClenshawCurtisQuadrature : public Quadrature<T, ClenshawCurtisQuadrature<T> > {
public:
/*
Default Ctor
*/
ClenshawCurtisQuadrature():
Quadrature<T,ClenshawCurtisQuadrature<T> >("ClenshawCurtisQuadrature",39,"ClenshawCurtisRuleData-20.txt",false) {
}//Ctor
/*
Ctor used with a different number of points = 2n-1. Should implement ComputeAbsWeights
*/
ClenshawCurtisQuadrature(size_t n, string fName, bool init):
Quadrature<T,ClenshawCurtisQuadrature<T> >("ClenshawCurtisQuadrature",2*n-1,fName,init) {}
size_t GetOrder() const {return 2*this->npts_-1;}
/*
Implement the abscissa, weights and err weights in the range (0,1)
*/
void ComputeAbsWeights(string fName) throw(string) {
//code adapted from John Burkart 2009
//the abs, weights are for (-1,1) range so rescaling (0,1) is necessary
//the error weights are adapted from Mathematica prescription
//Todo: contact him to see how to use his code?
int order = this->npts_;
if (order < 3){
string msg = this->name_ + " ComputeAbsWeights: Too few points ERROR";
throw(msg);
}
this->absc_.resize(order);
this->absw_.resize(order);
this->errw_.resize(order);
//Compute the abcsissa locations and their weights
ComputeAbsWeights(order,this->absc_,this->absw_);
//Compute the error weights based on the following properties:
//the abscissa of computed with n points is a subset of the abscissa computed with 2n-1 points
// so the difference between the two integral approximations give an error estimate.
int subOrder = (order+1)/2;
vector<T> xSub(subOrder);
vector<T> wSub(subOrder);
ComputeAbsWeights(subOrder,xSub,wSub);
for(int i=0; i<order;i++){
this->errw_[i] = (i%2 == 0) ? this->absw_[i] - wSub[i/2] : this->absw_[i];
}
//Explicit rescaling (-1,1) -> (0,1)
this->RescaleAbsWeights();
}//ComputeAbsWeights
/*
ComputeAbsWeights more related to the original code
*/
void ComputeAbsWeights(int order, vector<T>& x, vector<T>& w) throw(string) {
double b;
int i;
int j;
// double pi = 3.141592653589793;
double pi = M_PI; //JEC use cmath definition of PI
double theta;
//Todo this kind of error should be tracked before
if ( order < 1 ) {
cerr << "\n";
cerr << "CLENSHAW_CURTIS_COMPUTE - Fatal error!\n";
cerr << " Illegal value of ORDER = " << order << "\n";
exit (EXIT_FAILURE);
}
if ( order == 1 ) {
x[0] = 0.0;
w[0] = 2.0;
} else {
for ( i = 0; i < order; i++ ) {
x[i] = cos ( ( double ) ( order - 1 - i ) * pi
/ ( double ) ( order - 1 ) );
}
x[0] = -1.0;
if ( ( order % 2 ) == 1 ) {
x[(order-1)/2] = 0.0;
}
x[order-1] = +1.0;
for ( i = 0; i < order; i++ ) {
theta = ( double ) ( i ) * pi / ( double ) ( order - 1 );
w[i] = 1.0;
for ( j = 1; j <= ( order - 1 ) / 2; j++ ) {
if ( 2 * j == ( order - 1 ) ) {
b = 1.0;
} else {
b = 2.0;
}
w[i] = w[i] - b * cos ( 2.0 * ( double ) ( j ) * theta )
/ ( double ) ( 4 * j * j - 1 );
}
}
w[0] = w[0] / ( double ) ( order - 1 );
for ( i = 1; i < order - 1; i++ ) {
w[i] = 2.0 * w[i] / ( double ) ( order - 1 );
}
w[order-1] = w[order-1] / ( double ) ( order - 1 );
}
}//ComputeAbsWeights
/*
Destructor
*/
~ClenshawCurtisQuadrature() {}
}; //ClenshawCurtisQuadrature
//-------------------------------------------------------------------------------------
template <typename T>
class GaussKronrodQuadrature : public Quadrature<T, GaussKronrodQuadrature<T> > {
public:
GaussKronrodQuadrature():
Quadrature<T,GaussKronrodQuadrature<T> >("GaussKronrodQuadrature",41,"GaussKronrodRuleData-20.txt",false) {
}//Ctor
/*
Ctor used with a different number of points = 2n+1. Should implement ComputeAbsWeights
*/
GaussKronrodQuadrature(size_t n, string fName, bool init):
Quadrature<T, GaussKronrodQuadrature<T> >("GaussKronrodQuadrature",2*n+1,fName,init) {}
size_t GetOrder() const {return 2*this->npts_+1;}
/*
Implement the abscissa, weights and err weights in the range (0,1)
*/
void ComputeAbsWeights(string fName) throw(string) {
//code adapted from John Burkart 19/3/2009
//the abs, weights are for (-1,1) range so rescaling (0,1) is necessary
//the error weights are adapted from Mathematica prescription
//Todo: contact him to see how to use his code?
int npts = this->npts_;
this->absc_.resize(npts);
this->absw_.resize(npts);
this->errw_.resize(npts);
int order= (npts-1)/2;
if (order < 3){
string msg = this->name_ + " ComputeAbsWeights: Too few points ERROR";
throw(msg);
}
double* x = new double[order+1];
double* w1 = new double[order+1];
double* w2 = new double[order+1];
double eps = 0.000001;
kronrod(order, eps, x , w1, w2);
for (int i=0;i<(order+1);i++){
// cout << "["<<i<<"]"<<fixed << showpoint << setprecision(6) << x[i] << "\t" << w1[i] << "\t" << w2[i] << endl;
this->absc_[i] = -x[i];
this->absw_[i] = w1[i];
this->errw_[i] = w1[i]-w2[i];
int j=npts-1-i;
this->absc_[j] = x[i];
this->absw_[j] = this->absw_[i];
this->errw_[j] = this->errw_[i];
}
delete [] w2;
delete [] w1;
delete [] x;
//Explicit rescaling (-1,1) -> (0,1)
this->RescaleAbsWeights();
}//ComputeAbsWeights