在Matlab的delsq演示中，用numpy方法对向量映射进行ndarray处理？

Question

在Matlab中，你可以很容易地设置一个带编号的网格，并通过各种索引技巧将操作符映射到“矢量化”的表单上。更明确地编写这个并不难，但我想知道是否有一个numpy方法来实现这一点？也许是vianumpy.nditer？在

这段代码（参见mathworks示例）使用存储在矩阵网格G中的顺序从向量u到矩阵u：

U = G;
U(G>0) = full(u(G(G>0)));

参考文献：

http://www.mathworks.com/products/matlab/examples.html?file=/products/demos/shipping/matlab/delsqdemo.html

Answer 1

我不知道你所说的“标签网格”或“编号网格”是什么意思。在

看起来演示依赖于两个函数，numgrid和{}。你可以在MATLAB中查看代码。但在网上我发现

http://chmielowski.eu/POLITECHNIKA/Dydaktyka/AUTOMATYKA/AutoLab/Matlab/TOOLBOX/MATLAB/DEMOS/NUMGRID.M

http://chmielowski.eu/POLITECHNIKA/Dydaktyka/AUTOMATYKA/AutoLab/Matlab/TOOLBOX/MATLAB/DEMOS/DELSQ.M

日期是1991年，所以它们是用老式的MATLAB编写的，没有更新的“矢量化”或任何花哨的映射技巧（我可以看到）。唯一特别的是使用稀疏矩阵（delsq返回稀疏矩阵）。它生成的曲面已经是Mathworks标志很长时间了。在

所以将这些转换成{}应该是直接的。要生成稀疏矩阵，您需要scipy.sparse。MATLABsparse(i,j,s)=>；scipy.sparse.csr_matrix(s, (i,j))。在

但如果您更关心索引步骤：

U = G;
U(G>0) = full(u(G(G>0)));

对于匹配的data数组，以下内容是等效的

倍频程（MATLAB）：

^{pr2}$

numpy

U=np.zeros(data.shape)
u = np.arange(0,76)*.1+10
U[data>0]=u[data[data>0]]
contour(U)

Answer 2

下面复制了演示。G中的编号是不同的，但数字只是标签（标签网格让我困惑）。在

import numpy as np
from scipy import sparse
from scipy.sparse import linalg
import matplotlib.pyplot as plt

def numgrid(n):
    """
    NUMGRID Number the grid points in a two dimensional region.
    G = NUMGRID('R',n) numbers the points on an n-by-n grid in
    an L-shaped domain made from 3/4 of the entire square.
    adapted from C. Moler, 7-16-91, 12-22-93.
    Copyright (c) 1984-94 by The MathWorks, Inc.
    """
    x = np.ones((n,1))*np.linspace(-1,1,n)
    y = np.flipud(x.T)
    G = (x > -1) & (x < 1) & (y > -1) & (y < 1) & ( (x > 0) | (y > 0))
    G = np.where(G,1,0) # boolean to integer
    k = np.where(G)
    G[k] = 1+np.arange(len(k[0]))
    return G

def delsq(G):
    """
    DELSQ  Construct five-point finite difference Laplacian.
    delsq(G) is the sparse form of the two-dimensional,
    5-point discrete negative Laplacian on the grid G.
    adapted from  C. Moler, 7-16-91.
    Copyright (c) 1984-94 by The MathWorks, Inc.
    """
    [m,n] = G.shape
    # Indices of interior points
    G1 = G.flatten()
    p = np.where(G1)[0]
    N = len(p)
    # Connect interior points to themselves with 4's.
    i = G1[p]-1
    j = G1[p]-1
    s = 4*np.ones(p.shape)

    # for k = north, east, south, west
    for k in [-1, m, 1, -m]:
       # Possible neighbors in k-th direction
       Q = G1[p+k]
       # Index of points with interior neighbors
       q = np.where(Q)[0]
       # Connect interior points to neighbors with -1's.
       i = np.concatenate([i, G1[p[q]]-1])
       j = np.concatenate([j,Q[q]-1])
       s = np.concatenate([s,-np.ones(q.shape)])
    # sparse matrix with 5 diagonals
    return sparse.csr_matrix((s, (i,j)),(N,N))

if __name__ == '__main__':
    print numgrid(6)
    print delsq(numgrid(6)).todense()
    G = numgrid(32)
    D = delsq(G)
    N = D.shape[0]
    rhs = np.ones((N,1))
    u = linalg.spsolve(D, rhs) # vector solution
    U = np.zeros(G.shape) # map u back onto G space
    U[G>0] = u[G[G>0]-1]
    plt.contour(U)
    plt.show()

结果：

Answer 3

也许我误解了，但是你的问题似乎是关于boolean indexing and assignment，而不是与使用delsq的特定示例有关，等等

为了扩展我的评论，听起来你只是想：（使用x而不是u来避免混淆大写）

U = G.copy()
U[G > 0] = x

作为一个更具体的例子，假设您有一个数组（让我们使用这个例子中的一个，ipython的%paste x命令使这个过程更简单）：

^{pr2}$

这就产生了：

[[  0   0   0   0   0   0   0   0   0   0   0   0]
 [  0 385 417 399 394 104 474 328 396 129 230   0]
 [  0 189 440 298 120 278 115 148 142 454 405   0]
 [  0 490 484 240 312 339 192 212 287 468 225   0]
 [  0 348 300 486 485 162 159 258 418 335 110   0]
 [  0 495 496 392 364 215 122 213 222 412 122   0]
 [  0   0   0   0   0   0 442 188 325 248 225   0]
 [  0   0   0   0   0   0 431 233 141 307 339   0]
 [  0   0   0   0   0   0 259 325 102 131 333   0]
 [  0   0   0   0   0   0 458 381 127 333 170   0]
 [  0   0   0   0   0   0 454 229 162 216 192   0]
 [  0   0   0   0   0   0   0   0   0   0   0   0]]

另一方面，matlab示例处理的是稀疏数组。scipy.sparse提供了多种稀疏数组类型，但并非所有类型都支持这种类型的索引。但是，如果处理稀疏数组，只需sparse_array.data[:] = newvals来更改值，同时保持稀疏结构。在