python np.dstackの使い方について

Question

先日、以下のような質問をしました
https://teratail.com/questions/152409

この質問に対して寄せていただいた回答をもとに、8*8の行列で各成分が0～2πとなるものは以下のようにあらわすことができました。

import numpy as np
S = np.random.uniform(0, 2 * np.pi, (Nx, Ny))
print(S)

↓

[[4.41650775 2.71406146 4.41416128 3.14201439 2.69117847 4.27302197
  2.13271308 3.04400582]
 [1.12690231 3.42268815 4.46682612 4.73870683 5.92194469 2.70762728
  5.3183137  3.73280308]
 [2.39379053 4.57429553 2.89689483 1.63613895 2.62612769 2.82020993
  1.662676   2.52406824]
 [2.0699564  4.10302686 3.45921419 5.79703282 3.02967725 5.82700575
  2.46686931 0.37984799]
 [3.08676912 5.76371844 3.64893616 0.86066251 2.55529371 2.63367659
  3.29609631 1.03792125]
 [4.88785811 6.02503479 3.8441244  5.05772728 1.00185371 3.58283383
  1.9105011  4.4521828 ]
 [2.28177066 5.43348847 3.7151592  0.27352182 5.39425746 5.95659263
  2.02327713 4.24771053]
 [2.62177829 5.79892895 3.47362596 5.06351588 3.27111045 3.12798965
  0.02599374 3.94512495]]

これを用いて次は周期的境界条件を課して、先日の質問の中でもあった以下の値を計算したいと考えています。
S(i,j)とすると
E =-cos(S(i,j)-S(i+1,j))-cos(S(i,j)-S(i-1,j))-cos(S(i,j)-S(i,j+1))-cos(S(i,j)-S(i,j-1)))

ただし、ここで周期的境界条件というものを課して列１，８及び行１，８で反対側に折り返すという操作を行うこととします。

以上を組み込んで以下のようなコードを考えました。

import matplotlib.pyplot as plt
import numpy as np


def Ecalc(Nx,Ny):
    for i in range (-1,Nx):
        for j in range(0,Ny):
            l=i
            m=j
            if l==-1:
                l=Nx
            if l==Nx:
                l=0
            if m==-1:
                m=Ny
            if m==Ny:
                m=0

    for i in range(0,Nx):
        for j in range(-1,Ny):
            l=i
            m=j
            if l==-1:
                l=Nx
            if l==Nx:
                l=0
            if m==-1:
                m=Ny
            if m==Ny:
                m=0

    S = np.random.uniform(0, 2 * np.pi, (Nx, Ny))
    print(S)

    
    E = np.empty_like(S)
    for i, j in np.dstack(Nx,Ny):
        E[i - 1, j - 1] = -np.cos(S[i, j] - S[i - 1, j]) \
                      - np.cos(S[i, j] - S[i + 1, j]) \
                      - np.cos(S[i, j] - S[i, j + 1]) \
                      - np.cos(S[i, j] - S[i, j - 1])
    print(E)
    Esum=np.sum(E)
    print(Esum)

if __name__ == '__main__':
    Ecalc(8,8)

↓

Traceback (most recent call last):
  File "C:\Users\user2\Desktop\testtesttesttest2.py", line 46, in <module>
    Ecalc(8,8)
  File "C:\Users\user2\Desktop\testtesttesttest2.py", line 37, in Ecalc
    for i, j in np.dstack(Nx,Ny):
TypeError: dstack() takes 1 positional argument but 2 were given

このエラーを解消し、周期的境界条件を満たした値Eを出力することが目標です。よろしくお願い致します。

追記

やりたいことについて詳しく書きます。
1.NxNyの行列を生成し、各成分を(i,j)であらわす。
2.行列の各要素は０～2np.piの乱数を取るものとする。この値をS[i,j]とする。
3.E[i,j]を以下のように定義し、すべてのS[i,j]に対して計算を行う。
E[i,j] =-cos(S(i,j)-S(i+1,j))-cos(S(i,j)-S(i-1,j))-cos(S(i,j)-S(i,j+1))-cos(S(i,j)-S(i,j-1)))
しかし、このままでは行列の端では計算ができなくなる（例えばE[1.1]を計算するとき、S[0,1]やS[1,0]が定義されていないので計算ができない）ので、以下のような周期的境界条件を課します。
～～周期的境界条件～～
例えばNx=8,Ny=8の場合、S[0,1]=S[8,1]、S[1,0]=S[1,8]といった具合です。

Comments

[hayataka2049]

なぜdstackを使おうと思ったのでしょうか？　コードを読む限りあまり適切とも思えませんし、どういう動作をさせたいのかという意図がわからないと回答がつきづらいと思います

[astromelt0416]

前回の質問の中で、回答してくださった方の中にdstackを使っていらっしゃる方がいて、それに従ってやってみたらうまくいったので使っていた次第です。実際、dstack自体は3次元配列の結合の際に用いられるものとのことで、正直私も理解せず使ってしまっていました。私がやりたいことは質問文中にあります、Eを計算することです。その過程で行列要素のナンバリングと各要素i:1~Nx,j:1~Nyに対する繰り返し処理が必要となりました。質問文にも追記しておきます。

Answer 1

エラーメッセージで「位置引数が1つのはずが、2つ渡された」となっており、実際、リファレンス numpy.dstack() を確認すると、

np.dstack(Nx,Ny)

ではなく、

np.dstack([Nx,Ny])

とするよう書いてあります。

追記

前回の回答で dstack() を使用したのは、「各要素をナンバリングした行列を作りたい」と書いてあったので、np.mgrid で各要素の x 座標と y 座標の一覧を作成し、それを結合するために使いました。

indices = np.dstack(np.mgrid[1:9, 1:9]).reshape(-1, 2)
print(indices)

[[1 1]
 [1 2]
 [1 3]
 [1 4]
 [1 5]
 [1 6]
 [1 7]
...

ただ、各要素を iteratate するのが目的だったのであれば、np.ndenumerate() を使ったほうがシンプルです。

a = np.random.randn(5, 5)
print(a)

for index, value in np.ndenumerate(a):
    print(index, value)

[[1 1]
 [1 2]
 [1 3]
 [1 4]
 [1 5]
 [1 6]
 [1 7]
 [1 8]
 [2 1]
...

追記

境界条件って端の場合は反対側の値を使うということですよね。
であれば、padding='wrap' で np.pad() しておけばいいです。

どういうことかは以下を見ていただければわかるかと思います。

a = np.random.randint(0, 5, (5, 5))
print(a)
# [[3 3 2 2 0]
#  [0 4 1 2 1]
#  [1 4 4 1 4]
#  [4 4 2 4 0]
#  [4 1 4 3 3]]

b = np.pad(a, (1, 1), 'wrap')
print(b)
# [[3 4 1 4 3 3 4]
#  [0 3 3 2 2 0 3]
#  [1 0 4 1 2 1 0]
#  [4 1 4 4 1 4 1]
#  [0 4 4 2 4 0 4]
#  [3 4 1 4 3 3 4]
#  [0 3 3 2 2 0 3]]

np.set_printoptions(linewidth=120)

Nx, Ny = 5, 5
S = np.random.uniform(0, 2 * np.pi, (5, 5))
print(S)

E = np.empty_like(S)
S_pad = np.pad(S, (1, 1), 'wrap')
print(S_pad)

for i, j in np.dstack(np.mgrid[1:Nx+1, 1:Ny+1]).reshape(-1, 2):
    E[i - 1, j - 1] = -np.cos(S_pad[i, j] - S_pad[i - 1, j]) \
                      - np.cos(S_pad[i, j] - S_pad[i + 1, j]) \
                      - np.cos(S_pad[i, j] - S_pad[i, j + 1]) \
                      - np.cos(S_pad[i, j] - S_pad[i, j - 1])
print(E)

追記

話を簡単にするため (3, 3) 行列で考えます。

Nx, Ny = 3, 3
S = np.random.randint(0, 5, (Nx, Ny))
print(S)
# [[0 2 4]
#  [0 0 3]
#  [2 0 2]]

S_pad は (Nx, Ny) の配列の周囲に1ずつパディングを入れたものです。
パディングの値は padding='wrap' を指定することで、例えば、S[0, 0] であれば、S[2, 2] の値になるようにしています。
これにより、E[0, 0] = -cos(S[0, 0] - S[-1, 0]) を計算する場合、条件分岐が不要になります。

S_pad = np.pad(a, (1, 1), 'wrap')
print(S_pad)
# [[2 2 0 2 2]
#  [4 0 2 4 0]
#  [3 0 0 3 0]
#  [2 2 0 2 2]
#  [4 0 2 4 0]]

S_pad の各要素のインデックスは以下のようになっています。

(0, 4) (1, 4) (2, 4) (3, 4) (4, 4) (5, 4) (6, 4) 
(0, 4) (1, 4) (2, 4) (3, 4) (4, 4) (5, 4) (6, 4) 
(0, 4) (1, 4) (2, 4) (3, 4) (4, 4) (5, 4) (6, 4) 
(0, 4) (1, 4) (2, 4) (3, 4) (4, 4) (5, 4) (6, 4) 
(0, 4) (1, 4) (2, 4) (3, 4) (4, 4) (5, 4) (6, 4) 
(0, 4) (1, 4) (2, 4) (3, 4) (4, 4) (5, 4) (6, 4) 
(0, 4) (1, 4) (2, 4) (3, 4) (4, 4) (5, 4) (6, 4) 

このうち、元の行列 S の値がある以下のインデックスだけ iterate したいわけです。

(1, 4) (2, 4) (3, 4) (4, 4) (5, 4) 
(1, 4) (2, 4) (3, 4) (4, 4) (5, 4) 
(1, 4) (2, 4) (3, 4) (4, 4) (5, 4) 
(1, 4) (2, 4) (3, 4) (4, 4) (5, 4) 
(1, 4) (2, 4) (3, 4) (4, 4) (5, 4) 

なので、目標は以下の配列を作ることです。

[(1, 1), (1, 2), (1, 3), (1, 4), ...]

np.mgrid は返り値として次の2つの値を返します。

X, Y = np.mgrid[1:Nx + 1, 1:Ny + 1]
print(X)
# 各要素の行
# [[1 1 1]
#  [2 2 2]
#  [3 3 3]]

# 各要素の列
print(Y)
# [[1 2 3]
#  [1 2 3]
#  [1 2 3]]

よって、これを np.dstack() で結合すると、各要素のインデックスの一覧が作れます。

XY = np.dstack([X, Y])
print(XY)
# [[[1 1] [1 2] [1 3]]
#  [[2 1] [2 2] [2 3]]
#  [[3 1] [3 2] [3 3]]]

これは (3, 3, 2) の numpy 配列なので、(9, 2) として、reshape() します。
(3, 3, 2) を (9, 2) に reshape() する際に、(?, 2) と2が決まっていれば、? の部分は9と確定するので、この場合、-1 と指定しておくことができます。

print(XY.reshape(-1, 2))  # XY.reshape(9, 2) と一緒
# [[1 1]
#  [1 2]
#  [1 3]
#  [2 1]
#  [2 2]
#  [2 3]
#  [3 1]
#  [3 2]
#  [3 3]]

Comments

[astromelt0416]

回答ありがとうございます。ご指摘の通り、np.dstack([Nx,Ny])としたのですが、次は以下のようなValueErrorが出てしまっています。
ValueError: not enough values to unpack (expected 2, got 1)
この場合、np.dstackを使っていること自体が間違いなのでしょうか・・・

[tiitoi]

前回の回答では、np.dstack(np.mgrid[1:9, 1:9]).reshape(-1, 2) となっていたはずです。
1つずつ要素を走査していくのが目的であれば np.ndenumerate() を使ったほうがシンプルかと思います。

[tiitoi]

境界条件の判定は padding='wrap' で np.pad() しておけば不要です。サンプルコードを追記しました。

[astromelt0416]

意図を汲んでいただきありがとうございます。追記の方、参考にさせていただきました。想像通りの結果が得られたので、一度先に進んでみたいと思います。ありがとうございました。

[tiitoi]

想定どおり動いたようでよかったです。
ただ、コードでわからない部分があれば、よくわからないけど動くという状態で使うのではなく、理解しておいたほうがよいですよ。
コードで不明な点があれば補足します。

[astromelt0416]

おっしゃる通りです。お言葉に甘えて少し追質問させていただきます。
for i, j in np.dstack(np.mgrid[1:Nx, 1:Ny]).reshape(-1, 2):　の部分で.reshape(-1, 2)を除くと、ValueError: too many values to unpack (expected 2)が出るのを確認しました。先ほどもこのエラーで躓いていた記憶なのですが、この.reshape(-1, 2)はいったいどのような意味なのでしょうか。

[tiitoi]

先程のコードに一部誤りがあったので、訂正しました。
× np.mgrid[1:Nx, 1:Ny]
○ np.mgrid[1:Nx+1, 1:Ny+1]

np.dstack(np.mgrid[1:Nx, 1:Ny]).reshape(-1, 2) はループしたい添字の一覧を作っています。回答に追記しました。

meshgrid は最初とっつきにくいと思いますので、以下のページなども合わせて参照されるといいかと思います。

https://deepage.net/features/numpy-meshgrid.html

[astromelt0416]

回答ありがとうございました。meshgridの雰囲気はつかめたと思います。参考ページと合わせて、いろいろと試してみたいと思います。