From c083c59c5bd4dbf022b0adf0652e2688461634d0 Mon Sep 17 00:00:00 2001
From: JWu <j.wuttke@fz-juelich.de>
Date: Fri, 26 Oct 2012 15:56:23 +0200
Subject: [PATCH] expr.cpp splitted -> node.cpp
---
pub/src/Makefile.am | 1 +
pub/src/Makefile.in | 13 +-
pub/src/expr.cpp | 745 +------------------------------------------
pub/src/node.cpp | 758 ++++++++++++++++++++++++++++++++++++++++++++
pub/src/node.h | 4 +-
5 files changed, 770 insertions(+), 751 deletions(-)
create mode 100644 pub/src/node.cpp
diff --git a/pub/src/Makefile.am b/pub/src/Makefile.am
index 463b6d39..08c340a7 100644
--- a/pub/src/Makefile.am
+++ b/pub/src/Makefile.am
@@ -68,6 +68,7 @@ mem.cpp \
mem.h \
mystd.cpp \
mystd.h \
+node.cpp \
node.h \
olf.cpp \
olf.h \
diff --git a/pub/src/Makefile.in b/pub/src/Makefile.in
index 8171f321..28d40e73 100644
--- a/pub/src/Makefile.in
+++ b/pub/src/Makefile.in
@@ -68,11 +68,12 @@ am_frida_OBJECTS = axis.$(OBJEXT) calc.$(OBJEXT) commands.$(OBJEXT) \
edif.$(OBJEXT) expr.$(OBJEXT) file_in.$(OBJEXT) \
file_ops.$(OBJEXT) file_out.$(OBJEXT) frida2.$(OBJEXT) \
func.$(OBJEXT) integrate.$(OBJEXT) list.$(OBJEXT) \
- manip.$(OBJEXT) mem.$(OBJEXT) mystd.$(OBJEXT) olf.$(OBJEXT) \
- opr.$(OBJEXT) plot.$(OBJEXT) reg.$(OBJEXT) rng.$(OBJEXT) \
- rssm.$(OBJEXT) special1.$(OBJEXT) special2.$(OBJEXT) \
- strg.$(OBJEXT) xax_lex.$(OBJEXT) xax_yacc.$(OBJEXT) \
- var.$(OBJEXT) yaml_out.$(OBJEXT) zentry.$(OBJEXT)
+ manip.$(OBJEXT) mem.$(OBJEXT) mystd.$(OBJEXT) node.$(OBJEXT) \
+ olf.$(OBJEXT) opr.$(OBJEXT) plot.$(OBJEXT) reg.$(OBJEXT) \
+ rng.$(OBJEXT) rssm.$(OBJEXT) special1.$(OBJEXT) \
+ special2.$(OBJEXT) strg.$(OBJEXT) xax_lex.$(OBJEXT) \
+ xax_yacc.$(OBJEXT) var.$(OBJEXT) yaml_out.$(OBJEXT) \
+ zentry.$(OBJEXT)
frida_OBJECTS = $(am_frida_OBJECTS)
frida_LDADD = $(LDADD)
DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir)
@@ -268,6 +269,7 @@ mem.cpp \
mem.h \
mystd.cpp \
mystd.h \
+node.cpp \
node.h \
olf.cpp \
olf.h \
@@ -402,6 +404,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/manip.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/mem.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/mystd.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/node.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/olf.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/opr.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/plot.Po@am__quote@
diff --git a/pub/src/expr.cpp b/pub/src/expr.cpp
index 41f23cfd..7fee6eac 100644
--- a/pub/src/expr.cpp
+++ b/pub/src/expr.cpp
@@ -1,12 +1,10 @@
//**************************************************************************//
//* FRIDA: flexible rapid interactive data analysis *//
-//* expr.cpp: evaluate generic expressions *//
+//* expr: evaluate generic expressions *//
//* (C) Joachim Wuttke 1990-, v2(C++) 2001- *//
//* http://apps.jcns.fz-juelich.de/frida *//
//**************************************************************************//
-#include <cstdlib> // for boost/shared_ptr
-#include <cmath>
#include <iostream>
#include <vector>
@@ -16,16 +14,8 @@
#include "olf.h"
#include "mem.h"
#include "zentry.h"
-#include "reg.h"
-#include "func.h"
#include "var.h"
#include "expr.h"
-#include "node.h"
-#include "curve.h"
-
-extern bool allow_slow_conv; // set in curve.cpp
-
-#define DEBUG 0
//***************************************************************************//
//* CContext *//
@@ -325,736 +315,3 @@ double tree_conv_eval( double tt, void *data )
mydata->arg->tree_val( mydata->res, mydata->cv_ctx );
return mydata->res.to_r(0) * mydata->sv->intpol( tt );
}
-
-
-//***************************************************************************//
-//* CNodeFun: function/operator node *//
-//***************************************************************************//
-
-CNodeFun::CNodeFun( const class CFunc *_fun, PTree _a0 )
- : fun(_fun), narg(1)
-{
- if ( fun->narg!=1 )
- throw string( "BUG: CNodeFun(1) with conflicting #arg" );
- arg[0] = _a0;
-}
-CNodeFun::CNodeFun( const class CFunc *_fun, PTree _a0, PTree _a1 )
- : fun(_fun), narg(2)
-{
- if ( fun->narg!=2 )
- throw string( "BUG: CNodeFun(2) with conflicting #arg" );
- arg[0] = _a0;
- arg[1] = _a1;
-}
-CNodeFun::CNodeFun( const class CFunc *_fun, PTree _a0, PTree _a1, PTree _a2 )
- : fun(_fun), narg(3)
-{
- if ( fun->narg!=3 )
- throw string( "BUG: CNodeFun(3) with conflicting #arg" );
- arg[0] = _a0;
- arg[1] = _a1;
- arg[2] = _a2;
-}
-
-//! A function of 1,2,.. arguments (including operators like +,-,.. ).
-
-void CNodeFun::tree_val( CResult& ret, const CContext& ctx ) const
-{
- CResult a[maxarg];
- // evaluate arguments
- for ( uint iarg=0; iarg<narg; ++iarg ) {
- arg[iarg]->tree_val( a[iarg], ctx );
- }
- // scalar arguments ?
- bool is_scalar = true;
- for ( uint iarg=0; iarg<narg; ++iarg ) {
- if ( !a[iarg].vectorial )
- continue;
- is_scalar = false;
- }
- // some input has errors ?
- bool err_input = false;
- for ( uint iarg=0; iarg<narg; ++iarg ) {
- if ( a[iarg].has_err() ){
- err_input = true;
- break;
- }
- }
- // size of vectorial arguments
- int n=0;
- if( !is_scalar ) {
- for ( uint iarg=0; iarg<narg; ++iarg ) {
- if ( !a[iarg].vectorial )
- continue;
- if ( n==0 )
- n = a[iarg].v.size();
- else if ( a[iarg].v.size() != n )
- throw "vector arguments have different size";
- }
- ret.preset_v( n );
- }
- // now evaluate the function
- if ( narg==1 ) {
- if( ctx.want_error && fun->d1 && err_input ) {
- if ( is_scalar ) {
- (*(fun->d1))( ret.r, ret.dr, a[0].r, a[0].dr );
- } else {
- ret.dv.resize( n );
- for ( int i=0; i<n; ++i )
- (*(fun->d1))( ret.v[i], ret.dv[i], a[0].v[i], a[0].dv[i] );
- }
- } else {
- if ( is_scalar )
- ret.set_r( (*(fun->f1))( a[0].r ) );
- else
- for ( int i=0; i<n; ++i )
- ret.v[i] = (*(fun->f1))( a[0].v[i] );
- }
- } else if ( narg==2 ) {
- if( ctx.want_error && ( ( fun->d2 && err_input ) ||
- fun->txt=="+-" ) ) {
- if ( is_scalar ) {
- (*(fun->d2))( ret.r, ret.dr, a[0].r, a[0].dr, a[1].r, a[1].dr );
- } else {
- ret.dv.resize( n );
- for ( int i=0; i<n; ++i )
- (*(fun->d2))( ret.v[i], ret.dv[i],
- a[0].to_r(i), a[0].to_dr(i),
- a[1].to_r(i), a[1].to_dr(i) );
- }
- } else {
- if ( !a[0].vectorial && !a[1].vectorial )
- ret.set_r( (*(fun->f2))( a[0].r, a[1].r ) );
- else if ( !a[0].vectorial && a[1].vectorial ) {
- for ( int i=0; i<n; ++i )
- ret.v[i] = (*(fun->f2))( a[0].r, a[1].v[i] );
- } else if ( a[0].vectorial && !a[1].vectorial ) {
- for ( int i=0; i<n; ++i )
- ret.v[i] = (*(fun->f2))( a[0].v[i], a[1].r );
- } else if ( a[0].vectorial && a[1].vectorial ) {
- for ( int i=0; i<n; ++i )
- ret.v[i] = (*(fun->f2))( a[0].v[i], a[1].v[i] );
- }
- }
- } else if ( narg==3 ) {
- if( ctx.want_error && fun->d3 && err_input ) {
- if ( is_scalar ) {
- (*(fun->d3))( ret.r, ret.dr, a[0].r, a[0].dr, a[1].r, a[1].dr,
- a[2].r, a[2].dr );
- } else {
- ret.dv.resize( n );
- for ( int i=0; i<n; ++i )
- (*(fun->d3))(
- ret.v[i], ret.dv[i],
- a[0].to_r(i), a[0].to_dr(i),
- a[1].to_r(i), a[1].to_dr(i),
- a[2].to_r(i), a[2].to_dr(i) );
- }
- } else {
- if ( is_scalar )
- ret.set_r( (*(fun->f3))( a[0].r, a[1].r, a[2].r ) );
- else {
- for ( int i=0; i<n; ++i )
- ret.v[i] = (*(fun->f3))(
- a[0].to_r(i), a[1].to_r(i), a[2].to_r(i) );
- }
- }
- }
- if( DEBUG ) {
- cout << "tree_op:\n";
- for ( uint iarg=0; iarg<narg; ++iarg )
- cout << " arg: " << a[iarg].result_info() << "\n";
- cout << " res: " << ret.result_info() << "\n";
- }
-}
-
-bool CNodeFun::has_dummy() const
-{
- for ( uint iarg=0; iarg<narg; ++iarg )
- if( arg[iarg]->has_dummy() )
- return true;
- return false;
-}
-
-void CNodeFun::kdep_exec( bool *kd ) const
-{
- for ( uint iarg=0; iarg<narg; ++iarg )
- arg[iarg]->kdep_exec( kd );
-}
-
-void CNodeFun::npar_exec( uint *np ) const
-{
- for ( uint iarg=0; iarg<narg; ++iarg )
- arg[iarg]->npar_exec( np );
-}
-
-void CNodeFun::set_coord( CCoord& ret, uint k ) const {
- CCoord r[maxarg];
- for ( uint iarg=0; iarg<narg; ++iarg )
- arg[iarg]->set_coord( r[iarg], k );
- if ( narg==1 )
- ret = fun->coord( r+0);
- else if ( narg==2 )
- ret = fun->coord( r+0, r+1 );
- else if ( narg==3 )
- ret = fun->coord( r+0, r+1, r+2 );
-}
-
-string CNodeFun::tree_info() const
-{
- string ret = fun->txt + "(";
- for ( uint iarg=0; iarg<narg; ++iarg )
- ret += arg[iarg]->tree_info() + (iarg<narg-1 ? "," : ")" );
- return ret;
-}
-
-//***************************************************************************//
-//* CNodeVal: end node containing a numeric value *//
-//***************************************************************************//
-
-CNodeVal::CNodeVal( double _val )
- : val( _val ) {}
-
-void CNodeVal::set_coord( CCoord& ret, uint k ) const {
- ret = CCoord(strg(val), ""); }
-string CNodeVal::tree_info() const { return "val[" + strg(val) + "]"; }
-
-//***************************************************************************//
-//* CNodeIva: indexed variable node *//
-//***************************************************************************//
-
-void CNodeIva::npar_exec( uint *np ) const
-{
- if ( var->typ == CVar::_CP )
- *np = max( *np, var->num+1 );
-}
-
-//! Evaluate reference in given context. Result returned in ret.
-
-void CNodeIva::tree_val( CResult& ret, const CContext& ctx ) const
-{
- if ( (var->typ == CVar::_I || var->typ == CVar::_J ||
- var->typ == CVar::_K) && ( ref->ti || ref->tj || ref->tk ) )
- throw "invalid cross reference: index cannot be indexed";
-
- // Get k from context:
- uint k = ref->get_k( ctx );
-
- // Return k-dependent reference (unless already done above):
- if ( var->typ == CVar::_K ) {
- ret.set_r( k );
- return;
- }
-
- // Proceed with file reference f[k]:
- POlo f = NOlm::MOM[k];
-
- if ( var->typ == CVar::_NJ ) {
- ret.set_r( f->nJ() );
- return;
- } else if ( var->typ == CVar::_R ) {
- if ( var->num>= f->RPar.size() )
- throw "invalid reference " + var->var_info();
- ret.set_r( f->RPar[var->num].val );
- return;
- }
-
- // We need j:
- uint j = ref->get_j( ctx, f->nJ() );
-
- // Return k-j-dependent reference (unless done above):
- if ( var->typ == CVar::_J ) {
- ret.set_r( ctx.j );
- return;
- } else if ( var->typ == CVar::_Z ) {
- if ( var->num>= f->nZ() )
- throw "invalid reference " + var->var_info();
- ret.set_r( f->V[j]->z[var->num] );
- return;
- }
-
- // References for data file
- POld fd = P2D( f );
- if ( fd ) {
- PSpec sj = fd->VS(j);
- if ( var->pointwise() ) {
- if ( var->typ == CVar::_DY && !sj->dy.size() )
- throw string( "dy not set" );
- if ( ctx.dim==CContext::_1 ) {
- uint i;
- if ( ref->ti ) {
- ref->ti->tree_uival( &i, ctx );
- } else {
- i = ctx.i;
- if ( i==(uint)-1 )
- throw string( "missing index i" );
- }
- if( i>= sj->size() )
- throw "i=" + strg(i) + " exceeds scan size" +
- strg(sj->size());
- if ( var->typ == CVar::_I ) {
- ret.set_r( i );
- } else if ( var->typ == CVar::_X ) {
- ret.set_r( sj->x[i] );
- } else if ( var->typ == CVar::_Y ) {
- ret.set_r( sj->y[i] );
- if ( ctx.want_error && sj->dy.size() )
- ret.dr = sj->dy[i];
- } else if ( var->typ == CVar::_DY ) {
- ret.set_r( sj->dy[i] );
- }
- } else if ( ctx.dim==CContext::_VI ) {
- if ( ref->ti ) {
- CContext myctx = ctx;
- myctx.dim = CContext::_1;
- ret.preset_v( ctx.nv );
- if ( var->typ == CVar::_Y && ctx.want_error &&
- sj->dy.size() )
- ret.dv.resize( ctx.nv );
- for( uint ii=0; ii<ctx.nv; ++ii ){
- myctx.i = ii;
- uint i;
- ref->ti->tree_uival( &i, myctx );
- if( i >= sj->size() )
- throw "i=" + strg(i) + " exceeds scan size" +
- strg(sj->size());
- if ( var->typ == CVar::_I ) {
- ret.v[ii] = i;
- } else if ( var->typ == CVar::_X ) {
- ret.v[ii] = sj->x[i];
- } else if ( var->typ == CVar::_Y ) {
- ret.v[ii] = sj->y[i];
- if ( ret.dv.size() )
- ret.dv[ii] = sj->dy[i];
- } else if ( var->typ == CVar::_DY ) {
- ret.v[ii] = sj->dy[i];
- }
- }
- } else {
- if ( ctx.nv != sj->size() )
- throw "ref_val " + var->var_info() + ": requested vec(" +
- strg(ctx.nv) + "), found vec(" +
- strg( sj->size() );
- if ( var->typ == CVar::_I ) {
- ret.preset_v( ctx.nv );
- for( uint i=0; i<ctx.nv; ++i )
- ret.v[i] = i;
- } else if ( var->typ == CVar::_X ) {
- ret.set_v( sj->x );
- } else if ( var->typ == CVar::_Y ) {
- ret.set_v( sj->y );
- if ( ctx.want_error )
- ret.dv = sj->dy;
- } else if ( var->typ == CVar::_DY ) {
- ret.set_v( sj->dy );
- }
- }
- } else
- throw string( "BUG: invalid context::dim" );
- } else if ( var->typ == CVar::_NI ) {
- ret.set_r( sj->size() );
- } else
- throw "invalid ref(" + var->var_info() + ") in data file";
- return;
- }
-
- // References for curve file:
- POlc fc = P2C( f );
- if ( fc ) {
- PCurve cj = fc->VC(j);
- if ( var->typ == CVar::_CP ) {
- if ( var->num>= fc->nPar() )
- throw "invalid p ref(" + var->var_info() + ") in curve file";
- ret.set_r( cj->P[var->num] );
- } else if ( var->typ == CVar::_CQ ) {
- if ( var->num >= CCurve::mQuality )
- throw "invalid fm ref(" + var->var_info() + ") in curve file";
- ret.set_r( cj->Quality[var->num] );
- } else
- throw "invalid ref(" + var->var_info() + ") in curve file";
- return;
- }
-}
-
-//! Description of reference (for use in ??)
-
-void CNodeIva::set_coord( CCoord& ret, uint k_in ) const
-{
- if ( var->typ == CVar::_K ) {
- ret = CCoord("k","");
- return;
- } else if ( var->typ == CVar::_J ) {
- ret = CCoord("j","");
- return;
- } else if ( var->typ == CVar::_I ) {
- ret = CCoord("i","");
- return;
- }
-
- uint k;
- if (ref->tk) {
- // get k from reference :
- CContext ctx( k_in );
- ref->tk->tree_uival(&k, ctx);
- if (k>=NOlm::MOM.size())
- throw "set_coord: file index " + strg(k) + "invalid, maximum is " +
- strg(NOlm::MOM.size()-1);
- } else {
- k = k_in;
- }
-
- POlo f = NOlm::MOM[k];
- if (!f)
- throw "set_coord: ref where not refering to a file";
-
- if ( var->typ == CVar::_Z ) {
- if (var->num>=f->ZCo.size())
- throw "set_coord: z" + strg(var->num) + " undefined";
- ret = f->ZCo[var->num];
- } else if ( var->typ == CVar::_R ) {
- if (var->num>=f->RPar.size())
- throw "set_coord: r" + strg(var->num) + " undefined";
- ret = f->RPar[var->num].Co;
- } else if ( var->typ == CVar::_NJ ) {
- ret = CCoord("#spectra", "");
- } else {
- POld fd = P2D( f );
- POlc fc = P2C( f );
- if ( fd ) {
- if ( var->typ == CVar::_X )
- ret = fd->xco;
- else if ( var->typ == CVar::_Y )
- ret = fd->yco;
- else if ( var->typ == CVar::_DY )
- ret = CCoord( "d"+fd->yco.name, fd->yco.unit );
- else if ( var->typ == CVar::_NI )
- ret = CCoord("#points", "");
- else
- throw "reference " + var->var_info() +
- " not allowed in data file";
- } else if ( fc ) {
- if ( var->typ == CVar::_CP ) {
- if( var->num>=fc->nPar() )
- throw "invalid reference p" + strg(var->num);
- ret = fc->PCo[var->num];
- } else if ( var->typ == CVar::_CQ )
- ret = CCoord("cq"+strg(var->num), ""); // fit quality indicator
- else if ( var->typ == CVar::_C )
- ret = CCoord(fc->expr,"");
- else
- throw "reference " + var->var_info() +
- " not allowed in curve file";
- } else
- throw "BUG: ref->set_coord unexpected else";
- }
-}
-
-//***************************************************************************//
-//* CNodeRgr: register reference node *//
-//***************************************************************************//
-
-CNodeRgr::CNodeRgr( PRgr& _rgr )
- : rgr( _rgr ) {}
-
-void CNodeRgr::tree_val( CResult& ret, const CContext& ctx ) const {
- ret.set_r( NReg::lastresult ); }
-void CNodeRgr::set_coord( CCoord& ret, uint k ) const {
- ret = CCoord( "$", "" ); } // primitiver gehts nicht
-
-//***************************************************************************//
-//* CNodeDummy: function argument (t) node *//
-//***************************************************************************//
-
-//***************************************************************************//
-//* CNodeCev: curve evaluation node *//
-//***************************************************************************//
-
-CNodeCev::CNodeCev( PRef& _ref, PTree _arg )
- : ref( _ref ), arg( _arg )
-{
- if( ref->ti )
- throw "curve reference must not depend on index i";
-}
-
-//! A curve evaluation: fc[k,j](arg).
-
-void CNodeCev::tree_val( CResult& ret, const CContext& ctx ) const
-{
- // get k from context:
- uint k = ref->get_k( ctx );
-
- // k points to a curve:
- POlc fc = P2C( NOlm::MOM[k] );
- if( !fc )
- throw "file " + strg(k) + " is not a curve";
-
- // get j from context:
- uint j = ref->get_j( ctx, fc->nJ() );
-
- // evaluate and pack function argument:
- CResult farg;
- arg->tree_val( farg, ctx );
-
- // evaluate the curve:
- if ( farg.vectorial ) {
- ret.preset_v( farg.v.size() );
- fc->curve_val_vec( &(ret.v), farg.v, k, j );
- } else
- ret.set_r( fc->curve_val_sca( farg.r, k, j ) );
-}
-
-void CNodeCev::set_coord( CCoord& ret, uint k ) const {
- ret = CCoord( "c", "" ); } // primitiver gehts nicht
-
-
-//***************************************************************************//
-//* CNodeConv1: convolution or dirac, with argument: shift *//
-//***************************************************************************//
-
-void CNodeConv1::tree_val( CResult& ret, const CContext& ctx ) const
-{
- // Shift:
- CResult res_shift;
- shift->tree_val( res_shift, ctx );
- if ( res_shift.vectorial )
- throw "shift is not a scalar";
- double theshift = res_shift.r;
-
- // Return values are requested for n values of _DUMMY ( = t ):
- if ( !(ctx.vt) )
- throw "BUG: convolution tree_val with *vt=0";
- uint n = ctx.vt->size();
- ret.preset_v( n );
- ret.dv.clear();
-
- // Resolution is given by context 'ctx':
- uint kconv, jconv;
- PSpec sv;
- NCurveFile::setConv( &sv, &kconv, &jconv, ctx.k, ctx.j );
-
- if ( kconv==-1 ) {
- copy_theory( ret, ctx, theshift );
- } else {
- double conv_norm; // normalization constant (integral of resolution)
- double conv_step; // step in x/t, if equidistant
- string info = "resolution " + strg(kconv) + ", spectrum " + strg(jconv);
- conv_norm = sv->norm( info );
- conv_step = sv->has_step( ctx.vt_step ) ? ctx.vt_step : 0;
- if( !conv_step && !allow_slow_conv )
- throw "grid is not equidistant; slow convolution is not allowed";
- convolve( ret, ctx, theshift, sv, conv_norm, conv_step );
- }
-}
-
-//***************************************************************************//
-//* CNodeConv2: convolution with two arguments: theory, shift *//
-//***************************************************************************//
-
-void CNodeConv2::kdep_exec( bool *kd ) const
-{
- theory->kdep_exec( kd );
- shift->kdep_exec( kd );
-}
-
-void CNodeConv2::npar_exec( uint *np ) const
-{
- theory->npar_exec( np );
- shift->npar_exec( np );
-}
-
-
-//***************************************************************************//
-//* CNodeConv: convoluted function node *//
-//***************************************************************************//
-
-void CNodeConv::copy_theory( CResult& ret, const CContext& ctx,
- double theshift ) const
-{
- CResult res_theory;
- theory->tree_val( res_theory, ctx );
- for ( int i=0; i<ctx.vt->size(); ++i )
- ret.v[i] = res_theory.to_r(i);
-}
-
-void CNodeConv::convolve( CResult& ret, const CContext& ctx,
- double theshift, const PSpec& sv,
- double conv_norm, double conv_step ) const
-{
- uint n = ctx.vt->size();
- uint nv = sv->size();
- CResult res_theory;
- CContext cv_ctx(ctx);
- vector<double> tt; // theory grid
- cv_ctx.vt = &tt;
- // Convolute theory with resolution.
- if( conv_step ){ // accelerated version with equidistant grids
- // precompute theory on unified grid
- int ntt = n + nv - 1;
- tt.resize(ntt);
- double tt0 = (*(ctx.vt))[0] - sv->x[0] - theshift - (nv-1)*conv_step;
- for ( int ii=0; ii<ntt; ++ii ) // ii=(nv-1)+i-iv
- tt[ii] = tt0 + ii*conv_step;
- theory->tree_val( res_theory, cv_ctx );
- // now convolve
- for ( int i=0; i<n; ++i ){
- ret.v[i] = 0;
- for ( int iv=0; iv<nv; ++iv )
- ret.v[i] += res_theory.to_r(nv-1+i-iv) * sv->y[iv];
- ret.v[i] *= conv_step / conv_norm;
- }
- } else { // simplest version: non-equidistant x_out and x_res.
- for ( int i=0; i<n; ++i )
- ret.v[i] = 0;
- tt.resize(n);
- for ( int iv=0; iv<nv; ++iv ) {
- double dx; // variable step
- 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.vt))[i] - sv->x[iv] - theshift;
- theory->tree_val( res_theory, cv_ctx );
- for ( int i=0; i<n; ++i )
- ret.v[i] += res_theory.to_r(i) * sv->y[iv] *
- dx / conv_norm;
- }
- }
-}
-
-
-//***************************************************************************//
-//* CNodePConv: convolution with primitive *//
-//***************************************************************************//
-
-void CNodePConv::copy_theory( CResult& ret, const CContext& ctx,
- double theshift ) const
-{
- throw "non-convolution not implemented for pconv";
- // to implement this efficiently, I would like to have CResult::operator+
- // for which we should go to shared pointers ...
-}
-
-void CNodePConv::convolve( CResult& ret, const CContext& ctx,
- double theshift, const PSpec& sv,
- double conv_norm, double conv_step ) const
-{
- uint n = ctx.vt->size();
- uint nv = sv->size();
- CContext cv_ctx(ctx);
- vector<double> tt;
- cv_ctx.vt = &tt;
- // Convolute theory with resolution.
- if( conv_step ){ // accelerated version with equidistant grids
- // precompute primitive on unified grid
- int ntt = n + nv;
- tt.resize(ntt);
- double tt0 = (*(ctx.vt))[0] - sv->x[0] - theshift - (nv-0.5)*conv_step;
- for ( int ii=0; ii<ntt; ++ii ) // ii=(nv-1)+i-iv
- tt[ii] = tt0 + ii*conv_step;
- CResult res_theory;
- theory->tree_val( res_theory, cv_ctx );
- // now convolve resolution with differential of primitive
- for ( int i=0; i<n; ++i ){
- ret.v[i] = 0;
- for ( int iv=0; iv<nv; ++iv )
- ret.v[i] += (res_theory.to_r(nv +i-iv)-
- res_theory.to_r(nv-1+i-iv)) * sv->y[iv];
- ret.v[i] /= conv_norm;
- }
- } else { // simplest version: non-equidistant x_out and x_res.
-/* DOES NOT WORK
- double tt0;
- CResult res_theory;
- tt.resize(nv);
- if( nv<2 )
- throw "non-equidistant convolution requires nv>=2";
- for ( int i=0; i<n; ++i ){
- tt0 = (*(ctx.vt))[i] - theshift;
- for ( int iv=0; iv<nv; ++iv )
- tt[iv] = tt0 - sv->x[iv];
- theory->tree_val( res_theory, cv_ctx );
- ret.v[i] = 0;
- for ( int iv=0; iv<nv-1; ++iv ) {
- ret.v[i] += (sv->y[iv+1]+sv->y[iv])*
- (sv->x[iv+1]-sv->x[iv])*
- (res_theory.to_r(iv+1) - res_theory.to_r(iv) );
- }
- ret.v[i] /= 2*conv_norm;
- }
-*/
- for ( int i=0; i<n; ++i )
- ret.v[i] = 0;
- tt.resize(n);
- CResult res_theory_low, res_theory_hig;
- double dx;
- for ( int iv=0; iv<nv; ++iv ) {
- 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.vt))[i] - sv->x[iv] - theshift - dx/2;
- theory->tree_val( res_theory_low, cv_ctx );
- }
- for ( int i=0; i<n; ++i ) {
- tt[i] = (*(ctx.vt))[i] - sv->x[iv] - theshift + dx/2;
- theory->tree_val( res_theory_hig, cv_ctx );
- }
- for ( int i=0; i<n; ++i )
- ret.v[i] += (res_theory_hig.to_r(i)-res_theory_low.to_r(i)) *
- sv->y[iv] * dx;
- }
- for ( int i=0; i<n; ++i )
- ret.v[i] /= conv_norm;
- }
-}
-
-
-//***************************************************************************//
-//* CNodeDirac: node that copies the resolution function *//
-//***************************************************************************//
-
-void CNodeDirac::kdep_exec( bool *kd ) const
-{
- shift->kdep_exec( kd );
-}
-
-void CNodeDirac::npar_exec( uint *np ) const
-{
- shift->npar_exec( np );
-}
-
-//! Convolution or Dirac function: conv(f, shift), dirac(shift).
-
-void CNodeDirac::convolve( CResult& ret, const CContext& ctx,
- double theshift, const PSpec& sv,
- double conv_norm, double conv_step ) const
-{
- uint n = ctx.vt->size();
- // Delta function copies resolution.
- vector<double> tt(n), yy(n);
- // shift in t
- for ( int i=0; i<n; ++i )
- tt[i] = (*(ctx.vt))[i]-theshift;
- // interpolation
- sv->intpol( tt, yy );
- // normalization
- for ( int i=0; i<n; ++i )
- ret.v[i] = yy[i] / conv_norm;
-}
-
-void CNodeDirac::copy_theory( CResult& ret, const CContext& ctx,
- double theshift ) const
-{
- // Naive representation of delta function.
- for ( int i=0; i<ctx.vt->size(); ++i )
- ret.v[i] = fabs((*(ctx.vt))[i]-theshift)<1e-2 ? 1e2 : 0;
-}
-
-void CNodeDirac::set_coord( CCoord& ret, uint k ) const {
- ret = CCoord("resol", ""); }
diff --git a/pub/src/node.cpp b/pub/src/node.cpp
new file mode 100644
index 00000000..015b6a74
--- /dev/null
+++ b/pub/src/node.cpp
@@ -0,0 +1,758 @@
+//**************************************************************************//
+//* FRIDA: flexible rapid interactive data analysis *//
+//* node: nodes within expression trees *//
+//* (C) Joachim Wuttke 1990-, v2(C++) 2001- *//
+//* http://apps.jcns.fz-juelich.de/frida *//
+//**************************************************************************//
+
+#include <iostream>
+#include <vector>
+
+#include "defs.h"
+#include "strg.h"
+#include "olf.h"
+#include "mem.h"
+#include "zentry.h"
+#include "reg.h"
+#include "func.h"
+#include "var.h"
+#include "expr.h"
+#include "node.h"
+#include "curve.h"
+
+extern bool allow_slow_conv; // set in curve.cpp
+
+#define DEBUG 0
+
+
+//***************************************************************************//
+//* CNodeFun: function/operator node *//
+//***************************************************************************//
+
+CNodeFun::CNodeFun( const class CFunc *_fun, PTree _a0 )
+ : fun(_fun), narg(1)
+{
+ if ( fun->narg!=1 )
+ throw string( "BUG: CNodeFun(1) with conflicting #arg" );
+ arg[0] = _a0;
+}
+CNodeFun::CNodeFun( const class CFunc *_fun, PTree _a0, PTree _a1 )
+ : fun(_fun), narg(2)
+{
+ if ( fun->narg!=2 )
+ throw string( "BUG: CNodeFun(2) with conflicting #arg" );
+ arg[0] = _a0;
+ arg[1] = _a1;
+}
+CNodeFun::CNodeFun( const class CFunc *_fun, PTree _a0, PTree _a1, PTree _a2 )
+ : fun(_fun), narg(3)
+{
+ if ( fun->narg!=3 )
+ throw string( "BUG: CNodeFun(3) with conflicting #arg" );
+ arg[0] = _a0;
+ arg[1] = _a1;
+ arg[2] = _a2;
+}
+
+//! A function of 1,2,.. arguments (including operators like +,-,.. ).
+
+void CNodeFun::tree_val( CResult& ret, const CContext& ctx ) const
+{
+ CResult a[maxarg];
+ // evaluate arguments
+ for ( uint iarg=0; iarg<narg; ++iarg ) {
+ arg[iarg]->tree_val( a[iarg], ctx );
+ }
+ // scalar arguments ?
+ bool is_scalar = true;
+ for ( uint iarg=0; iarg<narg; ++iarg ) {
+ if ( !a[iarg].vectorial )
+ continue;
+ is_scalar = false;
+ }
+ // some input has errors ?
+ bool err_input = false;
+ for ( uint iarg=0; iarg<narg; ++iarg ) {
+ if ( a[iarg].has_err() ){
+ err_input = true;
+ break;
+ }
+ }
+ // size of vectorial arguments
+ int n=0;
+ if( !is_scalar ) {
+ for ( uint iarg=0; iarg<narg; ++iarg ) {
+ if ( !a[iarg].vectorial )
+ continue;
+ if ( n==0 )
+ n = a[iarg].v.size();
+ else if ( a[iarg].v.size() != n )
+ throw "vector arguments have different size";
+ }
+ ret.preset_v( n );
+ }
+ // now evaluate the function
+ if ( narg==1 ) {
+ if( ctx.want_error && fun->d1 && err_input ) {
+ if ( is_scalar ) {
+ (*(fun->d1))( ret.r, ret.dr, a[0].r, a[0].dr );
+ } else {
+ ret.dv.resize( n );
+ for ( int i=0; i<n; ++i )
+ (*(fun->d1))( ret.v[i], ret.dv[i], a[0].v[i], a[0].dv[i] );
+ }
+ } else {
+ if ( is_scalar )
+ ret.set_r( (*(fun->f1))( a[0].r ) );
+ else
+ for ( int i=0; i<n; ++i )
+ ret.v[i] = (*(fun->f1))( a[0].v[i] );
+ }
+ } else if ( narg==2 ) {
+ if( ctx.want_error && ( ( fun->d2 && err_input ) ||
+ fun->txt=="+-" ) ) {
+ if ( is_scalar ) {
+ (*(fun->d2))( ret.r, ret.dr, a[0].r, a[0].dr, a[1].r, a[1].dr );
+ } else {
+ ret.dv.resize( n );
+ for ( int i=0; i<n; ++i )
+ (*(fun->d2))( ret.v[i], ret.dv[i],
+ a[0].to_r(i), a[0].to_dr(i),
+ a[1].to_r(i), a[1].to_dr(i) );
+ }
+ } else {
+ if ( !a[0].vectorial && !a[1].vectorial )
+ ret.set_r( (*(fun->f2))( a[0].r, a[1].r ) );
+ else if ( !a[0].vectorial && a[1].vectorial ) {
+ for ( int i=0; i<n; ++i )
+ ret.v[i] = (*(fun->f2))( a[0].r, a[1].v[i] );
+ } else if ( a[0].vectorial && !a[1].vectorial ) {
+ for ( int i=0; i<n; ++i )
+ ret.v[i] = (*(fun->f2))( a[0].v[i], a[1].r );
+ } else if ( a[0].vectorial && a[1].vectorial ) {
+ for ( int i=0; i<n; ++i )
+ ret.v[i] = (*(fun->f2))( a[0].v[i], a[1].v[i] );
+ }
+ }
+ } else if ( narg==3 ) {
+ if( ctx.want_error && fun->d3 && err_input ) {
+ if ( is_scalar ) {
+ (*(fun->d3))( ret.r, ret.dr, a[0].r, a[0].dr, a[1].r, a[1].dr,
+ a[2].r, a[2].dr );
+ } else {
+ ret.dv.resize( n );
+ for ( int i=0; i<n; ++i )
+ (*(fun->d3))(
+ ret.v[i], ret.dv[i],
+ a[0].to_r(i), a[0].to_dr(i),
+ a[1].to_r(i), a[1].to_dr(i),
+ a[2].to_r(i), a[2].to_dr(i) );
+ }
+ } else {
+ if ( is_scalar )
+ ret.set_r( (*(fun->f3))( a[0].r, a[1].r, a[2].r ) );
+ else {
+ for ( int i=0; i<n; ++i )
+ ret.v[i] = (*(fun->f3))(
+ a[0].to_r(i), a[1].to_r(i), a[2].to_r(i) );
+ }
+ }
+ }
+ if( DEBUG ) {
+ cout << "tree_op:\n";
+ for ( uint iarg=0; iarg<narg; ++iarg )
+ cout << " arg: " << a[iarg].result_info() << "\n";
+ cout << " res: " << ret.result_info() << "\n";
+ }
+}
+
+bool CNodeFun::has_dummy() const
+{
+ for ( uint iarg=0; iarg<narg; ++iarg )
+ if( arg[iarg]->has_dummy() )
+ return true;
+ return false;
+}
+
+void CNodeFun::kdep_exec( bool *kd ) const
+{
+ for ( uint iarg=0; iarg<narg; ++iarg )
+ arg[iarg]->kdep_exec( kd );
+}
+
+void CNodeFun::npar_exec( uint *np ) const
+{
+ for ( uint iarg=0; iarg<narg; ++iarg )
+ arg[iarg]->npar_exec( np );
+}
+
+void CNodeFun::set_coord( CCoord& ret, uint k ) const {
+ CCoord r[maxarg];
+ for ( uint iarg=0; iarg<narg; ++iarg )
+ arg[iarg]->set_coord( r[iarg], k );
+ if ( narg==1 )
+ ret = fun->coord( r+0);
+ else if ( narg==2 )
+ ret = fun->coord( r+0, r+1 );
+ else if ( narg==3 )
+ ret = fun->coord( r+0, r+1, r+2 );
+}
+
+string CNodeFun::tree_info() const
+{
+ string ret = fun->txt + "(";
+ for ( uint iarg=0; iarg<narg; ++iarg )
+ ret += arg[iarg]->tree_info() + (iarg<narg-1 ? "," : ")" );
+ return ret;
+}
+
+//***************************************************************************//
+//* CNodeVal: end node containing a numeric value *//
+//***************************************************************************//
+
+CNodeVal::CNodeVal( double _val )
+ : val( _val ) {}
+
+void CNodeVal::set_coord( CCoord& ret, uint k ) const {
+ ret = CCoord(strg(val), ""); }
+string CNodeVal::tree_info() const { return "val[" + strg(val) + "]"; }
+
+//***************************************************************************//
+//* CNodeIva: indexed variable node *//
+//***************************************************************************//
+
+void CNodeIva::npar_exec( uint *np ) const
+{
+ if ( var->typ == CVar::_CP )
+ *np = max( *np, var->num+1 );
+}
+
+//! Evaluate reference in given context. Result returned in ret.
+
+void CNodeIva::tree_val( CResult& ret, const CContext& ctx ) const
+{
+ if ( (var->typ == CVar::_I || var->typ == CVar::_J ||
+ var->typ == CVar::_K) && ( ref->ti || ref->tj || ref->tk ) )
+ throw "invalid cross reference: index cannot be indexed";
+
+ // Get k from context:
+ uint k = ref->get_k( ctx );
+
+ // Return k-dependent reference (unless already done above):
+ if ( var->typ == CVar::_K ) {
+ ret.set_r( k );
+ return;
+ }
+
+ // Proceed with file reference f[k]:
+ POlo f = NOlm::MOM[k];
+
+ if ( var->typ == CVar::_NJ ) {
+ ret.set_r( f->nJ() );
+ return;
+ } else if ( var->typ == CVar::_R ) {
+ if ( var->num>= f->RPar.size() )
+ throw "invalid reference " + var->var_info();
+ ret.set_r( f->RPar[var->num].val );
+ return;
+ }
+
+ // We need j:
+ uint j = ref->get_j( ctx, f->nJ() );
+
+ // Return k-j-dependent reference (unless done above):
+ if ( var->typ == CVar::_J ) {
+ ret.set_r( ctx.j );
+ return;
+ } else if ( var->typ == CVar::_Z ) {
+ if ( var->num>= f->nZ() )
+ throw "invalid reference " + var->var_info();
+ ret.set_r( f->V[j]->z[var->num] );
+ return;
+ }
+
+ // References for data file
+ POld fd = P2D( f );
+ if ( fd ) {
+ PSpec sj = fd->VS(j);
+ if ( var->pointwise() ) {
+ if ( var->typ == CVar::_DY && !sj->dy.size() )
+ throw string( "dy not set" );
+ if ( ctx.dim==CContext::_1 ) {
+ uint i;
+ if ( ref->ti ) {
+ ref->ti->tree_uival( &i, ctx );
+ } else {
+ i = ctx.i;
+ if ( i==(uint)-1 )
+ throw string( "missing index i" );
+ }
+ if( i>= sj->size() )
+ throw "i=" + strg(i) + " exceeds scan size" +
+ strg(sj->size());
+ if ( var->typ == CVar::_I ) {
+ ret.set_r( i );
+ } else if ( var->typ == CVar::_X ) {
+ ret.set_r( sj->x[i] );
+ } else if ( var->typ == CVar::_Y ) {
+ ret.set_r( sj->y[i] );
+ if ( ctx.want_error && sj->dy.size() )
+ ret.dr = sj->dy[i];
+ } else if ( var->typ == CVar::_DY ) {
+ ret.set_r( sj->dy[i] );
+ }
+ } else if ( ctx.dim==CContext::_VI ) {
+ if ( ref->ti ) {
+ CContext myctx = ctx;
+ myctx.dim = CContext::_1;
+ ret.preset_v( ctx.nv );
+ if ( var->typ == CVar::_Y && ctx.want_error &&
+ sj->dy.size() )
+ ret.dv.resize( ctx.nv );
+ for( uint ii=0; ii<ctx.nv; ++ii ){
+ myctx.i = ii;
+ uint i;
+ ref->ti->tree_uival( &i, myctx );
+ if( i >= sj->size() )
+ throw "i=" + strg(i) + " exceeds scan size" +
+ strg(sj->size());
+ if ( var->typ == CVar::_I ) {
+ ret.v[ii] = i;
+ } else if ( var->typ == CVar::_X ) {
+ ret.v[ii] = sj->x[i];
+ } else if ( var->typ == CVar::_Y ) {
+ ret.v[ii] = sj->y[i];
+ if ( ret.dv.size() )
+ ret.dv[ii] = sj->dy[i];
+ } else if ( var->typ == CVar::_DY ) {
+ ret.v[ii] = sj->dy[i];
+ }
+ }
+ } else {
+ if ( ctx.nv != sj->size() )
+ throw "ref_val " + var->var_info() +
+ ": requested vec(" + strg(ctx.nv) +
+ "), found vec(" + strg( sj->size() );
+ if ( var->typ == CVar::_I ) {
+ ret.preset_v( ctx.nv );
+ for( uint i=0; i<ctx.nv; ++i )
+ ret.v[i] = i;
+ } else if ( var->typ == CVar::_X ) {
+ ret.set_v( sj->x );
+ } else if ( var->typ == CVar::_Y ) {
+ ret.set_v( sj->y );
+ if ( ctx.want_error )
+ ret.dv = sj->dy;
+ } else if ( var->typ == CVar::_DY ) {
+ ret.set_v( sj->dy );
+ }
+ }
+ } else
+ throw string( "BUG: invalid context::dim" );
+ } else if ( var->typ == CVar::_NI ) {
+ ret.set_r( sj->size() );
+ } else
+ throw "invalid ref(" + var->var_info() + ") in data file";
+ return;
+ }
+
+ // References for curve file:
+ POlc fc = P2C( f );
+ if ( fc ) {
+ PCurve cj = fc->VC(j);
+ if ( var->typ == CVar::_CP ) {
+ if ( var->num>= fc->nPar() )
+ throw "invalid p ref(" + var->var_info() + ") in curve file";
+ ret.set_r( cj->P[var->num] );
+ } else if ( var->typ == CVar::_CQ ) {
+ if ( var->num >= CCurve::mQuality )
+ throw "invalid fm ref(" + var->var_info() + ") in curve file";
+ ret.set_r( cj->Quality[var->num] );
+ } else
+ throw "invalid ref(" + var->var_info() + ") in curve file";
+ return;
+ }
+}
+
+//! Description of reference (for use in ??)
+
+void CNodeIva::set_coord( CCoord& ret, uint k_in ) const
+{
+ if ( var->typ == CVar::_K ) {
+ ret = CCoord("k","");
+ return;
+ } else if ( var->typ == CVar::_J ) {
+ ret = CCoord("j","");
+ return;
+ } else if ( var->typ == CVar::_I ) {
+ ret = CCoord("i","");
+ return;
+ }
+
+ uint k;
+ if (ref->tk) {
+ // get k from reference :
+ CContext ctx( k_in );
+ ref->tk->tree_uival(&k, ctx);
+ if (k>=NOlm::MOM.size())
+ throw "set_coord: file index " + strg(k) + "invalid, maximum is " +
+ strg(NOlm::MOM.size()-1);
+ } else {
+ k = k_in;
+ }
+
+ POlo f = NOlm::MOM[k];
+ if (!f)
+ throw "set_coord: ref where not refering to a file";
+
+ if ( var->typ == CVar::_Z ) {
+ if (var->num>=f->ZCo.size())
+ throw "set_coord: z" + strg(var->num) + " undefined";
+ ret = f->ZCo[var->num];
+ } else if ( var->typ == CVar::_R ) {
+ if (var->num>=f->RPar.size())
+ throw "set_coord: r" + strg(var->num) + " undefined";
+ ret = f->RPar[var->num].Co;
+ } else if ( var->typ == CVar::_NJ ) {
+ ret = CCoord("#spectra", "");
+ } else {
+ POld fd = P2D( f );
+ POlc fc = P2C( f );
+ if ( fd ) {
+ if ( var->typ == CVar::_X )
+ ret = fd->xco;
+ else if ( var->typ == CVar::_Y )
+ ret = fd->yco;
+ else if ( var->typ == CVar::_DY )
+ ret = CCoord( "d"+fd->yco.name, fd->yco.unit );
+ else if ( var->typ == CVar::_NI )
+ ret = CCoord("#points", "");
+ else
+ throw "reference " + var->var_info() +
+ " not allowed in data file";
+ } else if ( fc ) {
+ if ( var->typ == CVar::_CP ) {
+ if( var->num>=fc->nPar() )
+ throw "invalid reference p" + strg(var->num);
+ ret = fc->PCo[var->num];
+ } else if ( var->typ == CVar::_CQ )
+ ret = CCoord("cq"+strg(var->num), ""); // fit quality indicator
+ else if ( var->typ == CVar::_C )
+ ret = CCoord(fc->expr,"");
+ else
+ throw "reference " + var->var_info() +
+ " not allowed in curve file";
+ } else
+ throw "BUG: ref->set_coord unexpected else";
+ }
+}
+
+//***************************************************************************//
+//* CNodeRgr: register reference node *//
+//***************************************************************************//
+
+CNodeRgr::CNodeRgr( PRgr& _rgr )
+ : rgr( _rgr ) {}
+
+void CNodeRgr::tree_val( CResult& ret, const CContext& ctx ) const {
+ ret.set_r( NReg::lastresult ); }
+void CNodeRgr::set_coord( CCoord& ret, uint k ) const {
+ ret = CCoord( "$", "" ); } // primitiver gehts nicht
+
+//***************************************************************************//
+//* CNodeDummy: function argument (t) node *//
+//***************************************************************************//
+
+//***************************************************************************//
+//* CNodeCev: curve evaluation node *//
+//***************************************************************************//
+
+CNodeCev::CNodeCev( PRef& _ref, PTree _arg )
+ : ref( _ref ), arg( _arg )
+{
+ if( ref->ti )
+ throw "curve reference must not depend on index i";
+}
+
+//! A curve evaluation: fc[k,j](arg).
+
+void CNodeCev::tree_val( CResult& ret, const CContext& ctx ) const
+{
+ // get k from context:
+ uint k = ref->get_k( ctx );
+
+ // k points to a curve:
+ POlc fc = P2C( NOlm::MOM[k] );
+ if( !fc )
+ throw "file " + strg(k) + " is not a curve";
+
+ // get j from context:
+ uint j = ref->get_j( ctx, fc->nJ() );
+
+ // evaluate and pack function argument:
+ CResult farg;
+ arg->tree_val( farg, ctx );
+
+ // evaluate the curve:
+ if ( farg.vectorial ) {
+ ret.preset_v( farg.v.size() );
+ fc->curve_val_vec( &(ret.v), farg.v, k, j );
+ } else
+ ret.set_r( fc->curve_val_sca( farg.r, k, j ) );
+}
+
+void CNodeCev::set_coord( CCoord& ret, uint k ) const {
+ ret = CCoord( "c", "" ); } // primitiver gehts nicht
+
+
+//***************************************************************************//
+//* CNodeConv1: convolution or dirac, with argument: shift *//
+//***************************************************************************//
+
+void CNodeConv1::tree_val( CResult& ret, const CContext& ctx ) const
+{
+ // Shift:
+ CResult res_shift;
+ shift->tree_val( res_shift, ctx );
+ if ( res_shift.vectorial )
+ throw "shift is not a scalar";
+ double theshift = res_shift.r;
+
+ // Return values are requested for n values of _DUMMY ( = t ):
+ if ( !(ctx.vt) )
+ throw "BUG: convolution tree_val with *vt=0";
+ uint n = ctx.vt->size();
+ ret.preset_v( n );
+ ret.dv.clear();
+
+ // Resolution is given by context 'ctx':
+ uint kconv, jconv;
+ PSpec sv;
+ NCurveFile::setConv( &sv, &kconv, &jconv, ctx.k, ctx.j );
+
+ if ( kconv==-1 ) {
+ copy_theory( ret, ctx, theshift );
+ } else {
+ double conv_norm; // normalization constant (integral of resolution)
+ double conv_step; // step in x/t, if equidistant
+ string info = "resolution " + strg(kconv) + ", spectrum " + strg(jconv);
+ conv_norm = sv->norm( info );
+ conv_step = sv->has_step( ctx.vt_step ) ? ctx.vt_step : 0;
+ if( !conv_step && !allow_slow_conv )
+ throw "grid is not equidistant; slow convolution is not allowed";
+ convolve( ret, ctx, theshift, sv, conv_norm, conv_step );
+ }
+}
+
+//***************************************************************************//
+//* CNodeConv2: convolution with two arguments: theory, shift *//
+//***************************************************************************//
+
+void CNodeConv2::kdep_exec( bool *kd ) const
+{
+ theory->kdep_exec( kd );
+ shift->kdep_exec( kd );
+}
+
+void CNodeConv2::npar_exec( uint *np ) const
+{
+ theory->npar_exec( np );
+ shift->npar_exec( np );
+}
+
+
+//***************************************************************************//
+//* CNodeConv: convoluted function node *//
+//***************************************************************************//
+
+void CNodeConv::copy_theory( CResult& ret, const CContext& ctx,
+ double theshift ) const
+{
+ CResult res_theory;
+ theory->tree_val( res_theory, ctx );
+ for ( int i=0; i<ctx.vt->size(); ++i )
+ ret.v[i] = res_theory.to_r(i);
+}
+
+void CNodeConv::convolve( CResult& ret, const CContext& ctx,
+ double theshift, const PSpec& sv,
+ double conv_norm, double conv_step ) const
+{
+ uint n = ctx.vt->size();
+ uint nv = sv->size();
+ CResult res_theory;
+ CContext cv_ctx(ctx);
+ vector<double> tt; // theory grid
+ cv_ctx.vt = &tt;
+ // Convolute theory with resolution.
+ if( conv_step ){ // accelerated version with equidistant grids
+ // precompute theory on unified grid
+ int ntt = n + nv - 1;
+ tt.resize(ntt);
+ double tt0 = (*(ctx.vt))[0] - sv->x[0] - theshift - (nv-1)*conv_step;
+ for ( int ii=0; ii<ntt; ++ii ) // ii=(nv-1)+i-iv
+ tt[ii] = tt0 + ii*conv_step;
+ theory->tree_val( res_theory, cv_ctx );
+ // now convolve
+ for ( int i=0; i<n; ++i ){
+ ret.v[i] = 0;
+ for ( int iv=0; iv<nv; ++iv )
+ ret.v[i] += res_theory.to_r(nv-1+i-iv) * sv->y[iv];
+ ret.v[i] *= conv_step / conv_norm;
+ }
+ } else { // simplest version: non-equidistant x_out and x_res.
+ for ( int i=0; i<n; ++i )
+ ret.v[i] = 0;
+ tt.resize(n);
+ for ( int iv=0; iv<nv; ++iv ) {
+ double dx; // variable step
+ 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.vt))[i] - sv->x[iv] - theshift;
+ theory->tree_val( res_theory, cv_ctx );
+ for ( int i=0; i<n; ++i )
+ ret.v[i] += res_theory.to_r(i) * sv->y[iv] *
+ dx / conv_norm;
+ }
+ }
+}
+
+
+//***************************************************************************//
+//* CNodePConv: convolution with primitive *//
+//***************************************************************************//
+
+void CNodePConv::copy_theory( CResult& ret, const CContext& ctx,
+ double theshift ) const
+{
+ throw "non-convolution not implemented for pconv";
+ // to implement this efficiently, I would like to have CResult::operator+
+ // for which we should go to shared pointers ...
+}
+
+void CNodePConv::convolve( CResult& ret, const CContext& ctx,
+ double theshift, const PSpec& sv,
+ double conv_norm, double conv_step ) const
+{
+ uint n = ctx.vt->size();
+ uint nv = sv->size();
+ CContext cv_ctx(ctx);
+ vector<double> tt;
+ cv_ctx.vt = &tt;
+ // Convolute theory with resolution.
+ if( conv_step ){ // accelerated version with equidistant grids
+ // precompute primitive on unified grid
+ int ntt = n + nv;
+ tt.resize(ntt);
+ double tt0 = (*(ctx.vt))[0] - sv->x[0] - theshift - (nv-0.5)*conv_step;
+ for ( int ii=0; ii<ntt; ++ii ) // ii=(nv-1)+i-iv
+ tt[ii] = tt0 + ii*conv_step;
+ CResult res_theory;
+ theory->tree_val( res_theory, cv_ctx );
+ // now convolve resolution with differential of primitive
+ for ( int i=0; i<n; ++i ){
+ ret.v[i] = 0;
+ for ( int iv=0; iv<nv; ++iv )
+ ret.v[i] += (res_theory.to_r(nv +i-iv)-
+ res_theory.to_r(nv-1+i-iv)) * sv->y[iv];
+ ret.v[i] /= conv_norm;
+ }
+ } else { // simplest version: non-equidistant x_out and x_res.
+/* DOES NOT WORK
+ double tt0;
+ CResult res_theory;
+ tt.resize(nv);
+ if( nv<2 )
+ throw "non-equidistant convolution requires nv>=2";
+ for ( int i=0; i<n; ++i ){
+ tt0 = (*(ctx.vt))[i] - theshift;
+ for ( int iv=0; iv<nv; ++iv )
+ tt[iv] = tt0 - sv->x[iv];
+ theory->tree_val( res_theory, cv_ctx );
+ ret.v[i] = 0;
+ for ( int iv=0; iv<nv-1; ++iv ) {
+ ret.v[i] += (sv->y[iv+1]+sv->y[iv])*
+ (sv->x[iv+1]-sv->x[iv])*
+ (res_theory.to_r(iv+1) - res_theory.to_r(iv) );
+ }
+ ret.v[i] /= 2*conv_norm;
+ }
+*/
+ for ( int i=0; i<n; ++i )
+ ret.v[i] = 0;
+ tt.resize(n);
+ CResult res_theory_low, res_theory_hig;
+ double dx;
+ for ( int iv=0; iv<nv; ++iv ) {
+ 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.vt))[i] - sv->x[iv] - theshift - dx/2;
+ theory->tree_val( res_theory_low, cv_ctx );
+ }
+ for ( int i=0; i<n; ++i ) {
+ tt[i] = (*(ctx.vt))[i] - sv->x[iv] - theshift + dx/2;
+ theory->tree_val( res_theory_hig, cv_ctx );
+ }
+ for ( int i=0; i<n; ++i )
+ ret.v[i] += (res_theory_hig.to_r(i)-res_theory_low.to_r(i)) *
+ sv->y[iv] * dx;
+ }
+ for ( int i=0; i<n; ++i )
+ ret.v[i] /= conv_norm;
+ }
+}
+
+
+//***************************************************************************//
+//* CNodeDirac: node that copies the resolution function *//
+//***************************************************************************//
+
+void CNodeDirac::kdep_exec( bool *kd ) const
+{
+ shift->kdep_exec( kd );
+}
+
+void CNodeDirac::npar_exec( uint *np ) const
+{
+ shift->npar_exec( np );
+}
+
+//! Convolution or Dirac function: conv(f, shift), dirac(shift).
+
+void CNodeDirac::convolve( CResult& ret, const CContext& ctx,
+ double theshift, const PSpec& sv,
+ double conv_norm, double conv_step ) const
+{
+ uint n = ctx.vt->size();
+ // Delta function copies resolution.
+ vector<double> tt(n), yy(n);
+ // shift in t
+ for ( int i=0; i<n; ++i )
+ tt[i] = (*(ctx.vt))[i]-theshift;
+ // interpolation
+ sv->intpol( tt, yy );
+ // normalization
+ for ( int i=0; i<n; ++i )
+ ret.v[i] = yy[i] / conv_norm;
+}
+
+void CNodeDirac::copy_theory( CResult& ret, const CContext& ctx,
+ double theshift ) const
+{
+ // Naive representation of delta function.
+ for ( int i=0; i<ctx.vt->size(); ++i )
+ ret.v[i] = fabs((*(ctx.vt))[i]-theshift)<1e-2 ? 1e2 : 0;
+}
+
+void CNodeDirac::set_coord( CCoord& ret, uint k ) const {
+ ret = CCoord("resol", ""); }
diff --git a/pub/src/node.h b/pub/src/node.h
index e3665490..ae3b3a79 100644
--- a/pub/src/node.h
+++ b/pub/src/node.h
@@ -1,7 +1,7 @@
//**************************************************************************//
//* FRIDA: flexible rapid interactive data analysis *//
-//* node.h: subclasses for expression *//
-//* (C) Joachim Wuttke 2012- *//
+//* node: nodes within expression trees *//
+//* (C) Joachim Wuttke 1990-, v2(C++) 2001- *//
//* http://apps.jcns.fz-juelich.de/frida *//
//**************************************************************************//
--
GitLab