python - 一个 "ij"meshgrid和一个长meshgrid的对应关系

Question

考虑一个矩阵 Z包含 z = z(a,m,e) 的基于网格的结果. Z有形状 (len(aGrid), len(mGrid), len(eGrid)) . Z[0,1,2]包含 z(a=aGrid[0], m=mGrid[1], e=eGrid[2]) .然而，我们可能已经从对象的状态空间中删除了一些元素(例如和简单，(a,m,e : a > 3)。假设有效状态空间的大小是x。

我have been suggested将此对象转换为对象的代码 Z2形状(x, 3) . Z2 中的每一行对应一个元素i来自 Z2 : (aGrid[a[i]], mGrid[m[i]], eGrid[e[i]]) .

# first create Z, a mesh grid based matrix that has some invalid states (we set them to NaN)
aGrid = np.arange(0, 10, dtype=float)
mGrid = np.arange(100, 110, dtype=float)
eGrid = np.arange(1000, 1200, dtype=float)
A,M,E = np.meshgrid(aGrid, mGrid, eGrid, indexing='ij')
Z = A
Z[Z > 3] = np.NaN #remove some states from being "allowed"

# now, translate them from shape (len(aGrid), len(mGrid), len(eGrid)) to 
grids = [A,M,E]
grid_bc = np.broadcast_arrays(*grids)
Z2 = np.column_stack([g.ravel() for g in grid_bc])
Z2[np.isnan(Z.ravel())] = np.nan
Z3 = Z2[~np.isnan(Z2)]

通过一些计算，我得到一个矩阵 V4形状为 Z3但包含 4 列。

我得到了

• Z2 (如上)
• Z3 (如上)
• V4这是一个矩阵形状 ( Z3.shape[0], Z3.shape[1]+1 ):它附加了一个额外的列
• (如有必要，我仍然可以访问网格 A,M,E )

我需要重新创建

• V , 这是包含 V4 的值(最后一列)的矩阵, 但会变回 Z1 的形状.

也就是说，如果V4中有一行上面写着 (aGrid[0], mGrid[1], eGrid[2], v1) , 然后是 V 的值在 V[0,1,2] = v1等 V4 中的所有行,

效率是关键。

Answer 1

给定您的原始问题条件，按如下方式重新创建，并修改为 A 是 Z 的副本:

aGrid = np.arange(0, 10, dtype=float)
mGrid = np.arange(100, 110, dtype=float)
eGrid = np.arange(1000, 1200, dtype=float)
A,M,E = np.meshgrid(aGrid, mGrid, eGrid, indexing='ij')
Z = A.copy()
Z[Z > 3] = np.NaN 

grids = [A,M,E]
grid_bc = np.broadcast_arrays(*grids)
Z2 = np.column_stack([g.ravel() for g in grid_bc])
Z2[np.isnan(Z.ravel())] = np.nan
Z3 = Z2[~np.isnan(Z2)]

函数可以定义如下，从稀疏的二维 # 数据点 x # dims + 1 矩阵重新创建密集的 N 维矩阵。该函数的第一个参数是上述二维矩阵，最后一个(可选)参数是每个维度的网格索引:

import numpy as np
def map_array_to_index(uniq_arr):
    return np.vectorize(dict(map(reversed, enumerate(uniq_arr))).__getitem__)

def recreate(arr, *coord_arrays):
    if len(coord_arrays) != arr.shape[1] - 1:
      coord_arrays = map(np.unique, arr.T[0:-1])
    lookups = map(map_array_to_index, coord_arrays)
    new_array = np.nan * np.ones(map(len, coord_arrays))
    new_array[tuple(l(c) for c, l in zip(arr.T[0:-1], lookups))] = arr[:, -1]
    new_grids = np.meshgrid(*coord_arrays, indexing='ij')
    return new_array, new_grids

给定一个二维矩阵 V4，上面用从 Z 派生的值定义，

V4 = np.column_stack([g.ravel() for g in grid_bc] + [Z.ravel()])

可以按如下方式重新创建 Z:

V4_orig_form, V4_grids = recreate(V4, aGrid, mGrid, eGrid)

所有非 NaN 值都正确测试是否相等:

np.all(Z[~np.isnan(Z)] == V4_orig_form[~np.isnan(V4_orig_form)])

该函数在没有传入 aGrid、mGrid、eGrid 的情况下也可以工作，但在这种情况下，它不会包含输入数组相应列中不存在的任何坐标。