From e153cd417e207ceee081434815e97d68ee54ab87 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 5 Oct 2021 10:19:18 -0300 Subject: [PATCH] PYTHON: Rewrite interpolations. Signed-off-by: Grzegorz Kowal --- python/amunpy/src/amunpy/interpolation.py | 255 ++++++++++------------ 1 file changed, 113 insertions(+), 142 deletions(-) diff --git a/python/amunpy/src/amunpy/interpolation.py b/python/amunpy/src/amunpy/interpolation.py index 52ea89d..c1cf00d 100644 --- a/python/amunpy/src/amunpy/interpolation.py +++ b/python/amunpy/src/amunpy/interpolation.py @@ -22,14 +22,10 @@ ================================================================================ - module: AMUN + module: INTERPOLATION - Python module with subroutines to read AMUN code HDF5 files. - - The only requirements for this package are: - - - numpy - - scipy (optional, for data interpolation) + Support module for Amun snapshots for the block prolongation with different + types of interpolation. -------------------------------------------------------------------------------- """ @@ -43,154 +39,129 @@ except ImportError: def rebin(a, newshape): - ''' - Subroutine changes the size of the input array to to new shape, - by copying cells or averaging them. - ''' - assert len(a.shape) == len(newshape) + ''' + Subroutine changes the size of the input array to to new shape, + by copying cells or averaging them. + ''' + assert len(a.shape) == len(newshape) - m = a.ndim - 1 - if (a.shape[m] > newshape[m]): - if a.ndim == 3: - nn = [newshape[0], int(a.shape[0] / newshape[0]), - newshape[1], int(a.shape[1] / newshape[1]), - newshape[2], int(a.shape[2] / newshape[2])] - return a.reshape(nn).mean(5).mean(3).mean(1) + m = a.ndim - 1 + if (a.shape[m] > newshape[m]): + if a.ndim == 3: + nn = [newshape[0], a.shape[0] // newshape[0], \ + newshape[1], a.shape[1] // newshape[1], \ + newshape[2], a.shape[2] // newshape[2]] + return a.reshape(nn).mean(5).mean(3).mean(1) + else: + nn = [newshape[0], a.shape[0] // newshape[0], \ + newshape[1], a.shape[1] // newshape[1]] + return a.reshape(nn).mean(3).mean(1) else: - nn = [newshape[0], int(a.shape[0] / newshape[0]), - newshape[1], int(a.shape[1] / newshape[1])] - return a.reshape(nn).mean(3).mean(1) - else: - for n in range(a.ndim): - a = np.repeat(a, int(newshape[n] / a.shape[n]), axis=n) - return(a) + for n in range(a.ndim): + a = np.repeat(a, newshape[n] // a.shape[n], axis=n) + return(a) -def interpolate(a, newshape, nghosts, method='rebin', order=3): - ''' - Subroutine rescales the block by interpolating its values. - ''' - if (method == 'rebin' or not scipy_available): +def interpolate(a, newshape, nghosts=0, method=None, order=1): + ''' + Subroutine rescales the block by interpolating its values. + ''' + if method == None or method == 'rebin' or not scipy_available: - # calculate the indices in order to remove the ghost zones - # - if a.ndim == 3: - ib = nghosts - jb = nghosts - kb = nghosts - ie = a.shape[2] - nghosts - je = a.shape[1] - nghosts - ke = a.shape[0] - nghosts - if (a.shape[0] == 1): - kb = 0 - ke = 1 + ng = nghosts + if a.ndim == 3: + return rebin(a[ng:-ng,ng:-ng,ng:-ng], newshape) + else: + return rebin(a[ng:-ng,ng:-ng], newshape) + + elif method == 'zoom': + + zf = (newshape[1] // (a.shape[1] - 2 * nghosts)) + ng = zf * nghosts + if a.ndim == 3: + return zoom(a, zf, order=order, grid_mode=True, mode='nearest')[ng:-ng,ng:-ng,ng:-ng] + else: + return zoom(a, zf, order=order, grid_mode=True, mode='nearest')[ng:-ng,ng:-ng] + + elif method in [ 'monotonic', 'pchip' ]: + + dims = np.arange(a.ndim) + q = a + + for n in dims: + + d2 = np.roll(q,-1, axis=0) + np.roll(q, 1, axis=0) - 2.0 * q + q = q - d2 / 24.0 + + d = np.array(q.shape) + + xo = (np.arange(0.5, a.shape[n]) - nghosts) / (a.shape[n] - 2 * nghosts) + xn = np.arange(0.5, newshape[n]) / newshape[n] + + u = q.reshape([d[0], q.size // d[0]]) + f = np.zeros([newshape[n], q.size // d[0]]) + for i in range(q.size // d[0]): + f[:,i] = pchip_interpolate(xo, u[:,i], xn) + + d[0] = newshape[n] + f = f.reshape(d) + + q = f.transpose(np.roll(dims, -1)) + + return q + + elif method == 'spline': + + dims = np.arange(a.ndim) + q = a + + for n in dims: + + d2 = np.roll(q,-1, axis=0) + np.roll(q, 1, axis=0) - 2.0 * q + q = q - d2 / 24.0 + + d = np.array(q.shape) + + xo = (np.arange(0.5, a.shape[n]) - nghosts) / (a.shape[n] - 2 * nghosts) + xn = np.arange(0.5, newshape[n]) / newshape[n] + + u = q.reshape([d[0], q.size // d[0]]) + f = np.zeros([newshape[n], q.size // d[0]]) + for i in range(q.size // d[0]): + t = splrep(xo, u[:,i], k=5, s=0.0) + f[:,i] = splev(xn, t) + + d[0] = newshape[n] + f = f.reshape(d) + + q = f.transpose(np.roll(dims, -1)) + + return q - return rebin(a[kb:ke,jb:je,ib:ie], newshape) else: - ib = nghosts - jb = nghosts - ie = a.shape[1] - nghosts - je = a.shape[0] - nghosts - return rebin(a[jb:je,ib:ie], newshape) + dims = np.arange(a.ndim) + q = a - elif (method == 'zoom'): + for n in dims: - fc = int(newshape[1] / (a.shape[1] - 2 * nghosts)) - ib, ie = fc * nghosts, - fc * nghosts - jb, je = fc * nghosts, - fc * nghosts - if a.ndim == 3: - kb, ke = fc * nghosts, - fc * nghosts - if (a.shape[0] == 1): - kb = 0 - ke = 1 + d2 = np.roll(q,-1, axis=0) + np.roll(q, 1, axis=0) - 2.0 * q + q = q - d2 / 24.0 - return zoom(a, fc, order=order, grid_mode=True, mode='nearest')[kb:ke,jb:je,ib:ie] - else: - return zoom(a, fc, order=order, grid_mode=True, mode='nearest')[jb:je,ib:ie] + d = np.array(q.shape) - elif (method == 'monotonic' or method == 'pchip'): + xo = (np.arange(0.5, a.shape[n]) - nghosts) / (a.shape[n] - 2 * nghosts) + xn = np.arange(0.5, newshape[n]) / newshape[n] - dims = np.arange(a.ndim) + u = q.reshape([d[0], q.size // d[0]]) + f = np.zeros([newshape[n], q.size // d[0]]) + for i in range(q.size // d[0]): + t = interp1d(xo, u[:,i], kind=method) + f[:,i] = t(xn) - q = a + d[0] = newshape[n] + f = f.reshape(d) - for n in dims: + q = f.transpose(np.roll(dims, -1)) - d2 = np.roll(q,-1, axis=0) + np.roll(q, 1, axis=0) - 2.0 * q - q = q - d2 / 24.0 - - d = np.array(q.shape) - - xo = (np.arange(0.5, a.shape[n]) - nghosts) / (a.shape[n] - 2 * nghosts) - xn = np.arange(0.5, newshape[n]) / newshape[n] - - u = q.reshape([d[0], int(q.size / d[0])]) - f = np.zeros([newshape[n], int(q.size / d[0])]) - for i in range(int(q.size / d[0])): - f[:,i] = pchip_interpolate(xo, u[:,i], xn) - - d[0] = newshape[n] - f = f.reshape(d) - - q = f.transpose(np.roll(dims, -1)) - - return q - - elif (method == 'spline'): - - dims = np.arange(a.ndim) - - q = a - - for n in dims: - - d2 = np.roll(q,-1, axis=0) + np.roll(q, 1, axis=0) - 2.0 * q - q = q - d2 / 24.0 - - d = np.array(q.shape) - - xo = (np.arange(0.5, a.shape[n]) - nghosts) / (a.shape[n] - 2 * nghosts) - xn = np.arange(0.5, newshape[n]) / newshape[n] - - u = q.reshape([d[0], int(q.size / d[0])]) - f = np.zeros([newshape[n], int(q.size / d[0])]) - for i in range(int(q.size / d[0])): - t = splrep(xo, u[:,i], k=5, s=0.0) - f[:,i] = splev(xn, t) - - d[0] = newshape[n] - f = f.reshape(d) - - q = f.transpose(np.roll(dims, -1)) - - return q - - else: - - dims = np.arange(a.ndim) - - q = a - - for n in dims: - - d2 = np.roll(q,-1, axis=0) + np.roll(q, 1, axis=0) - 2.0 * q - q = q - d2 / 24.0 - - d = np.array(q.shape) - - xo = (np.arange(0.5, a.shape[n]) - nghosts) / (a.shape[n] - 2 * nghosts) - xn = np.arange(0.5, newshape[n]) / newshape[n] - - u = q.reshape([d[0], int(q.size / d[0])]) - f = np.zeros([newshape[n], int(q.size / d[0])]) - for i in range(int(q.size / d[0])): - t = interp1d(xo, u[:,i], kind=method) - f[:,i] = t(xn) - - d[0] = newshape[n] - f = f.reshape(d) - - q = f.transpose(np.roll(dims, -1)) - - return q + return q