From 66dd58a94b2da65eebf8d54734af406966a15351 Mon Sep 17 00:00:00 2001 From: Marco Cammarata Date: Thu, 1 Dec 2016 11:15:17 +0100 Subject: [PATCH] trying to get the gauss convolution to work again... --- mcutils/mcutils.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/mcutils/mcutils.py b/mcutils/mcutils.py index 58db3fc..0c10db7 100644 --- a/mcutils/mcutils.py +++ b/mcutils/mcutils.py @@ -288,7 +288,7 @@ def convolveQuad(x,y,xres,yres,useAutoCrop=True, x_integral = xresInv+x[i] ytemp = interpolate(x,y,x_integral,fill_value=fill_value,axis=axis)/area # do integration - res = integrate.simps(ytemp*yresInv,x=xresInv+x[i],axis=axis) + res = integrate.simps(ytemp[0]*yresInv,x=xresInv+x[i],axis=axis) colon = (slice(None),) setElement(out,i,res,axis=axis) return out @@ -326,6 +326,14 @@ def convolve(x,y,xres,yres,useAutoCrop=True,approximantOrder=None,fill_value=0., ## prepare output out = np.empty_like(y) + + ## fill up NaN + for i in range(y.shape[0]): + isok = np.isfinite(y[i]) + xp = x[isok] + fp = y[i,isok] + y[i] = np.interp(x, xp, fp) + nP = len(x) f = interpolator(x,y,fill_value=fill_value,axis=axis) for i in range(nP): @@ -335,10 +343,10 @@ def convolve(x,y,xres,yres,useAutoCrop=True,approximantOrder=None,fill_value=0., #ytemp = interpolate(x,y,x_integral,fill_value=fill_value,axis=axis)/area #ytemp = f(x_integral)/area # do integration - res = integrate.simps( np.transpose(ytemp.T*yresInv),x=xresInv+x[i],axis=axis) + res = integrate.simps( ytemp*yresInv,x=xresInv+x[i],axis=axis) colon = (slice(None),) setElement(out,i,res,axis=axis) - return out + return out.T def fftconvolve(x,y,yres,xres=None,normalize=False):