From bc86645f860a0604209fc5b56caae1233c092cf4 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (o)" <j.wuttke@fz-juelich.de>
Date: Tue, 2 Dec 2014 19:23:39 +0100
Subject: [PATCH] quadratic error propagation in convolution
---
pub/lib/node.cpp | 24 ++++++++++++++++--------
1 file changed, 16 insertions(+), 8 deletions(-)
diff --git a/pub/lib/node.cpp b/pub/lib/node.cpp
index 9c897729..9654ca37 100644
--- a/pub/lib/node.cpp
+++ b/pub/lib/node.cpp
@@ -662,11 +662,11 @@ void CNodeConv::convolve( CResult& ret, const CContext& ctx,
for ( int iv=0; iv<nv; ++iv ) {
ret.v[i] += res_theory.to_r(nv-1+i-iv) * sv->y[iv];
if( with_error )
- ret.dv[i] += res_theory.to_r(nv-1+i-iv) * sv->dy[iv];
+ ret.dv[i] += res_theory.to_r(nv-1+i-iv) * SQR(sv->dy[iv]);
}
ret.v[i] *= conv_step / conv_norm;
if( with_error )
- ret.dv[i] *= conv_step / conv_norm;
+ ret.dv[i] = sqrt(ret.dv[i]) * conv_step / conv_norm;
}
} else { // simplest version: non-equidistant x_out and x_res.
ret.v.assign( n, 0. );
@@ -685,11 +685,16 @@ void CNodeConv::convolve( CResult& ret, const CContext& ctx,
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;
+ ret.v[i] += res_theory.to_r(i) * sv->y[iv] * dx;
if( with_error )
- ret.dv[i] += res_theory.to_r(i) * sv->dy[iv] *dx/conv_norm;
+ ret.dv[i] += res_theory.to_r(i) * SQR(sv->dy[iv]);
}
}
+ for ( int i=0; i<n; ++i ) {
+ ret.v[i] /= conv_norm;
+ if( with_error )
+ ret.dv[i] = sqrt(ret.dv[i])/conv_norm;
+ }
}
}
@@ -743,11 +748,11 @@ void CNodePConv::convolve( CResult& ret, const CContext& ctx,
res_theory.to_r(nv-1+i-iv);
ret.v[i] += igral * sv->y[iv];
if( with_error )
- ret.dv[i] += igral * sv->dy[iv];
+ ret.dv[i] += igral * SQR(sv->dy[iv]);
}
ret.v[i] /= conv_norm;
if( with_error )
- ret.dv[i] /= conv_norm;
+ ret.dv[i] = sqrt(ret.dv[i]) / conv_norm;
}
} else { // simplest version: non-equidistant x_out and x_res.
ret.v.assign( n, 0. );
@@ -775,11 +780,14 @@ void CNodePConv::convolve( CResult& ret, const CContext& ctx,
double igral = res_theory_hig.to_r(i)-res_theory_low.to_r(i);
ret.v[i] += igral * sv->y[iv] * dx;
if ( with_error )
- ret.dv[i] += igral * sv->dy[iv] * dx;
+ ret.dv[i] += igral * SQR(sv->dy[iv]) * dx;
}
}
- for ( int i=0; i<n; ++i )
+ for ( int i=0; i<n; ++i ) {
ret.v[i] /= conv_norm;
+ if ( with_error )
+ ret.dv[i] = sqrt(ret.dv[i]) / conv_norm;
+ }
}
}
--
GitLab