yeeah! pconv works on first attempt (equidistant case, tested with atan)

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

test/test1-fit.inp

 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