quadratic error propagation in convolution

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;
+        }
     }
 }