在Python中在非均匀网格上移动数组



我想知道Numpy或SciPy中是否有允许在非均匀网格上移动数组的Python功能。我创建了一个最小的例子来说明这个过程,但这似乎在这个最小的例子中不起作用:

import numpy as np
import matplotlib.pyplot as pyt
def roll_arrays( a, shift_values,x_grid ):
#from scipy.interpolate import interp1d

x_max         = np.amax(x_grid)    
total_items   = a.shape[0]  
the_ddtype    = a.dtype 

result = np.zeros( (a.shape[0], a.shape[1] ), dtype=the_ddtype )    


for k in range( total_items ):        
edge_val_left  = a[k,0]
edge_val_right = a[k,-1]     

#extend grid to edges with boundary values (flat extrapolation)
extended_boundary = np.abs( shift_values[k] )#positive or negative depending on shift 

if( shift_values[k] != 0.0 ):

x0_right          = np.linspace( x_max +1e-3, x_max + 1e-3 + extended_boundary, 10 )
x0_left           = np.linspace( -x_max - 1e-3 -extended_boundary, -x_max - 1e-3, 10 )
if( shift_values[k]>0.0 ):
#we fill left values
x_dense_grid  = np.concatenate( ( x0_left, x_grid + shift_values[k] ) ) 
ynew          = np.concatenate(  ( edge_val_left*np.ones( 10 ), a[k,:] )  )            

elif( shift_values[k]<0.0 ):
x_dense_grid  = np.concatenate( ( x_grid + shift_values[k], x0_right ) )               
ynew          = np.concatenate(  ( a[k,:], edge_val_right*np.ones( 10 ) )  ) 


###
#return on the original grid                       
f_interp    = np.interp( x_grid, x_dense_grid, ynew )                
result[k,:] = f_interp  

else:
#no shift
result[k,:] = a[k,:]


return result

x_geom     = np.array( [ 100*( 1.5**(-0.5*k) ) for k in range(1000)] )
x_geom_neg =-( x_geom )
x_geom = np.concatenate( (np.array([0.0]), np.flip(x_geom)) )
x_geom = np.concatenate( (x_geom_neg, x_geom) )
shifts = np.array([-1.0,-2.0,1.0])
f      = np.array( [ k**2/( x_geom**2 + k**4 ) for k in range(1,shifts.shape[0]+1)  ] )
fs     = roll_arrays( f, shifts, x_geom)
pyt.plot( x_geom, f[0,:], marker='.' )
pyt.plot( x_geom, fs[0,:], marker='.' )

print("done")

注意;x_ grid";在这种情况下,是对数间隔的。有没有办法利用Scipy/Numpy做到这一点?通过插值方法或类似方法。

EDIT:我注意到,删除关于边界偏移的if,elif,else语句(在进行平面外推的情况下(似乎可以解决这个问题;但我仍然认为,对于Python中已经存在的东西来说,这是一个过于天真的实现;因此问题仍然存在。

如果我正确理解这个问题,np.interp将按照你的意愿执行(默认情况下,它会复制边缘的值(:

def roll_arrays(a, shift_values, x_grid):
total_items = a.shape[0]
result = np.zeros_like(a)
for k in range(total_items):
if shift_values[k] != 0.0:
# shift the x values
x_grid_shifted = x_grid + shift_values[k]
# interpolate back to the original grid
f_interp = np.interp(x_grid, x_grid_shifted, a[k, :])
result[k, :] = f_interp
else:
# no shift
result[k, :] = a[k, :]
return result

对于问题的示例输入,这将给出非常接近的内容

fs_expected = np.array([k ** 2 / ((x_geom - shift) ** 2 + k ** 4) for k, shift in enumerate(shifts, start=1)])

最新更新