From 9737dd9807ade2761de9f34ad729ebf9e1334b7c Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (office)" <j.wuttke@fz-juelich.de>
Date: Thu, 28 Aug 2008 15:02:07 +0200
Subject: [PATCH] - various integration routines at choice - func' cauchy(x,p)
= p/pi/(x^2+p^2)
---
src/curve.cpp | 93 +++++++++++++++++++++++++++++++++++++++++---------
src/expr.cpp | 2 +-
src/func.cpp | 3 ++
src/func.h | 3 +-
src/readln.cpp | 6 ++--
src/readln.h | 2 +-
6 files changed, 86 insertions(+), 23 deletions(-)
diff --git a/src/curve.cpp b/src/curve.cpp
index ee3a1de6..cf256604 100644
--- a/src/curve.cpp
+++ b/src/curve.cpp
@@ -13,6 +13,7 @@
#include <algorithm>
#include "olm.h" // includes also curve.h
+#include "readln.h"
#include "asi.h"
#include "../../lmfit/lmmin.h"
#include "xax.h"
@@ -465,8 +466,9 @@ namespace NCurveFile
double NumericIntegral( const COlc *fc, const CCurve *c,
const int k, const int j );
// parameters for numeric integration:
- double numint_bound_low=0, numint_bound_hig=0;
+ string numint_bound_low="!", numint_bound_hig="!";
// tuning parameters:
+ int numint_mode=4;
double numint_epsabs=1e-10, numint_epsrel=1e-10;
};
@@ -476,44 +478,56 @@ namespace NCurveFile
void NCurveFile::SetIntTuningPars( string which )
{
if ( which=="?" )
- printf( "Numeric integration setup:\n epsabs=%g\n epsrel=%g\n",
- numint_epsabs, numint_epsrel );
+ printf( "Numeric integration setup:\n epsabs=%g\n epsrel=%g\n"
+ " mode: %d\n",
+ numint_epsabs, numint_epsrel, numint_mode );
else if ( which=="a" )
numint_epsabs = iask( "Espilon abs", numint_epsabs );
else if ( which=="r" )
numint_epsrel = iask( "Epsilon rel", numint_epsrel );
+ else if ( which=="m" )
+ cout << "Numeric integration modes (ignored for infinite intervals):\n"
+ " 0: non-adaptive Gauss-Kronod-87 for very smooth functions\n"
+ " 1: adaptive Gauss-Kronod-61 for smooth functions\n"
+ " ...\n"
+ " 6: adaptive Gauss-Kronod-15 for not so smooth functions\n"
+ " 7: adaptive for functions with singularities\n";
+ numint_mode = iask( "Option", numint_mode );
}
//! Outer wrapper: dialog, loop over files and parameters
void NCurveFile::Integrate(void)
{
- cout << "Generate integral file:\n"
- " <list of numbers> extract curve parameters\n"
- " * extract all curve parameters\n"
- " q0 | q1 | q2 extract fit metrics\n"
- " i numeric integral\n";
+ if( NRead::empty() )
+ cout << "Generate integral file:\n"
+ " <list of numbers> extract curve parameters\n"
+ " * extract all curve parameters\n"
+ " q0 | q1 | q2 extract fit metrics\n"
+ " i numeric integral\n";
static string mode = "*";
mode = wask( "Option", mode );
if ( mode==" " )
return;
else if ( mode=="i" ){
- numint_bound_low = dask( "integrate from", numint_bound_low );
- numint_bound_hig = dask( "integrate to", numint_bound_hig );
+ if( NRead::empty() )
+ cout << "integration bounds: use ! for +/- infinity\n";
+ numint_bound_low = wask( "Integrate from", numint_bound_low );
+ numint_bound_hig = wask( "Integrate to", numint_bound_hig );
}
NOlm::IterateC fiter;
const COlc* fin;
while( (fin=fiter()) ){
- if( isdigit(mode[0]) || mode=="*" ){
+ if( isdigit(mode[0]) || mode=="*" ){ // extract parameter(s)
// decode list of parameter numbers
CList parsel( mode, CLimits( 0, fin->nPar()-1 ) );
for( uint iv=0; iv<parsel.n(); ++iv ){
int iP = parsel.V[iv];
IntegrateFile( fin, fiter.SelNo(), "p", iP, "p"+strg(iP) );
}
- } else
+ } else // other integrals
IntegrateFile( fin, fiter.SelNo(), mode, -1, mode );
}
}
@@ -559,6 +573,9 @@ void NCurveFile::IntegrateFile( const COlc *fin, const int k,
val = fin->VC[j].Quality[2];
else if ( mode=="i" )
val = NumericIntegral( fin, &(fin->VC[j]), k, j );
+ else
+ throw string("invalid mode");
+
// print out ?
if (!savable) {
cout << fin->name << ": " << name << " = " << val << "\n";
@@ -601,9 +618,51 @@ double NCurveFile::NumericIntegral( const COlc *fc, const CCurve *c,
F.params = &data;
gsl_integration_workspace *work = gsl_integration_workspace_alloc( nwork );
double val, err;
-
- gsl_integration_qag( &F, numint_bound_low, numint_bound_hig,
- numint_epsabs, numint_epsrel, nwork,
- GSL_INTEG_GAUSS31, work, &val, &err );
- cout << "val="<<val<<", err="<<err<<"\n";
+ int rule;
+ size_t neval;
+
+ double low, hig;
+
+ if( numint_bound_low!="!" && !mystd::any2dbl( numint_bound_low, &low ) )
+ throw string( "invalid lower integration bound" );
+ if( numint_bound_hig!="!" && !mystd::any2dbl( numint_bound_hig, &hig ) )
+ throw string( "invalid upper integration bound" );
+
+ if ( numint_bound_low=="!" && numint_bound_hig=="!" )
+ gsl_integration_qagi( &F,
+ numint_epsabs, numint_epsrel, nwork,
+ work, &val, &err );
+ else if ( numint_bound_low=="!" )
+ gsl_integration_qagil( &F, hig,
+ numint_epsabs, numint_epsrel, nwork,
+ work, &val, &err );
+ else if ( numint_bound_hig=="!" )
+ gsl_integration_qagiu( &F, low,
+ numint_epsabs, numint_epsrel, nwork,
+ work, &val, &err );
+ else if ( numint_mode==0 )
+ gsl_integration_qng( &F, low, hig,
+ numint_epsabs, numint_epsrel,
+ &val, &err, &neval );
+ else if ( numint_mode>0 && numint_mode<7 ){
+ switch( numint_mode ){
+ case 1: rule = GSL_INTEG_GAUSS61; break;
+ case 2: rule = GSL_INTEG_GAUSS51; break;
+ case 3: rule = GSL_INTEG_GAUSS41; break;
+ case 4: rule = GSL_INTEG_GAUSS31; break;
+ case 5: rule = GSL_INTEG_GAUSS21; break;
+ case 6: rule = GSL_INTEG_GAUSS15; break;
+ }
+ gsl_integration_qag( &F, low, hig,
+ numint_epsabs, numint_epsrel, nwork,
+ rule, work, &val, &err );
+ }
+ else if ( numint_mode==7 )
+ gsl_integration_qags( &F, low, hig,
+ numint_epsabs, numint_epsrel, nwork,
+ work, &val, &err );
+ else
+ throw string( "invalid case in NumericIntegral" );
+
+ return val;
}
diff --git a/src/expr.cpp b/src/expr.cpp
index abcdbc58..3289d58e 100644
--- a/src/expr.cpp
+++ b/src/expr.cpp
@@ -763,7 +763,7 @@ void NCalc::Calculator()
// if there is an expression on the stack, then
// evaluate just this single expression;
// otherwise start a loop and request input
- bool expr_on_stack = NRead::nonempty();
+ bool expr_on_stack = !NRead::empty();
while (1) {
if( expr_on_stack ){
diff --git a/src/func.cpp b/src/func.cpp
index 494ee6d9..949d7b79 100644
--- a/src/func.cpp
+++ b/src/func.cpp
@@ -106,6 +106,8 @@ double func_min (double v, double a) { return v<a ? v : a; }
double func_max (double v, double a) { return v>a ? v : a; }
double func_gauss (double x, double s) { return s==0 ? 0 :
exp(-x*x/(2*s*s))/sqrt(twopi)/s; }
+double func_cauchy (double x, double s) { return s==0 ? 0 :
+ 2*s/twopi/(x*x+s*s); }
double func_ran (double v, double a) {
// a random number between v and a :
@@ -207,6 +209,7 @@ void NFunctions::Register(void) {
fmap.insert(make_pair(string("min"), CFunc("min", func_min)));
fmap.insert(make_pair(string("max"), CFunc("max", func_max)));
fmap.insert(make_pair(string("gauss"),CFunc("gauss", func_gauss)));
+ fmap.insert(make_pair(string("cauchy"),CFunc("cauchy", func_cauchy)));
fmap.insert(make_pair(string("ran"), CFunc("ran", func_ran)));
diff --git a/src/func.h b/src/func.h
index 369dcea6..eddf3de2 100644
--- a/src/func.h
+++ b/src/func.h
@@ -108,9 +108,10 @@ extern double func_in( double, vector<double> );
extern double func_min(double,double);
extern double func_max(double,double);
extern double func_gauss(double,double);
+extern double func_cauchy(double,double);
extern double func_ran(double,double);
-extern double func_condass(double,double,double);
+extern double func_condass(double,double,double); // conditional assignement
extern double func_q4w(double,double,double);
diff --git a/src/readln.cpp b/src/readln.cpp
index e5741546..13645966 100644
--- a/src/readln.cpp
+++ b/src/readln.cpp
@@ -48,7 +48,7 @@ namespace NRead {
int iofmac( const string in, int *out );
int dofmac( const string in, double *out );
int wofmac( const string in, string *out );
- bool nonempty();
+ bool empty();
// internals, used only in this file:
bool initialized=0;
@@ -170,9 +170,9 @@ int NRead::wofmac( const string in, string *out )
//! Anything left on input stack?
-bool NRead::nonempty()
+bool NRead::empty()
{
- return halfline!="";
+ return halfline=="";
}
diff --git a/src/readln.h b/src/readln.h
index e3c3166c..5408e4fc 100644
--- a/src/readln.h
+++ b/src/readln.h
@@ -13,5 +13,5 @@ namespace NRead {
int iofmac(const string in, int *out);
int dofmac(const string in, double *out);
int wofmac(const string in, string *out);
- bool nonempty( void );
+ bool empty( void );
};
--
GitLab