From c4ced3eea5a8ca43c8e5fe4b360b712bd02e3503 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (laptop)" <j.wuttke@fz-juelich.de>
Date: Wed, 21 Jan 2009 09:55:32 +0100
Subject: [PATCH] towards conv-integration
---
src/curve.cpp | 71 ++++++--------------------------------------------
src/curve.h | 2 +-
src/expr.cpp | 29 +++++++++++++++++++--
src/frida2.cpp | 2 +-
src/mystd.cpp | 61 +++++++++++++++++++++++++++++++++++++++++++
src/mystd.h | 6 +++++
src/readln.cpp | 4 +--
src/readln.h | 7 ++---
8 files changed, 108 insertions(+), 74 deletions(-)
diff --git a/src/curve.cpp b/src/curve.cpp
index b2112038..7fa99817 100644
--- a/src/curve.cpp
+++ b/src/curve.cpp
@@ -291,6 +291,9 @@ void NCurveFile::SetConvolution( bool onoff )
void NCurveFile::SetConvTuningPars( string which )
{
+ // per-file ???
+ throw string( "convolution tuning parameters not operational" );
+
static int mode = -1;
static double epsabs = 1e-10, epsrel = 1e-10;
COlc *fc;
@@ -617,9 +620,9 @@ void NCurveFile::SetIntTuningPars( string which )
//! Outer wrapper: dialog, loop over files and parameters
-void NCurveFile::Integrate(void)
+void NCurveFile::IntegralProperty(void)
{
- if( NRead::empty() )
+ if( NRead::stack_empty() )
cout << "Generate integral file:\n"
" <list of numbers> extract curve parameters\n"
" * extract all curve parameters\n"
@@ -630,7 +633,7 @@ void NCurveFile::Integrate(void)
if ( mode==" " )
return;
else if ( mode=="i" ){
- if( NRead::empty() )
+ if( NRead::stack_empty() )
cout << "integration bounds: -inf,inf are allowed\n";
numint_bound_low = dask( "Integrate from", numint_bound_low );
numint_bound_hig = dask( "Integrate to", numint_bound_hig );
@@ -732,64 +735,6 @@ void NCurveFile::IntegrateFile( const COlc *fin, const int k,
};
-//! Generic integration routine.
-
-#include <gsl/gsl_integration.h>
-
-double Integrate( double (*func_) (double, void*), void *data_,
- const double low, const double hig,
- const uint mode, const double epsabs, const double epsrel )
-{
- gsl_function F;
- F.function = func_;
- F.params = data_;
-#define nwork 20000
- gsl_integration_workspace *work = gsl_integration_workspace_alloc( nwork );
- double val, err;
- int rule;
- size_t neval;
-
- if ( isinf(low)==-1 && isinf(hig)==+1 )
- gsl_integration_qagi( &F,
- epsabs, epsrel, nwork,
- work, &val, &err );
- else if ( isinf(low)==-1 )
- gsl_integration_qagil( &F, hig,
- epsabs, epsrel, nwork,
- work, &val, &err );
- else if ( isinf(hig)==+1 )
- gsl_integration_qagiu( &F, low,
- epsabs, epsrel, nwork,
- work, &val, &err );
- else if ( mode==0 )
- gsl_integration_qng( &F, low, hig,
- epsabs, epsrel,
- &val, &err, &neval );
- else if ( mode>0 && mode<7 ){
- switch( 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,
- epsabs, epsrel, nwork,
- rule, work, &val, &err );
- }
- else if ( mode==7 )
- gsl_integration_qags( &F, low, hig,
- epsabs, epsrel, nwork,
- work, &val, &err );
- else
- throw string( "invalid mode in NumericIntegral" );
-
- gsl_integration_workspace_free( work );
- return val;
-}
-
-
//! Numeric interation of one curve.
double NCurveFile::NumericIntegral(
@@ -802,8 +747,8 @@ double NCurveFile::NumericIntegral(
data.k = k;
data.j = j;
- return Integrate( &myeval, &data, low, hig,
- numint_mode, numint_epsabs, numint_epsrel );
+ return mystd::Integrate( &myeval, &data, low, hig,
+ numint_mode, numint_epsabs, numint_epsrel );
}
//***************************************************************************//
diff --git a/src/curve.h b/src/curve.h
index 92470fdd..ad026e0c 100644
--- a/src/curve.h
+++ b/src/curve.h
@@ -34,7 +34,7 @@ namespace NCurveFile {
void Fit( void );
void SetIntTuningPars( string which );
- void Integrate( void );
+ void IntegralProperty( void );
void CumulativeDistribution( void );
};
diff --git a/src/expr.cpp b/src/expr.cpp
index 644a716f..07fdc5b8 100644
--- a/src/expr.cpp
+++ b/src/expr.cpp
@@ -12,6 +12,7 @@
#include "mystd.h"
#include "olm.h" // for cross-references
#include "expr.h"
+#include "readln.h"
bool debug = false;
#define DEB(args) if ( debug ) printf args
@@ -339,6 +340,17 @@ string CRef::ref_info(void) const
//* class CTree *//
//***************************************************************************//
+typedef struct {
+ const COlc *fc;
+ uint k;
+ uint j;
+} ConvDatTyp;
+
+double tree_conv_eval( double tt, void *data )
+{
+ ConvDatTyp *mydata = (ConvDatTyp*) data;
+ return 1;
+}
//! Evaluate tree and convert result to integer, for use as index.
@@ -472,7 +484,7 @@ void CTree::tree_val( CTOut *ret, const CContext *ctx ) const
arg[0]->tree_val( &r, ctx );
for ( uint i=0; i<n; ++i )
ret->vd[i] = r.to_d(i);
- } else {
+ } else if ( NRead::iofmac("\\cv_mode")==-1 ) {
// Convolute theory with resolution.
// Simplest version: non-equidistant x_out and x_res.
CTOut r;
@@ -493,7 +505,20 @@ void CTree::tree_val( CTOut *ret, const CContext *ctx ) const
for ( uint i=0; i<n; ++i )
ret->vd[i] += r.to_d(i) * sv->y[iv] / conv_norm;
}
+ } else {
+ ConvDatTyp conv_data;
+ for ( uint i=0; i<n; ++i ) {
+ ret->vd[i] =
+ mystd::Integrate( &tree_conv_eval,
+ (void *) &conv_data,
+ sv->x[0], sv->x[nv-1],
+ NRead::iofmac( "\\cv_mode" ),
+ NRead::dofmac( "\\cv_epsabs" ),
+ NRead::dofmac( "\\cv_epsrel" ) ) /
+ conv_norm;
+ }
}
+
} else if ( typ == _DIRAC ) {
CTOut r;
// center of delta function given by arg[0]:
@@ -752,7 +777,7 @@ void NCalc::Calculator()
PTree T;
double v;
- if ( NRead::empty() ) {
+ if ( NRead::stack_empty() ) {
// interactive calculator
while (1) {
s = sask("calc> ");
diff --git a/src/frida2.cpp b/src/frida2.cpp
index c1b502c1..c99994f9 100644
--- a/src/frida2.cpp
+++ b/src/frida2.cpp
@@ -138,7 +138,7 @@ int main()
NCurveFile::SetConvolution( false );
} else if (cmd == "ci") {
- NCurveFile::Integrate();
+ NCurveFile::IntegralProperty();
} else if (cmd == "cp") {
NCurveFile::ShowPar();
diff --git a/src/mystd.cpp b/src/mystd.cpp
index 055289ea..9b9fdb9c 100644
--- a/src/mystd.cpp
+++ b/src/mystd.cpp
@@ -5,6 +5,7 @@
//* http://frida.sourceforge.net *//
//**************************************************************************//
+#include <stdlib.h> // otherwise buggy compiler error from gsl
#include <string.h>
#include <math.h>
#include <sys/time.h>
@@ -482,3 +483,63 @@ int mystd::freadln(FILE *fd, string *s)
*s += c;
}
}
+
+//**************************************************************************//
+//* GSL integration *//
+//**************************************************************************//
+
+#include <gsl/gsl_integration.h>
+
+double mystd::Integrate( double (*func_) (double, void*), void *data_,
+ const double low, const double hig,
+ const uint mode, const double epsabs,
+ const double epsrel )
+{
+ gsl_function F;
+ F.function = func_;
+ F.params = data_;
+#define nwork 20000
+ gsl_integration_workspace *work = gsl_integration_workspace_alloc( nwork );
+ double val, err;
+ int rule;
+ size_t neval;
+
+ if ( isinf(low)==-1 && isinf(hig)==+1 )
+ gsl_integration_qagi( &F,
+ epsabs, epsrel, nwork,
+ work, &val, &err );
+ else if ( isinf(low)==-1 )
+ gsl_integration_qagil( &F, hig,
+ epsabs, epsrel, nwork,
+ work, &val, &err );
+ else if ( isinf(hig)==+1 )
+ gsl_integration_qagiu( &F, low,
+ epsabs, epsrel, nwork,
+ work, &val, &err );
+ else if ( mode==0 )
+ gsl_integration_qng( &F, low, hig,
+ epsabs, epsrel,
+ &val, &err, &neval );
+ else if ( mode>0 && mode<7 ){
+ switch( 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,
+ epsabs, epsrel, nwork,
+ rule, work, &val, &err );
+ }
+ else if ( mode==7 )
+ gsl_integration_qags( &F, low, hig,
+ epsabs, epsrel, nwork,
+ work, &val, &err );
+ else
+ throw string( "invalid mode in NumericIntegral" );
+
+ gsl_integration_workspace_free( work );
+ return val;
+}
diff --git a/src/mystd.h b/src/mystd.h
index b73e5b0f..14e97244 100644
--- a/src/mystd.h
+++ b/src/mystd.h
@@ -58,4 +58,10 @@ namespace mystd {
const char* mod, string *fname=0);
int freadln(ifstream *f, string *s);
int freadln(FILE *fd, string *s);
+
+ // integration
+ double Integrate( double (*func_) (double, void*), void *data_,
+ const double low, const double hig,
+ const uint mode, const double epsabs,
+ const double epsrel );
}
diff --git a/src/readln.cpp b/src/readln.cpp
index adb4d3f7..94b4bb0e 100644
--- a/src/readln.cpp
+++ b/src/readln.cpp
@@ -48,7 +48,7 @@ namespace NRead {
int iofmac( const string& in );
double dofmac( const string& in );
string wofmac( const string& in );
- bool empty();
+ bool stack_empty();
// internals, used only in this file:
bool initialized=0;
@@ -173,7 +173,7 @@ string NRead::wofmac( const string& in )
//! Anything left on input stack?
-bool NRead::empty()
+bool NRead::stack_empty()
{
return halfline=="";
}
diff --git a/src/readln.h b/src/readln.h
index 5565b210..fddb8b1e 100644
--- a/src/readln.h
+++ b/src/readln.h
@@ -4,14 +4,11 @@ namespace NRead {
/** everything must be declared a second time in readln.cpp **/
string readln(const string&); // all read requests pass here ?
- string readwd(const string&); // all read requests pass here ?
+ string readwd(const string&);
void clear(void);
-// inline string readln(char *prompt) { return readln(string(prompt)); };
-// inline string readwd(char *prompt) { return readwd(string(prompt)); };
string exemac(const string in);
- inline string exemac(const char *in) { return exemac(string(in)); };
int iofmac( const string& in );
double dofmac( const string& in );
string wofmac( const string& in );
- bool empty( void );
+ bool stack_empty( void );
};
--
GitLab