2017-03-26 62 views
0

我有以下问题:迭代2D网格内插

我有两个网格,s0Bs1B和,两者20x20,并且在这些网格点评价函数s_0(s0, s1)

但是,对于某些角落,无法评估s_0(s0, s1),并且返回np.nan

我想有一个插值函数我评估的网格点,然后我可以在根解决问题中使用(固定点s_0和模拟s_1)。下面是代码的草图现在:

# mesh 1d grids: 
ss0 = np.array([ 0.38445455, 0.40388198, 0.42110923, 0.43626859, 0.44949237, 
    0.46091287, 0.47066239, 0.47887323, 0.48567771, 0.49120811, 
    0.49559675, 0.49897593, 0.50147794, 0.50323509, 0.50437969, 
    0.50504404, 0.50536044, 0.50546118, 0.50547859, 0.50554495]) 
ss1 = np.array([ 0.43490909, 0.50764658, 0.5719651 , 0.62838234, 0.67741596, 
    0.71958365, 0.75540309, 0.78539195, 0.81006791, 0.82994865, 
    0.84555185, 0.85739519, 0.86599635, 0.87187299, 0.87554281, 
    0.87752348, 0.87833267, 0.87848808, 0.87850736, 0.87890821]) 

s0B, s1B = np.meshgrid(ss0, ss1, indexing='ij') 
# generate interpolated functions 
idx = ~np.isnan(s0Max) 
if idx.sum() == 0: 
    # check if there are any points at all 
    raise notImplementedError 
s0MaxInterp = interpolate.interp2d(
    s0B[idx], s1B[idx], s0Max[idx], 
    fill_value=np.nan, kind='linear') 
idx = ~np.isnan(s1Max) 
s1MaxInterp = interpolate.interp2d(
    s0B[idx], s1B[idx], s1Max[idx], 
    fill_value=np.nan, kind='linear') 

def errorFixedPoint(x, s0MaxInterp, s1MaxInterp, grids): 
    s0, s1 = x 
    # I skip some checks whether s0, s1 are inside the interpolation range 
    return np.array([s0 - s0MaxInterp(s0, s1)[0], s1 - s1MaxInterp(s0, s1)[0]]) 

result = optimize.root(errorFixedPoint, np.array([0.5, 0.9]), 
          args=(s0MaxInterp, s1MaxInterp, (ss0, ss1)), method='lm') 

然而,这种尝试的时候,我得到一个RunTimeWarning

RuntimeWarning: No more knots can be added because the number of B-spline 
coefficients already exceeds the number of data points m. 
Probable causes: either s or m too small. (fp>s) 
    kx,ky=1,1 nx,ny=20,20 m=314 fp=0.000001 s=0.000000 
    warnings.warn(RuntimeWarning(_iermess2[ierm][0] + _mess)) 

插值有时显得相当没谱:

s0Max, grid and interpolated

但是,这可能是因为该函数在两个维度上都是非常非线性的:

s0Max in 2d

我读上,而不是interpolate2d建议scipy.interpolate.griddata等地。但是,我必须迭代插入(对于根查找过程),并且网格数据不支持。

我该如何解决这个问题?

对于重现,我贴的s0Maxs1Max这里的值:

nan = np.nan 
s0Max = np.array([[  nan,   nan,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan, 
     0.51494141, 0.51347656, 0.51210937, 0.51210937, 0.51210937, 
     0.51210937, 0.51210937, 0.51210937, 0.51210937, 0.51210937], 
     [  nan,   nan,   nan,   nan,   nan, 
     0.51337891, 0.5109375 , 0.51083984, 0.50253906, 0.50175781, 
     0.50429687, 0.50957031, 0.50859375, 0.50849609, 0.50839844, 
     0.50839844, 0.50839844, 0.50839844, 0.50839844, 0.50839844], 
     [  nan,   nan,   nan, 0.51347656, 0.50527344, 
     0.50009766, 0.49824219, 0.49746094, 0.49746094, 0.496875 , 
     0.49746094, 0.49941406, 0.50742187, 0.50732422, 0.50732422, 
     0.50732422, 0.50732422, 0.50732422, 0.50732422, 0.50732422], 
     [  nan,   nan, 0.51347656, 0.5015625 , 0.49863281, 
     0.49746094, 0.49628906, 0.49619141, 0.49570313, 0.49570313, 
     0.49619141, 0.49746094, 0.50722656, 0.50634766, 0.50625 , 
     0.50625 , 0.50625 , 0.50625 , 0.50625 , 0.50625 ], 
     [  nan,   nan, 0.50253906, 0.49824219, 0.496875 , 
     0.49619141, 0.49560547, 0.49511719, 0.49511719, 0.49511719, 
     0.49550781, 0.49677734, 0.50625 , 0.50625 , 0.50625 , 
     0.50625 , 0.50625 , 0.50625 , 0.50625 , 0.50625 ], 
     [  nan, 0.51210937, 0.49941406, 0.49746094, 0.49619141, 
     0.49560547, 0.49511719, 0.49453125, 0.49453125, 0.49453125, 
     0.49511719, 0.49628906, 0.50625 , 0.50527344, 0.50527344, 
     0.50527344, 0.50527344, 0.50527344, 0.50527344, 0.50527344], 
     [  nan, 0.50341797, 0.49804688, 0.49628906, 0.49570313, 
     0.49511719, 0.49453125, 0.49414062, 0.49404297, 0.49404297, 
     0.49453125, 0.49570313, 0.50527344, 0.50527344, 0.50527344, 
     0.50527344, 0.50527344, 0.50527344, 0.50527344, 0.50527344], 
     [  nan, 0.50078125, 0.49746094, 0.49619141, 0.49511719, 
     0.49453125, 0.49414062, 0.49404297, 0.49394531, 0.49394531, 
     0.49453125, 0.49570313, 0.50527344, 0.50527344, 0.50527344, 
     0.50527344, 0.50527344, 0.50527344, 0.50527344, 0.50527344], 
     [  nan, 0.49941406, 0.496875 , 0.49570313, 0.49511719, 
     0.49453125, 0.49404297, 0.49394531, 0.49355469, 0.49394531, 
     0.49414062, 0.49570313, 0.50527344, 0.50527344, 0.50517578, 
     0.50517578, 0.50507812, 0.50507812, 0.50507812, 0.50449219], 
     [  nan, 0.49873047, 0.49667969, 0.49560547, 0.49462891, 
     0.49414062, 0.49394531, 0.49355469, 0.49345703, 0.49355469, 
     0.49453125, 0.49628906, 0.50527344, 0.50517578, 0.50429687, 
     0.50429687, 0.50429687, 0.50429687, 0.50429687, 0.50429687], 
     [ 0.51347656, 0.49863281, 0.49628906, 0.49511719, 0.49453125, 
     0.49404297, 0.49394531, 0.49355469, 0.49355469, 0.49404297, 
     0.49570313, 0.50527344, 0.50429687, 0.50429687, 0.50429687, 
     0.50429687, 0.50429687, 0.50429687, 0.50429687, 0.50429687], 
     [ 0.51347656, 0.49814453, 0.49628906, 0.49511719, 0.49453125, 
     0.49404297, 0.49394531, 0.49355469, 0.49404297, 0.49570313, 
     0.50527344, 0.50429687, 0.50429687, 0.50429687, 0.50351563, 
     0.50351563, 0.50341797, 0.50341797, 0.50341797, 0.50341797], 
     [ 0.51347656, 0.49814453, 0.49619141, 0.49511719, 0.49453125, 
     0.49404297, 0.49394531, 0.49404297, 0.49511719, 0.50527344, 
     0.50429687, 0.50429687, 0.50351563, 0.50341797, 0.50341797, 
     0.50332031,   nan,   nan,   nan,   nan], 
     [ 0.51347656, 0.49814453, 0.49619141, 0.49511719, 0.49453125, 
     0.49404297, 0.49404297, 0.49453125, 0.49873047, 0.50507812, 
     0.50429687, 0.50351563, 0.50341797,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan], 
     [ 0.51347656, 0.49814453, 0.49619141, 0.49511719, 0.49453125, 
     0.49404297, 0.49404297, 0.49511719, 0.50527344, 0.50429687, 
     0.50429687, 0.50341797,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan], 
     [ 0.51347656, 0.49804688, 0.49619141, 0.49511719, 0.49453125, 
     0.49404297, 0.49453125, 0.49570313, 0.50527344, 0.50429687, 
     0.50429687, 0.50341797,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan], 
     [ 0.51347656, 0.49804688, 0.49619141, 0.49511719, 0.49453125, 
     0.49414062, 0.49453125, 0.49619141, 0.50527344, 0.50429687, 
     0.50429687, 0.50341797,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan], 
     [ 0.51347656, 0.49804688, 0.49619141, 0.49511719, 0.49453125, 
     0.49414062, 0.49453125, 0.49628906, 0.50527344, 0.50429687, 
     0.50351563, 0.50341797,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan], 
     [ 0.51347656, 0.49804688, 0.49619141, 0.49511719, 0.49453125, 
     0.49414062, 0.49453125, 0.49628906, 0.50527344, 0.50429687, 
     0.50351563, 0.50341797,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan], 
     [ 0.51347656, 0.49804688, 0.49619141, 0.49511719, 0.49453125, 
     0.49414062, 0.49453125, 0.49628906, 0.50527344, 0.50429687, 
     0.50351563, 0.50332031,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan]]) 

s1Max = np.array([[  nan,   nan,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan, 
     0.99894531, 0.99895996, 0.99895996, 0.99895996, 0.99895996, 
     0.99895996, 0.99895996, 0.99895996, 0.99895996, 0.99895996], 
     [  nan,   nan,   nan,   nan,   nan, 
     0.99895996, 0.99895996, 0.99895996, 0.98726563, 0.98429688, 
     0.99896118, 0.99895996, 0.99895996, 0.99895996, 0.99895996, 
     0.99895996, 0.99895996, 0.99895996, 0.99895996, 0.99895996], 
     [  nan,   nan,   nan, 0.99895996, 0.99896118, 
     0.97  , 0.95265625, 0.93984375, 0.93195313, 0.930625 , 
     0.93929688, 0.96421875, 0.99895996, 0.99895996, 0.99895996, 
     0.99895996, 0.99895996, 0.99895996, 0.99895996, 0.99895996], 
     [  nan,   nan, 0.99895996, 0.98117188, 0.95429688, 
     0.93351563, 0.91695313, 0.90453125, 0.89703125, 0.896875 , 
     0.908125 , 0.93789063, 0.99895996, 0.99895996, 0.99895996, 
     0.99895996, 0.99895996, 0.99895996, 0.99895996, 0.99895996], 
     [  nan,   nan, 0.98648438, 0.95296875, 0.9278125 , 
     0.9075 , 0.89101563, 0.87859375, 0.8715625 , 0.8725 , 
     0.88617188, 0.92085938, 0.99895996, 0.99895996, 0.99895996, 
     0.99895996, 0.99895996, 0.99895996, 0.99895996, 0.99895996], 
     [  nan, 0.99895996, 0.9640625 , 0.93265625, 0.90789063, 
     0.88757813, 0.87109375, 0.85875 , 0.85210938, 0.8540625 , 
     0.86992188, 0.90914063, 0.99895996, 0.99895996, 0.99895996, 
     0.99895996, 0.99895996, 0.99895996, 0.99895996, 0.99895996], 
     [  nan, 0.99023438, 0.94773438, 0.91703125, 0.89234375, 
     0.87203125, 0.85546875, 0.84320313, 0.83695313, 0.83976563, 
     0.85757813, 0.9009375 , 0.99895996, 0.99895996, 0.99895996, 
     0.99895996, 0.99895996, 0.99895996, 0.99895996, 0.99895996], 
     [  nan, 0.97570313, 0.93515625, 0.90476563, 0.88015625, 
     0.8596875 , 0.843125 , 0.8309375 , 0.825  , 0.82867188, 
     0.84820313, 0.8953125 , 0.99895996, 0.99895996, 0.99895996, 
     0.99895996, 0.99895996, 0.99895996, 0.99895996, 0.99895996], 
     [  nan, 0.96515625, 0.92546875, 0.89515625, 0.87046875, 
     0.85  , 0.83335938, 0.82140625, 0.81585938, 0.82046875, 
     0.84203125, 0.89390625, 0.99895996, 0.99895996, 0.99895996, 
     0.99895996, 0.99895996, 0.99895996, 0.99895996, 0.99895996], 
     [  nan, 0.95757813, 0.918125 , 0.88773438, 0.86296875, 
     0.8425 , 0.82609375, 0.8146875 , 0.81054688, 0.81820313, 
     0.84664063, 0.91648438, 0.99895996, 0.99895996, 0.99895996, 
     0.99895996, 0.99895996, 0.99895996, 0.99895996, 0.99895996], 
     [ 0.99896484, 0.953125 , 0.913125 , 0.88242188, 0.8575 , 
     0.83726563, 0.821875 , 0.81289062, 0.81414063, 0.83429688, 
     0.89601563, 0.99895996, 0.99895996, 0.99895996, 0.99895996, 
     0.99895996, 0.99895996, 0.99895996, 0.99895996, 0.99895996], 
     [ 0.99896484, 0.95078125, 0.90992188, 0.87867188, 0.85390625, 
     0.83453125, 0.82164063, 0.81859375, 0.83421875, 0.89546875, 
     0.99895996, 0.99895996, 0.99895996, 0.99896484, 0.99125 , 
     0.9909375 , 0.99078125, 0.99078125, 0.99078125, 0.99070313], 
     [ 0.99896484, 0.94945313, 0.9078125 , 0.87632813, 0.85195313, 
     0.83429688, 0.82554688, 0.83273438, 0.87804688, 0.99895996, 
     0.99895996, 0.99895996, 0.99117188, 0.99023438, 0.98945313, 
     0.98890625,   nan,   nan,   nan,   nan], 
     [ 0.99896484, 0.94859375, 0.90640625, 0.875  , 0.85125 , 
     0.835625 , 0.83203125, 0.8534375 , 0.9571875 , 0.99895996, 
     0.99895996, 0.99125 , 0.98992188,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan], 
     [ 0.99896484, 0.948125 , 0.905625 , 0.87429688, 0.85125 , 
     0.83757813, 0.83914063, 0.87632813, 0.99895996, 0.99895996, 
     0.99896484, 0.99046875,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan], 
     [ 0.99896484, 0.94789063, 0.90523438, 0.87398438, 0.85148438, 
     0.83929688, 0.8446875 , 0.89546875, 0.99895996, 0.99895996, 
     0.99896484, 0.98992188,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan], 
     [ 0.99896484, 0.94773438, 0.905  , 0.87382813, 0.85164063, 
     0.84023438, 0.84789063, 0.906875 , 0.99895996, 0.99895996, 
     0.99140625, 0.98960938,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan], 
     [ 0.99896484, 0.94773438, 0.90492188, 0.87382813, 0.85171875, 
     0.84054688, 0.84890625, 0.9109375 , 0.99895996, 0.99895996, 
     0.99132813, 0.98953125,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan], 
     [ 0.99896484, 0.94773438, 0.90492188, 0.87375 , 0.85171875, 
     0.840625 , 0.84914063, 0.91164063, 0.99895996, 0.99895996, 
     0.99132813, 0.98945313,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan], 
     [ 0.99896484, 0.94765625, 0.90492188, 0.87375 , 0.85179688, 
     0.84085938, 0.84984375, 0.91445313, 0.99895996, 0.99895996, 
     0.99132813, 0.989375 ,   nan,   nan,   nan, 
       nan,   nan,   nan,   nan,   nan]]) 

回答

0

根据这个优秀的答案here,使用interp2d不会为更复杂的数据。建议使用radial basis function interpolation

这似乎工作,并提供更好的插值结果,

import scipy.optimize as optimize 
from scipy.interpolate import Rbf 

s0MaxInterp = Rbf(s0B[idx], s1B[idx], s0Max[idx], function='cubic') 
s1MaxInterp = Rbf(s0B[idx], s1B[idx], s1Max[idx], function='cubic') 

def errorFixedPoint(x, s0MaxInterp, s1MaxInterp, grids): 
    s0, s1 = x 
    # I skip some checks whether s0, s1 are inside the interpolation range 
    return np.array([s0 - s0MaxInterp(s0, s1), s1 - s1MaxInterp(s0, s1)]) 

result = optimize.root(errorFixedPoint, np.array([0.5, 0.9]), 
         args=(s0MaxInterp, s1MaxInterp, (ss0, ss1)), method='lm') 

注意,我在errorFixedPoint功能删除[0]