NumPy的等效于Matlab中的magic()函数

Question

在Octave/Matlab中，我可以使用magic()函数获得一个幻方，例如：

magic(4)

  16    2    3   13
   5   11   10    8
   9    7    6   12
   4   14   15    1

定义：一个幻方（magic square）是一个N×N的数字网格，在其中每一行、每一列和主对角线上的数之和都等于相同的数字（等于N(N^2+1)/2）。

我该如何使用NumPy生成相同的幻方？

- Tom Hale

---

好的，请发送需要翻译的内容。 - m7913d

@m7913d，请将其作为答案发布，因为链接可能会失效... - MaxU - stand with Ukraine

@MaxU版权问题怎么样？ - m7913d

你看过 MATLAB 代码了吗？我怀疑它是一个 .m 文件（可读的 MATLAB）。Octave 版本是可读的（如果需要，我可以发布）。我的记忆是这是很久以前放入 MATLAB 中的，更多的是作为展示而不是什么复杂或有用的东西。 - hpaulj

所以答案是你不能这样做（除非你自己编写代码）。 - liyuan

Answer 1

这个实现遵循Matlab的实现，并且应该会给出完全相同的结果，但有一个例外：当n < 3时，它会抛出一个错误，而不是像Matlab一样返回一个非幻方[[1, 3], [4, 2]]，其中n=2。

通常情况下，有三种情况：奇数阶、可以被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)] 来检查正确性。

- user6655984

Answer 2

这里提供奇数和双偶数情况的快速实现方法。

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

奇数阶幻方的方法来自这里，其余内容来自如何构造偶阶幻方。对于单数偶阶的情况，我变得有些懒惰，但思路是类似的。

- percusse

Answer 3

我有同样的问题，这是我使用的代码：

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之间的整数，但你已经明白了...

- Mark

---

在numpy中，我们可以使用arr[x, y]对数组进行索引；但是你不需要循环，只需使用matrix = (matrix*10).astype(int)。 - hpaulj

np.random.permutation(16).reshape(4,4) 或者写成一行。 - Tom Hale

当然，这个随机矩阵不一定是一个幻方。 - user6655984