gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

PYTHON: Rewrite interpolations.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2021-10-05 10:19:18 -03:00
parent d714974b33
commit e153cd417e
1 changed files with 113 additions and 142 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

255
python/amunpy/src/amunpy/interpolation.py
View File

@ -22,14 +22,10 @@
================================================================================ ================================================================================
module: AMUN module: INTERPOLATION
Python module with subroutines to read AMUN code HDF5 files. Support module for Amun snapshots for the block prolongation with different
types of interpolation.
The only requirements for this package are:
- numpy
- scipy (optional, for data interpolation)
-------------------------------------------------------------------------------- --------------------------------------------------------------------------------
""" """
@ -43,154 +39,129 @@ except ImportError:
def rebin(a, newshape): def rebin(a, newshape):
''' '''
Subroutine changes the size of the input array to to new shape, Subroutine changes the size of the input array to to new shape,
by copying cells or averaging them. by copying cells or averaging them.
''' '''
assert len(a.shape) == len(newshape) assert len(a.shape) == len(newshape)
m = a.ndim - 1 m = a.ndim - 1
if (a.shape[m] > newshape[m]): if (a.shape[m] > newshape[m]):
if a.ndim == 3: if a.ndim == 3:
nn = [newshape[0], int(a.shape[0] / newshape[0]), nn = [newshape[0], a.shape[0] // newshape[0], \
newshape[1], int(a.shape[1] / newshape[1]), newshape[1], a.shape[1] // newshape[1], \
newshape[2], int(a.shape[2] / newshape[2])] newshape[2], a.shape[2] // newshape[2]]
return a.reshape(nn).mean(5).mean(3).mean(1) 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: else:
nn = [newshape[0], int(a.shape[0] / newshape[0]), for n in range(a.ndim):
newshape[1], int(a.shape[1] / newshape[1])] a = np.repeat(a, newshape[n] // a.shape[n], axis=n)
return a.reshape(nn).mean(3).mean(1) return(a)
else:
for n in range(a.ndim):
a = np.repeat(a, int(newshape[n] / a.shape[n]), axis=n)
return(a)
def interpolate(a, newshape, nghosts, method='rebin', order=3): def interpolate(a, newshape, nghosts=0, method=None, order=1):
''' '''
Subroutine rescales the block by interpolating its values. Subroutine rescales the block by interpolating its values.
''' '''
if (method == 'rebin' or not scipy_available): if method == None or method == 'rebin' or not scipy_available:
# calculate the indices in order to remove the ghost zones ng = nghosts
# if a.ndim == 3:
if a.ndim == 3: return rebin(a[ng:-ng,ng:-ng,ng:-ng], newshape)
ib = nghosts else:
jb = nghosts return rebin(a[ng:-ng,ng:-ng], newshape)
kb = nghosts
ie = a.shape[2] - nghosts elif method == 'zoom':
je = a.shape[1] - nghosts
ke = a.shape[0] - nghosts zf = (newshape[1] // (a.shape[1] - 2 * nghosts))
if (a.shape[0] == 1): ng = zf * nghosts
kb = 0 if a.ndim == 3:
ke = 1 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: 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)) d2 = np.roll(q,-1, axis=0) + np.roll(q, 1, axis=0) - 2.0 * q
ib, ie = fc * nghosts, - fc * nghosts q = q - d2 / 24.0
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
return zoom(a, fc, order=order, grid_mode=True, mode='nearest')[kb:ke,jb:je,ib:ie] d = np.array(q.shape)
else:
return zoom(a, fc, order=order, grid_mode=True, mode='nearest')[jb:je,ib:ie]
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 return 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