MATLABの魔法のnumpyの等価（）Ocatave/MATLABでは

Question

、私は例えば、魔法の人口正方行列を取得するために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を使用して同じものを生成するにはどうすればよいですか？

Answer 1

Matlabのようにn = 2のときに非マジックの正方形[[1, 3], [4, 2]]を返すのではなく、n <の場合はエラーをスローします。

いつものように、奇数、4で割り切れる、そして4で割っても割り切れない、最後のものが最も複雑です。

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)]を試してみてください：

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 

は私もこれをテストする機能を書きました。

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 

奇数の場合は、hereからだと私はHow to construct magic squares of even orderから休息を得ました。その後、私は単一の偶数の場合には怠け者になりましたが、アイデアは似ています。