From fbb6c62b7434f86e309313709f6b385f6e6c0178 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (laptop)" <j.wuttke@fz-juelich.de>
Date: Thu, 6 May 2010 23:22:09 +0200
Subject: [PATCH] important clarification: separate tree_val according to typ
---
pub/src/expr.cpp | 660 ++++++++++++++++++++++++-----------------------
pub/src/expr.h | 17 +-
pub/src/opr.cpp | 4 +-
3 files changed, 357 insertions(+), 324 deletions(-)
diff --git a/pub/src/expr.cpp b/pub/src/expr.cpp
index ce5cd369..a3219963 100644
--- a/pub/src/expr.cpp
+++ b/pub/src/expr.cpp
@@ -216,8 +216,7 @@ uint CRef::get_k( const CContext& ctx ) const
else
k = ctx.k;
if ( k==(uint)-1 )
- throw string( "k-dependent expression out of context or "
- "k out of bounds" );
+ throw string( "missing or invalid k" );
if ( k>=NOlm::MOM.size())
throw string( "index k=" + strg(k) + " out of bounds" );
@@ -236,21 +235,20 @@ uint CRef::get_j( const CContext& ctx ) const
j = ctx.j;
}
if ( j==(uint)-1 )
- throw string( "j-dependent expression out of context or "
- "j out of bounds" );
+ throw string( "missing or invalid j" );
return j;
}
//! Evaluate this reference in given context. Result returned in ret.
-void CRef::ref_val( CTOut *ret, const CContext& ctx ) const
+void CRef::ref_val( CTOut& ret, const CContext& ctx ) const
{
uint k = get_k( ctx );
// Return k-dependent reference:
if ( typ == _K ) {
- ret->set_u( k );
+ ret.set_u( k );
return;
}
@@ -262,15 +260,15 @@ void CRef::ref_val( CTOut *ret, const CContext& ctx ) const
// Return k-j-dependent reference:
if ( typ == _J ) {
- ret->set_u( ctx.j );
+ ret.set_u( ctx.j );
return;
} else if ( typ==_NJ ) {
- ret->set_u( f->nJ() );
+ ret.set_u( f->nJ() );
return;
} else if ( typ == _Z ) {
if ( num>= f->nZ() )
throw ( "invalid reference " + var_info() );
- ret->set_d( f->V[j]->z[num] );
+ ret.set_d( f->V[j]->z[num] );
return;
}
@@ -293,25 +291,25 @@ void CRef::ref_val( CTOut *ret, const CContext& ctx ) const
if( i>= sj->size() )
throw string( "i out of range" );
if ( typ == _I ) {
- ret->set_u( i );
+ ret.set_u( i );
} else if ( typ == _X ) {
- ret->set_d( sj->x[i] );
+ ret.set_d( sj->x[i] );
} else if ( typ == _Y ) {
- ret->set_d( sj->y[i] );
+ ret.set_d( sj->y[i] );
if ( ctx.want_error && sj->dy.size() )
- ret->dd = sj->dy[i];
+ ret.dd = sj->dy[i];
} else if ( typ == _DY ) {
- ret->set_d( sj->dy[i] );
+ ret.set_d( sj->dy[i] );
}
} else if ( ctx.dim_vi() ) {
- ret->preset_vd( ctx.nv );
+ ret.preset_vd( ctx.nv );
if ( ti ) { // is this also true when just "i" is given ??
CContext myctx = ctx;
myctx.reqDim = CContext::_1;
if ( ctx.want_error && sj->dy.size() )
- ret->dvd.resize( ctx.nv );
+ ret.dvd.resize( ctx.nv );
else
- ret->dvd.clear();
+ ret.dvd.clear();
for( uint ii=0; ii<ctx.nv; ++ii ){
myctx.i = ii;
uint i;
@@ -319,15 +317,15 @@ void CRef::ref_val( CTOut *ret, const CContext& ctx ) const
if( i >= sj->size() )
throw string( "i out of range" );
if ( typ == _I ) {
- ret->vd[ii] = i;
+ ret.vd[ii] = i;
} else if ( typ == _X ) {
- ret->vd[ii] = sj->x[i];
+ ret.vd[ii] = sj->x[i];
} else if ( typ == _Y ) {
- ret->vd[ii] = sj->y[i];
- if ( ret->dvd.size() )
- ret->dvd[ii] = sj->dy[i];
+ ret.vd[ii] = sj->y[i];
+ if ( ret.dvd.size() )
+ ret.dvd[ii] = sj->dy[i];
} else if ( typ == _DY ) {
- ret->vd[ii] = sj->dy[i];
+ ret.vd[ii] = sj->dy[i];
}
}
} else {
@@ -337,20 +335,20 @@ void CRef::ref_val( CTOut *ret, const CContext& ctx ) const
strg( sj->size() );
if ( typ == _I ) {
for( uint i=0; i<ctx.nv; ++i )
- ret->vd[i] = i;
+ ret.vd[i] = i;
} else if ( typ == _X ) {
- ret->vd = sj->x;
+ ret.vd = sj->x;
} else if ( typ == _Y ) {
- ret->vd = sj->y;
+ ret.vd = sj->y;
if ( ctx.want_error )
- ret->dvd = sj->dy;
+ ret.dvd = sj->dy;
} else if ( typ == _DY ) {
- ret->vd = sj->dy;
+ ret.vd = sj->dy;
}
}
}
} else if ( typ == _NI ) {
- ret->set_d( sj->size() );
+ ret.set_d( sj->size() );
} else
throw( "invalid ref(" + var_info() + ") in data file" );
return;
@@ -363,11 +361,11 @@ void CRef::ref_val( CTOut *ret, const CContext& ctx ) const
if ( typ == _FP ) {
if ( num>= fc->nPar() )
throw( "invalid p ref(" + var_info() + ") in curve file" );
- ret->set_d( cj->P[num] );
+ ret.set_d( cj->P[num] );
} else if ( typ == _FM ) {
if ( num >= CCurve::mQuality )
throw( "invalid fm ref(" + var_info() + ") in curve file" );
- ret->set_d( cj->Quality[num] );
+ ret.set_d( cj->Quality[num] );
} else
throw( "invalid ref(" + var_info() + ") in curve file" );
return;
@@ -377,16 +375,16 @@ void CRef::ref_val( CTOut *ret, const CContext& ctx ) const
//! Description of reference (for use in ??)
-void CRef::coord( CCoord *ret, uint k_in ) const
+void CRef::coord( CCoord& ret, uint k_in ) const
{
if ( typ == _K ) {
- *ret = CCoord("k","");
+ ret = CCoord("k","");
return;
} else if ( typ == _J ) {
- *ret = CCoord("j","");
+ ret = CCoord("j","");
return;
} else if ( typ == _I ) {
- *ret = CCoord("i","");
+ ret = CCoord("i","");
return;
}
@@ -411,30 +409,30 @@ void CRef::coord( CCoord *ret, uint k_in ) const
POlc fc = P2C( f );
if ( fd ) {
if ( typ == _X )
- *ret = fd->xco;
+ ret = fd->xco;
else if ( typ == _Y )
- *ret = fd->yco;
+ ret = fd->yco;
else if ( typ == _Z ) {
if (num>=fd->ZCo.size())
throw( "coord: " + strg(num) + " undefined" );
- *ret = fd->ZCo[num];
+ ret = fd->ZCo[num];
} else if ( typ == _NJ )
- *ret = CCoord("#spectra", "");
+ ret = CCoord("#spectra", "");
else if ( typ == _NI )
- *ret = CCoord("#points", "");
+ ret = CCoord("#points", "");
else
throw( "Reference " + var_info() + " not allowed in data file" );
} else if ( fc ) {
if ( typ == _FP )
- *ret = fc->PCo[num];
+ ret = fc->PCo[num];
else if ( typ == _FM )
- *ret = CCoord("fqi_"+strg(num), ""); // fit quality indicator
+ ret = CCoord("fqi_"+strg(num), ""); // fit quality indicator
else if ( typ == _Z )
- *ret = fc->ZCo[num];
+ ret = fc->ZCo[num];
else if ( typ == _FC )
- *ret = CCoord(fc->expr,"");
+ ret = CCoord(fc->expr,"");
else if ( typ == _NJ )
- *ret = CCoord("#spectra", "");
+ ret = CCoord("#spectra", "");
else
throw( "Reference " + var_info() + " not allowed in curve file" );
} else
@@ -553,7 +551,7 @@ double tree_conv_eval( double tt, void *data )
vector<double> t1(1);
t1[0] = mydata->t - tt;
mydata->cvctx.fargs = &t1;
- mydata->arg->tree_val( &(mydata->r), mydata->cvctx );
+ mydata->arg->tree_val( mydata->r, mydata->cvctx );
return mydata->r.to_d(0) * mydata->sv->intpol( tt );
}
@@ -563,7 +561,7 @@ double tree_conv_eval( double tt, void *data )
void CTree::tree_uival( uint *ret, const CContext& ctx ) const
{
CTOut val;
- tree_val( &val, ctx );
+ tree_val( val, ctx );
// since we have no integer data type, convert real by rounding
*ret = (uint) ( val.d + 0.5);
}
@@ -571,287 +569,317 @@ void CTree::tree_uival( uint *ret, const CContext& ctx ) const
//! Evaluate tree. Private. User calls pass through API implemented below.
-void CTree::tree_val( CTOut *ret, const CContext& ctx ) const
+void CTree::tree_val( CTOut& ret, const CContext& ctx ) const
{
- if ( typ == _OP ) {
- CTOut r[maxarg];
- // evaluate arguments
- for ( uint iarg=0; iarg<narg; ++iarg )
- arg[iarg]->tree_val( r+iarg, ctx );
- // scalar arguments ?
- bool is_scalar = true;
+ switch ( typ ) {
+ case _OP:
+ tree_op( ret, ctx );
+ break;
+ case _CONV: case _DIRAC:
+ tree_conv( ret, ctx );
+ break;
+ case _VAL:
+ ret.set_d( val );
+ break;
+ case _FARG:
+ if ( !(ctx.fargs) )
+ throw string( "BUG: tree_val with *fargs=0" );
+ ret.set_vd( *(ctx.fargs) );
+ break;
+ case _REF:
+ ref->ref_val( ret, ctx );
+ break;
+ case _FRF:
+ tree_frf( ret, ctx );
+ break;
+ default:
+ throw string( "BUG: unexpected tree type" );
+ }
+}
+
+
+//! A function of 1,2,.. arguments (including operators like +,-,.. ).
+
+void CTree::tree_op( CTOut& ret, const CContext& ctx ) const
+{
+ CTOut r[maxarg];
+ // evaluate arguments
+ for ( uint iarg=0; iarg<narg; ++iarg )
+ arg[iarg]->tree_val( r[iarg], ctx );
+ // scalar arguments ?
+ bool is_scalar = true;
+ for ( uint iarg=0; iarg<narg; ++iarg ) {
+ if ( r[iarg].is_scalar() )
+ continue;
+ if ( !r[iarg].is_vector() )
+ throw string( "invalid data type in func arg" );
+ is_scalar = false;
+ }
+ // size of vectorial arguments
+ int n=0;
+ if( !is_scalar ) {
for ( uint iarg=0; iarg<narg; ++iarg ) {
if ( r[iarg].is_scalar() )
continue;
- if ( !r[iarg].is_vector() )
- throw string( "invalid data type in func arg" );
- is_scalar = false;
+ if ( n==0 )
+ n = r[iarg].vd.size();
+ else if ( r[iarg].vd.size() != n )
+ throw string( "vector arguments have different size" );
}
- // size of vectorial arguments
- int n=0;
- if( !is_scalar ) {
- for ( uint iarg=0; iarg<narg; ++iarg ) {
- if ( r[iarg].is_scalar() )
- continue;
- if ( n==0 )
- n = r[iarg].vd.size();
- else if ( r[iarg].vd.size() != n )
- throw string( "vector arguments have different size" );
- }
- ret->preset_vd( n );
+ ret.preset_vd( n );
+ if( ctx.want_error )
+ ret.dvd = vector<double>( n, 0. );
+ else
+ ret.dvd.clear();
+ }
+ // now evaluate the function
+ if ( narg==1 ) {
+ if ( is_scalar )
if( ctx.want_error )
- ret->dvd = vector<double>( n, 0. );
+ ret.set_d( (*(fun->f1))( r[0].d ) );
else
- ret->dvd.clear();
+ ret.set_d( (*(fun->f1))( r[0].d ) );
+ else
+ for ( int i=0; i<n; ++i )
+ ret.vd[i] = (*(fun->f1))( r[0].vd[i] );
+ } else if ( narg==2 ) {
+ if ( r[0].is_scalar() && r[1].is_scalar() )
+ ret.set_d( (*(fun->f2))( r[0].d, r[1].d ) );
+ else if ( r[0].is_scalar() && r[1].is_vector() ) {
+ for ( int i=0; i<n; ++i )
+ ret.vd[i] = (*(fun->f2))( r[0].d, r[1].vd[i] );
+ } else if ( r[0].is_vector() && r[1].is_scalar() ) {
+ for ( int i=0; i<n; ++i )
+ ret.vd[i] = (*(fun->f2))( r[0].vd[i], r[1].d );
+ } else if ( r[0].is_vector() && r[1].is_vector() ) {
+ for ( int i=0; i<n; ++i )
+ ret.vd[i] = (*(fun->f2))( r[0].vd[i], r[1].vd[i] );
}
- // now evaluate the function
- if ( narg==1 ) {
- if ( is_scalar )
- if( ctx.want_error )
- ret->set_d( (*(fun->f1))( r[0].d ) );
- else
- ret->set_d( (*(fun->f1))( r[0].d ) );
- else
- for ( int i=0; i<n; ++i )
- ret->vd[i] = (*(fun->f1))( r[0].vd[i] );
- } else if ( narg==2 ) {
- if ( r[0].is_scalar() && r[1].is_scalar() )
- ret->set_d( (*(fun->f2))( r[0].d, r[1].d ) );
- else if ( r[0].is_scalar() && r[1].is_vector() ) {
- for ( int i=0; i<n; ++i )
- ret->vd[i] = (*(fun->f2))( r[0].d, r[1].vd[i] );
- } else if ( r[0].is_vector() && r[1].is_scalar() ) {
- for ( int i=0; i<n; ++i )
- ret->vd[i] = (*(fun->f2))( r[0].vd[i], r[1].d );
- } else if ( r[0].is_vector() && r[1].is_vector() ) {
- for ( int i=0; i<n; ++i )
- ret->vd[i] = (*(fun->f2))( r[0].vd[i], r[1].vd[i] );
- }
- } else if ( narg==3 ) {
- if ( is_scalar )
- ret->set_d( (*(fun->f3))( r[0].d, r[1].d, r[2].d ) );
- else {
- for ( int i=0; i<n; ++i )
- ret->vd[i] = (*(fun->f3))(
- r[0].to_d(i), r[1].to_d(i), r[2].to_d(i) );
- }
+ } else if ( narg==3 ) {
+ if ( is_scalar )
+ ret.set_d( (*(fun->f3))( r[0].d, r[1].d, r[2].d ) );
+ else {
+ for ( int i=0; i<n; ++i )
+ ret.vd[i] = (*(fun->f3))(
+ r[0].to_d(i), r[1].to_d(i), r[2].to_d(i) );
}
- } else if ( typ == _CONV || typ == _DIRAC ) {
- // Return values are requested for n values of _FARG ( = t ):
- if ( !(ctx.fargs) )
- throw string( "BUG: treeval(conv,dirac) with *fargs=0" );
- uint n = ctx.fargs->size();
- ret->preset_vd( n );
- ret->dvd.clear();
- // Resolution is given by context 'ctx':
- uint kconv = setConvK( ctx.k );
- uint jconv = kconv == -1 ? -1 : setConvJ( ctx.k, ctx.j, kconv );
- double conv_norm; // normalization constant (integral of resolution)
- double conv_step; // step in x/t, if equidistant
- PSpec sv;
- int nv; /* declaring this as unsigned lead to wrong wrapping */
- if( kconv!=-1 ){
- // Check whether return grid is equidistant
- if( n>1 ){
- conv_step = ((*(ctx.fargs))[n-1]-(*(ctx.fargs))[0]) / (n-1);
- if( conv_step<0 )
- conv_step = 0;
- else {
- for( int i=0; i<n-1; ++i ) {
- if( fabs((*(ctx.fargs))[i+1]-(*(ctx.fargs))[i]
- -conv_step)>1e-8*conv_step ){
- conv_step = 0; // not equidistant
- break;
- }
+ }
+}
+
+
+//! Convolution or Dirac function: conv(f, shift), dirac(shift).
+
+void CTree::tree_conv( CTOut& ret, const CContext& ctx ) const
+{
+ // Return values are requested for n values of _FARG ( = t ):
+ if ( !(ctx.fargs) )
+ throw string( "BUG: treeval(conv,dirac) with *fargs=0" );
+ uint n = ctx.fargs->size();
+ ret.preset_vd( n );
+ ret.dvd.clear();
+ // Resolution is given by context 'ctx':
+ uint kconv = setConvK( ctx.k );
+ uint jconv = kconv == -1 ? -1 : setConvJ( ctx.k, ctx.j, kconv );
+ double conv_norm; // normalization constant (integral of resolution)
+ double conv_step; // step in x/t, if equidistant
+ PSpec sv;
+ int nv; /* declaring this as unsigned lead to wrong wrapping */
+ if( kconv!=-1 ){
+ // Check whether return grid is equidistant
+ if( n>1 ){
+ conv_step = ((*(ctx.fargs))[n-1]-(*(ctx.fargs))[0]) / (n-1);
+ if( conv_step<0 )
+ conv_step = 0;
+ else {
+ for( int i=0; i<n-1; ++i ) {
+ if( fabs((*(ctx.fargs))[i+1]-(*(ctx.fargs))[i]
+ -conv_step)>1e-8*conv_step ){
+ conv_step = 0; // not equidistant
+ break;
}
}
- } else
- conv_step = 0;
- // Preset fv,sv,nv.
- POld fv = P2D( NOlm::MOM[kconv] );
- if ( !fv )
- throw "Resolution " + strg(kconv) + "is not a data file";
- sv = fv->VS(jconv);
- nv = sv->size();
- if ( nv < 2 )
- throw "Resolution " + strg(kconv) + ", spectrum " +
- strg(jconv) + ", has only " + strg(nv) + " points";
- // Calculate normalization integral.
- // At this occasion, assert that x is ordered,
- conv_norm = 0;
- double dx;
- // special treatment for first and last step:
- dx = sv->x[1] - sv->x[0];
- if ( dx <= 0 )
- throw "Resolution " + strg(kconv) + ", spectrum " +
- strg(jconv) + ", not ordered in x";
- conv_norm += sv->y[0] * dx;
- dx = sv->x[nv-1] - sv->x[nv-2];
+ }
+ } else
+ conv_step = 0;
+ // Preset fv,sv,nv.
+ POld fv = P2D( NOlm::MOM[kconv] );
+ if ( !fv )
+ throw "Resolution " + strg(kconv) + "is not a data file";
+ sv = fv->VS(jconv);
+ nv = sv->size();
+ if ( nv < 2 )
+ throw "Resolution " + strg(kconv) + ", spectrum " +
+ strg(jconv) + ", has only " + strg(nv) + " points";
+ // Calculate normalization integral.
+ // At this occasion, assert that x is ordered,
+ conv_norm = 0;
+ double dx;
+ // special treatment for first and last step:
+ dx = sv->x[1] - sv->x[0];
+ if ( dx <= 0 )
+ throw "Resolution " + strg(kconv) + ", spectrum " +
+ strg(jconv) + ", not ordered in x";
+ conv_norm += sv->y[0] * dx;
+ dx = sv->x[nv-1] - sv->x[nv-2];
+ if ( dx <= 0 )
+ throw "Resolution " + strg(kconv) + ", spectrum " +
+ strg(jconv) + ", not ordered in x";
+ conv_norm += sv->y[nv-1] * dx;
+ // the other steps:
+ for( int i=1; i<nv-1; ++i ) {
+ dx = ( sv->x[i+1] - sv->x[i-1] ) / 2;
if ( dx <= 0 )
throw "Resolution " + strg(kconv) + ", spectrum " +
strg(jconv) + ", not ordered in x";
- conv_norm += sv->y[nv-1] * dx;
- // the other steps:
- for( int i=1; i<nv-1; ++i ) {
- dx = ( sv->x[i+1] - sv->x[i-1] ) / 2;
- if ( dx <= 0 )
- throw "Resolution " + strg(kconv) + ", spectrum " +
- strg(jconv) + ", not ordered in x";
- conv_norm += sv->y[i] * dx;
- }
- // check result:
- if ( conv_norm <= 0 )
- throw "Resolution " + strg(kconv) + ", spectrum " +
- strg(jconv) + ", has non-positive norm " + strg(conv_norm);
- // Check whether x is equidistant
- if( conv_step ){
- for( int i=0; i<nv-1; ++i ) {
- if( fabs(sv->x[i+1]-sv->x[i]-conv_step)>1e-8*conv_step ){
- conv_step = 0; // not equidistant
- break;
- }
+ conv_norm += sv->y[i] * dx;
+ }
+ // check result:
+ if ( conv_norm <= 0 )
+ throw "Resolution " + strg(kconv) + ", spectrum " +
+ strg(jconv) + ", has non-positive norm " + strg(conv_norm);
+ // Check whether x is equidistant
+ if( conv_step ){
+ for( int i=0; i<nv-1; ++i ) {
+ if( fabs(sv->x[i+1]-sv->x[i]-conv_step)>1e-8*conv_step ){
+ conv_step = 0; // not equidistant
+ break;
}
}
- if( !conv_step && !allow_slow_conv )
- throw string(
- "not equidistant and slow convolution not allowed" );
}
- // The actual computation is rather different for _CONV and _DIRAC:
- if ( typ == _CONV ) {
- if( narg!=2 )
- throw string( "BUG: convolution requires 2 args" );
+ if( !conv_step && !allow_slow_conv )
+ throw string(
+ "not equidistant and slow convolution not allowed" );
+ }
+ // The actual computation is rather different for _CONV and _DIRAC:
+ if ( typ == _CONV ) {
+ if( narg!=2 )
+ throw string( "BUG: convolution requires 2 args" );
+ CTOut r;
+ arg[1]->tree_val( r, ctx );
+ if ( !r.is_scalar() )
+ throw string( "2nd arg of conv must be scalar shift" );
+ double shift = r.d;
+ if ( kconv==-1 ) {
+ // Copy convolutand (= unconvoluted theory).
CTOut r;
- arg[1]->tree_val( &r, ctx );
- if ( !r.is_scalar() )
- throw string( "2nd arg of conv must be scalar shift" );
- double shift = r.d;
- if ( kconv==-1 ) {
- // Copy convolutand (= unconvoluted theory).
- CTOut r;
- // unconvoluted theory given by arg[0]:
- arg[0]->tree_val( &r, ctx );
- for ( int i=0; i<n; ++i )
- ret->vd[i] = r.to_d(i);
- } else if ( NRead::iofmac("\\cv_mode")==-1 ) {
- CTOut r;
- CContext cvctx(ctx);
- vector<double> tt;
- cvctx.fargs = &tt;
- // Convolute theory with resolution.
- if( conv_step ){
- // Accelerated version with equidistant grids
- int ntt = n + nv - 1;
- tt.resize(ntt);
- for ( int ii=0; ii<ntt; ++ii ) // ii=(nv-1)+i-iv
- tt[ii] = (*(ctx.fargs))[0] - sv->x[0] - shift +
- (ii-(nv-1))*conv_step;
- arg[0]->tree_val( &r, cvctx );
- for ( int i=0; i<n; ++i ){
- ret->vd[i] = 0;
- for ( int iv=0; iv<nv; ++iv )
- ret->vd[i] += r.to_d(nv-1+i-iv) * sv->y[iv];
- ret->vd[i] *= conv_step / conv_norm;
- }
- } else {
- // Simplest version: non-equidistant x_out and x_res.
- for ( int i=0; i<n; ++i )
- ret->vd[i] = 0;
- tt.resize(n);
- for ( int iv=0; iv<nv; ++iv ) {
- double dx;
- if ( iv==0 )
- dx = sv->x[1] - sv->x[0];
- else if( iv==nv-1 )
- dx = sv->x[nv-1] - sv->x[nv-2];
- else
- dx = ( sv->x[iv+1] - sv->x[iv-1] ) / 2;
- for ( int i=0; i<n; ++i )
- tt[i] = (*(ctx.fargs))[i] - sv->x[iv] - shift;
- arg[0]->tree_val( &r, cvctx );
- for ( int i=0; i<n; ++i )
- ret->vd[i] += r.to_d(i) * sv->y[iv] *
- dx / conv_norm;
- }
+ // unconvoluted theory given by arg[0]:
+ arg[0]->tree_val( r, ctx );
+ for ( int i=0; i<n; ++i )
+ ret.vd[i] = r.to_d(i);
+ } else if ( NRead::iofmac("\\cv_mode")==-1 ) {
+ CTOut r;
+ CContext cvctx(ctx);
+ vector<double> tt;
+ cvctx.fargs = &tt;
+ // Convolute theory with resolution.
+ if( conv_step ){
+ // Accelerated version with equidistant grids
+ int ntt = n + nv - 1;
+ tt.resize(ntt);
+ for ( int ii=0; ii<ntt; ++ii ) // ii=(nv-1)+i-iv
+ tt[ii] = (*(ctx.fargs))[0] - sv->x[0] - shift +
+ (ii-(nv-1))*conv_step;
+ arg[0]->tree_val( r, cvctx );
+ for ( int i=0; i<n; ++i ){
+ ret.vd[i] = 0;
+ for ( int iv=0; iv<nv; ++iv )
+ ret.vd[i] += r.to_d(nv-1+i-iv) * sv->y[iv];
+ ret.vd[i] *= conv_step / conv_norm;
}
} else {
- // Convolute theory with resolution.
- // Generic integration.
- ConvDatTyp conv_data;
- conv_data.sv = sv;
- conv_data.nv = nv;
- conv_data.arg = arg[0];
- conv_data.cvctx = ctx;
- for ( int i=0; i<n; ++i ) {
- conv_data.t = (*(ctx.fargs))[i] - shift;
- 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;
+ // Simplest version: non-equidistant x_out and x_res.
+ for ( int i=0; i<n; ++i )
+ ret.vd[i] = 0;
+ tt.resize(n);
+ for ( int iv=0; iv<nv; ++iv ) {
+ double dx;
+ if ( iv==0 )
+ dx = sv->x[1] - sv->x[0];
+ else if( iv==nv-1 )
+ dx = sv->x[nv-1] - sv->x[nv-2];
+ else
+ dx = ( sv->x[iv+1] - sv->x[iv-1] ) / 2;
+ for ( int i=0; i<n; ++i )
+ tt[i] = (*(ctx.fargs))[i] - sv->x[iv] - shift;
+ arg[0]->tree_val( r, cvctx );
+ for ( int i=0; i<n; ++i )
+ ret.vd[i] += r.to_d(i) * sv->y[iv] *
+ dx / conv_norm;
}
}
-
- } else if ( typ == _DIRAC ) {
- CTOut r;
- // center of delta function given by arg[0]:
- arg[0]->tree_val( &r, ctx );
- if ( !r.is_scalar() )
- throw string( "dirac function requires scalar argument" );
- if ( kconv==-1 ) {
- // Naive representation of delta function.
- for ( int i=0; i<n; ++i )
- ret->vd[i] = fabs((*(ctx.fargs))[i]-r.to_d(i))<1e-2 ?
- 1e2 : 0;
- } else {
- // Delta function copies resolution.
- vector<double> tt(n), yy(n);
- // shift in t
- for ( int i=0; i<n; ++i )
- tt[i] = (*(ctx.fargs))[i]-r.d;
- // interpolation
- sv->intpol( tt, yy );
- // normalization
- for ( int i=0; i<n; ++i )
- ret->vd[i] = yy[i] / conv_norm;
+ } else {
+ // Convolute theory with resolution.
+ // Generic integration.
+ ConvDatTyp conv_data;
+ conv_data.sv = sv;
+ conv_data.nv = nv;
+ conv_data.arg = arg[0];
+ conv_data.cvctx = ctx;
+ for ( int i=0; i<n; ++i ) {
+ conv_data.t = (*(ctx.fargs))[i] - shift;
+ 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
- throw string( "FATAL: unexpected typ" );
- } else if ( typ == _VAL ) {
- ret->set_d( val );
- } else if ( typ == _FARG ) {
- if ( !(ctx.fargs) )
- throw string( "BUG: tree_val with *fargs=0" );
- ret->set_vd( *(ctx.fargs) );
- } else if ( typ == _REF ) {
- ref->ref_val( ret, ctx );
- } else if ( typ == _FRF ) {
- // evaluate function argument:
- CTOut rarg;
- arg[0]->tree_val( &rarg, ctx );
- vector<double> vt;
- if ( rarg.is_scalar() ) {
- vt.resize(1);
- vt[0] = rarg.d;
+ }
+
+ } else if ( typ == _DIRAC ) {
+ CTOut r;
+ // center of delta function given by arg[0]:
+ arg[0]->tree_val( r, ctx );
+ if ( !r.is_scalar() )
+ throw string( "dirac function requires scalar argument" );
+ if ( kconv==-1 ) {
+ // Naive representation of delta function.
+ for ( int i=0; i<n; ++i )
+ ret.vd[i] = fabs((*(ctx.fargs))[i]-r.to_d(i))<1e-2 ?
+ 1e2 : 0;
} else {
- vt.resize(rarg.vd.size());
- vt = rarg.vd;
+ // Delta function copies resolution.
+ vector<double> tt(n), yy(n);
+ // shift in t
+ for ( int i=0; i<n; ++i )
+ tt[i] = (*(ctx.fargs))[i]-r.d;
+ // interpolation
+ sv->intpol( tt, yy );
+ // normalization
+ for ( int i=0; i<n; ++i )
+ ret.vd[i] = yy[i] / conv_norm;
}
- // evaluate function:
- CContext cu_ctx;
- uint k = ref->get_k( ctx );
- uint j = ref->get_j( ctx );
- cu_ctx.request_vector_of_farg( k, j, &vt );
- POlc fc = P2C( NOlm::MOM[k] );
- fc->T->tree_val( ret, cu_ctx );
- if ( rarg.is_scalar() )
- ret->set_d( ret->vd[0] );
- } else if ( typ == _NOTREE )
- throw string( "BUG: tree_val(tree_typ=NOTREE)" );
- else
- throw string( "BUG: tree_val(tree_typ=UNDEF)" );
+ }
+}
+
+
+//! A curve evaluation: fc[k,j](arg).
+
+void CTree::tree_frf( CTOut& ret, const CContext& ctx ) const
+{
+ // evaluate and pack function argument:
+ CTOut rarg;
+ arg[0]->tree_val( rarg, ctx );
+ vector<double> vt;
+ if ( rarg.is_scalar() ) {
+ vt.resize(1);
+ vt[0] = rarg.d;
+ } else {
+ vt.resize(rarg.vd.size());
+ vt = rarg.vd;
+ }
+ // evaluate function:
+ CContext cu_ctx;
+ uint k = ref->get_k( ctx );
+ uint j = ref->get_j( ctx );
+ cu_ctx.request_vector_of_farg( k, j, &vt );
+ POlc fc = P2C( NOlm::MOM[k] );
+ fc->T->tree_val( ret, cu_ctx );
+ // transcribe result:
+ if ( rarg.is_scalar() )
+ ret.set_d( ret.vd[0] );
}
@@ -864,7 +892,7 @@ void CTree::tree_point_val(
ctx.request_scalar( k, j, i );
ctx.want_error = (bool)dret;
CTOut val;
- tree_val( &val, ctx );
+ tree_val( val, ctx );
if ( !val.is_scalar() )
throw "tree_point_val got " + val.info();
*ret = val.d;
@@ -880,7 +908,7 @@ void CTree::tree_vec_val(
CContext ctx;
ctx.request_vector_of_indx( k, j, ret->size() );
ctx.want_error = (bool)dret;
- tree_val( &val, ctx );
+ tree_val( val, ctx );
if( dret && !ctx.want_error )
dret->clear();
if ( val.is_vector() ){
@@ -900,7 +928,7 @@ void CTree::tree_vec_val(
//! Evaluate a curve, vectorial version. Used in fitting and in plotting.
void CTree::tree_curve_val_vec(
- vector<double> *ret, // output: the values y(t)
+ vector<double>* ret, // output: the values y(t)
const vector<double>& fargs, // input: the function arguments i
uint k, uint j // input: used to find the parameters
) const
@@ -911,13 +939,13 @@ void CTree::tree_curve_val_vec(
CContext ctx;
if ( !has_farg ) { // curve does not depend on t
ctx.request_scalar( k, j );
- tree_val( &val, ctx );
+ tree_val( val, ctx );
if ( !val.is_scalar() )
throw "tree_curve_val (sca) got " + val.info();
(*ret)[0] = val.d;
} else {
ctx.request_vector_of_farg( k, j, &fargs );
- tree_val( &val, ctx );
+ tree_val( val, ctx );
if ( !val.is_vector() )
throw "tree_curve_val (vec) got " + val.info();
*ret = val.vd;
@@ -982,30 +1010,30 @@ uint CTree::setConvJ( uint k, uint j, uint kconv )
//! Compose coordinate names and units according to the tree.
-void CTree::coord( CCoord *ret, uint k ) const
+void CTree::coord( CCoord& ret, uint k ) const
{
CCoord r[maxarg];
if ( typ == _OP ) {
for ( uint iarg=0; iarg<narg; ++iarg )
- arg[iarg]->coord( &(r[iarg]), k );
+ arg[iarg]->coord( r[iarg], k );
if ( narg==1 )
- *ret = fun->coord( r+0);
+ ret = fun->coord( r+0);
else if ( narg==2 )
- *ret = fun->coord( r+0, r+1 );
+ ret = fun->coord( r+0, r+1 );
else if ( narg==3 )
- *ret = fun->coord( r+0, r+1, r+2 );
+ ret = fun->coord( r+0, r+1, r+2 );
}
else if ( typ == _CONV )
arg[0]->coord( ret, k );
else if ( typ == _DIRAC )
- *ret = CCoord("dirac", "");
+ ret = CCoord("dirac", "");
else if ( typ == _VAL )
- *ret = CCoord(strg(val), "");
+ ret = CCoord(strg(val), "");
else if ( typ == _REF )
ref->coord(ret,k);
else if ( typ == _FRF ) {
// ref->coord(ret,k);
- *ret = CCoord( "f", "" ); // primitiver gehts nicht
+ ret = CCoord( "f", "" ); // primitiver gehts nicht
} else
throw string( "BUG: tree::coord(tree_typ=UNDEF)" );
}
diff --git a/pub/src/expr.h b/pub/src/expr.h
index 2d941ceb..c092b808 100644
--- a/pub/src/expr.h
+++ b/pub/src/expr.h
@@ -107,8 +107,8 @@ class CRef: public CVariable {
uint get_k( const CContext& ctx ) const;
uint get_j( const CContext& ctx ) const;
- void ref_val( CTOut* ret, const CContext& ctx ) const;
- void coord( CCoord *ret, uint k ) const;
+ void ref_val( CTOut& ret, const CContext& ctx ) const;
+ void coord( CCoord& ret, uint k ) const;
};
@@ -153,7 +153,10 @@ class CTree {
uint npar() const;
- void tree_val( CTOut *ret, const CContext& ctx ) const;
+ void tree_val ( CTOut& ret, const CContext& ctx ) const;
+ void tree_op ( CTOut& ret, const CContext& ctx ) const;
+ void tree_conv ( CTOut& ret, const CContext& ctx ) const;
+ void tree_frf ( CTOut& ret, const CContext& ctx ) const;
void tree_uival( uint *ret, const CContext& ctx ) const;
void tree_point_val(
double *ret, double *dret,
@@ -161,13 +164,15 @@ class CTree {
void tree_vec_val(
vector<double> *ret, vector<double> *dret,
uint k=(uint)-1, uint j=(uint)-1 ) const;
- void tree_curve_val_vec( vector<double> *ret, const vector<double>& farg,
- uint k, uint j ) const;
+ void tree_curve_val_vec(
+ vector<double> *ret, const vector<double>& farg,
+ uint k, uint j ) const;
double tree_curve_val_sca( double farg, uint k, uint j ) const;
- void coord( CCoord *ret, uint k=(uint)-1 ) const;
+ void coord( CCoord& ret, uint k=(uint)-1 ) const;
string tree_info() const;
private:
void npar_exec( uint *np ) const;
};
+
diff --git a/pub/src/opr.cpp b/pub/src/opr.cpp
index 614b4067..5d706d2d 100644
--- a/pub/src/opr.cpp
+++ b/pub/src/opr.cpp
@@ -174,7 +174,7 @@ void NOperate::Pointwise( string llabel )
throw string( "no pointwise operation on curve" );
CCoord co;
- T->coord( &co, k );
+ T->coord( co, k );
if ( lref.typ==CVariable::_X ) {
fout->xco = co;
} else if ( lref.typ==CVariable::_Y ) {
@@ -264,7 +264,7 @@ void NOperate::Integral(void)
fout->ZCo.pop_back();
fout->xco = fin->ZCo.back();
- T->coord( &(fout->yco), k );
+ T->coord( fout->yco, k );
PSpec sout( new CSpec );
sout->z = fin->V[0]->z;
--
GitLab