From 4606ebdbb3b5afd5254e9734ef237fbe1e7b75e3 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (laptop)" <j.wuttke@fz-juelich.de>
Date: Mon, 24 Sep 2012 23:06:16 +0200
Subject: [PATCH] yeeah! pconv works on first attempt (equidistant case, tested
with atan)
---
pub/src/expr.cpp | 50 +++++++++++++++++++++++++++++++++-------------
test/test1-fit.inp | 2 +-
2 files changed, 37 insertions(+), 15 deletions(-)
diff --git a/pub/src/expr.cpp b/pub/src/expr.cpp
index d883957c..eac26d73 100644
--- a/pub/src/expr.cpp
+++ b/pub/src/expr.cpp
@@ -940,20 +940,42 @@ void CNodePConv::convolve( CResult& ret, const CContext& ctx,
vector<double> tt;
cv_ctx.vt = &tt;
// Convolute theory with resolution.
- if( !conv_step )
- throw "pconv not implemented for non-equidistant grid";
- // both grids have step conv_step
- int ntt = n + nv - 1;
- tt.resize(ntt);
- for ( int ii=0; ii<ntt; ++ii ) // ii=(nv-1)+i-iv
- tt[ii] = (*(ctx.vt))[0] - sv->x[0] - theshift +
- (ii-(nv-1))*conv_step;
- theory->tree_val( res_theory, cv_ctx );
- 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;
+ if( conv_step ){ // accelerated version with equidistant grids
+ // precompute theory 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;
+ 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 +i-iv)-
+ 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.
+ throw "pconv non-eq conv not yet implemented";
+ 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;
+ }
}
}
diff --git a/test/test1-fit.inp b/test/test1-fit.inp
index a1f9bce2..9e5e7abb 100644
--- a/test/test1-fit.inp
+++ b/test/test1-fit.inp
@@ -1,7 +1,7 @@
fl dat-res
fl dat-qel
-1 cc p0+p1*resol(0)+p2*conv(cauchy(t,p3),0)
1 cc p0+p1*resol+p2*conv(cauchy(t,p3))
+1 cc p0+p1*resol+p2*pconv(atan(p3*t))
2,3 cv 0
2,3 cf
g2
--
GitLab