NumPy 相当于 Matlab 的 magic()
NumPy equivalent of Matlab's magic()
在 Ocatave / Matlab 中,我可以使用 magic()
得到幻方,例如,
magic(4)
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
定义:幻方是一个N×N的数字网格,其中每行、每列和主对角线的条目总和为相同的数字(等于N(N^2+1)/2
)。
如何使用 NumPy 生成相同的内容?
这里有一个针对奇数和双偶数情况的快速实现。
def magic_odd(n):
if n % 2 == 0:
raise ValueError('n must be odd')
return np.mod((np.arange(n)[:, None] + np.arange(n)) + (n-1)//2+1, n)*n + \
np.mod((np.arange(1, n+1)[:, None] + 2*np.arange(n)), n) + 1
def magic_double_even(n):
if n % 4 != 0:
raise ValueError('n must be a multiple of 4')
M = np.empty([n, n], dtype=int)
M[:, :n//2] = np.arange(1, n**2//2+1).reshape(-1, n).T
M[:, n//2:] = np.flipud(M[:, :n//2]) + (n**2//2)
M[1:n//2:2, :] = np.fliplr(M[1:n//2:2, :])
M[n//2::2, :] = np.fliplr(M[n//2::2, :])
return M
奇数来自 here and I got the rest from How to construct magic squares of even order。然后我偷懒了单个偶数的情况,但想法是相似的。
这个实现遵循 Matlab 的,应该给出完全相同的结果,但有以下例外:如果 n < 3,它会抛出一个错误,而不是 return 当 n=2 时,一个非幻方 [[1, 3], [4, 2]]
就像 Matlab 一样。
照例分三种情况:奇数能被4整除,偶数不能被4整除,最后一种最复杂。
def magic(n):
n = int(n)
if n < 3:
raise ValueError("Size must be at least 3")
if n % 2 == 1:
p = np.arange(1, n+1)
return n*np.mod(p[:, None] + p - (n+3)//2, n) + np.mod(p[:, None] + 2*p-2, n) + 1
elif n % 4 == 0:
J = np.mod(np.arange(1, n+1), 4) // 2
K = J[:, None] == J
M = np.arange(1, n*n+1, n)[:, None] + np.arange(n)
M[K] = n*n + 1 - M[K]
else:
p = n//2
M = magic(p)
M = np.block([[M, M+2*p*p], [M+3*p*p, M+p*p]])
i = np.arange(p)
k = (n-2)//4
j = np.concatenate((np.arange(k), np.arange(n-k+1, n)))
M[np.ix_(np.concatenate((i, i+p)), j)] = M[np.ix_(np.concatenate((i+p, i)), j)]
M[np.ix_([k, k+p], [0, k])] = M[np.ix_([k+p, k], [0, k])]
return M
我也写了一个函数来测试这个:
def test_magic(ms):
n = ms.shape[0]
s = n*(n**2+1)//2
columns = np.all(ms.sum(axis=0) == s)
rows = np.all(ms.sum(axis=1) == s)
diag1 = np.diag(ms).sum() == s
diag2 = np.diag(ms[::-1, :]).sum() == s
return columns and rows and diag1 and diag2
尝试 [test_magic(magic(n)) for n in range(3, 20)]
检查正确性。
我遇到了同样的问题,这是我使用的:
import numpy as np
matrix = np.random.random((15,15))
for x in range(15):
for y in range(15):
matrix[x][y] = int(matrix[x][y]*10)
我需要 0 到 10 之间的整数,但你明白了...
在 Ocatave / Matlab 中,我可以使用 magic()
得到幻方,例如,
magic(4)
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
定义:幻方是一个N×N的数字网格,其中每行、每列和主对角线的条目总和为相同的数字(等于N(N^2+1)/2
)。
如何使用 NumPy 生成相同的内容?
这里有一个针对奇数和双偶数情况的快速实现。
def magic_odd(n):
if n % 2 == 0:
raise ValueError('n must be odd')
return np.mod((np.arange(n)[:, None] + np.arange(n)) + (n-1)//2+1, n)*n + \
np.mod((np.arange(1, n+1)[:, None] + 2*np.arange(n)), n) + 1
def magic_double_even(n):
if n % 4 != 0:
raise ValueError('n must be a multiple of 4')
M = np.empty([n, n], dtype=int)
M[:, :n//2] = np.arange(1, n**2//2+1).reshape(-1, n).T
M[:, n//2:] = np.flipud(M[:, :n//2]) + (n**2//2)
M[1:n//2:2, :] = np.fliplr(M[1:n//2:2, :])
M[n//2::2, :] = np.fliplr(M[n//2::2, :])
return M
奇数来自 here and I got the rest from How to construct magic squares of even order。然后我偷懒了单个偶数的情况,但想法是相似的。
这个实现遵循 Matlab 的,应该给出完全相同的结果,但有以下例外:如果 n < 3,它会抛出一个错误,而不是 return 当 n=2 时,一个非幻方 [[1, 3], [4, 2]]
就像 Matlab 一样。
照例分三种情况:奇数能被4整除,偶数不能被4整除,最后一种最复杂。
def magic(n):
n = int(n)
if n < 3:
raise ValueError("Size must be at least 3")
if n % 2 == 1:
p = np.arange(1, n+1)
return n*np.mod(p[:, None] + p - (n+3)//2, n) + np.mod(p[:, None] + 2*p-2, n) + 1
elif n % 4 == 0:
J = np.mod(np.arange(1, n+1), 4) // 2
K = J[:, None] == J
M = np.arange(1, n*n+1, n)[:, None] + np.arange(n)
M[K] = n*n + 1 - M[K]
else:
p = n//2
M = magic(p)
M = np.block([[M, M+2*p*p], [M+3*p*p, M+p*p]])
i = np.arange(p)
k = (n-2)//4
j = np.concatenate((np.arange(k), np.arange(n-k+1, n)))
M[np.ix_(np.concatenate((i, i+p)), j)] = M[np.ix_(np.concatenate((i+p, i)), j)]
M[np.ix_([k, k+p], [0, k])] = M[np.ix_([k+p, k], [0, k])]
return M
我也写了一个函数来测试这个:
def test_magic(ms):
n = ms.shape[0]
s = n*(n**2+1)//2
columns = np.all(ms.sum(axis=0) == s)
rows = np.all(ms.sum(axis=1) == s)
diag1 = np.diag(ms).sum() == s
diag2 = np.diag(ms[::-1, :]).sum() == s
return columns and rows and diag1 and diag2
尝试 [test_magic(magic(n)) for n in range(3, 20)]
检查正确性。
我遇到了同样的问题,这是我使用的:
import numpy as np
matrix = np.random.random((15,15))
for x in range(15):
for y in range(15):
matrix[x][y] = int(matrix[x][y]*10)
我需要 0 到 10 之间的整数,但你明白了...