Is Seidel's algorithm on wikipedia page incorrect?

Question

I am trying to use the python program from https://en.wikipedia.org/wiki/Seidel%27s_algorithm .

from numpy import matrix
def apd(A, n: int):
    """Compute the shortest-paths lengths."""
    if all(print(i,j)or A[i][j] for i in range(n) for j in range(n) if i != j):
        return A
    Z=A**2
    B=matrix([[1 if i != j and (A[i][j] == 1 or Z[i][j] > 0) else 0 for j in range(n)]for i in range(n)])
    T=apd(B,n)
    X = T*A
    degree = [sum(A[i][j] for j in range(n)) for i in range(n)]
    D = matrix([[2 * T[i][j] if X[i][j] >= T[i][j] * degree[j] else 2 * T[i][j] - 1 for j in range(n)]for i in range(n)])
    return D
maatriks=matrix([
         [0,0,1,0,0],
         [0,0,1,0,1],
         [1,1,0,1,1],
         [0,0,1,0,0],
         [0,1,1,0,0]])
print(apd(maatriks,5))

But I gives error:

Traceback (most recent call last):
  File "[PATH]/Seidel.py", line 19, in <module>
    print(apd(maatriks,5))
  File "[PATH]/Seidel.py", line 4, in apd
    if all(print(i,j)or A[i][j] for i in range(n) for j in range(n) if i != j):
  File "[PATH]/Seidel.py", line 4, in <genexpr>
    if all(print(i,j)or A[i][j] for i in range(n) for j in range(n) if i != j):
  File "/usr/local/lib/python3.8/dist-packages/numpy/matrixlib/defmatrix.py", line 193, in __getitem__
    out = N.ndarray.__getitem__(self, index)
IndexError: index 1 is out of bounds for axis 0 with size 1

Am I using the program wrongly or is the program itself incorrect? Here on page 9 seems to be the exactly same algorithm as is implemented in python.

Answer 1

This implementation worked for me. The code provided on wikipedia has the indexing incorrect. when indexing a numpy matrix to get the i'th , j'th element you need to do A[i,j] not A[i][j]

from numpy import matrix


def apd(A, n: int):
    """Compute the shortest-paths lengths."""
    if all(A[i, j] for i in range(n) for j in range(n) if i != j):
        return A
    Z = A ** 2
    B = matrix(
        [
            [1 if i != j and (A[i, j] == 1 or Z[i, j] > 0) else 0 for j in range(n)]
            for i in range(n)
        ]
    )
    T = apd(B, n)
    X = T * A
    degree = [sum(A[i, j] for j in range(n)) for i in range(n)]
    D = matrix(
        [
            [
                2 * T[i, j] if X[i, j] >= T[i, j] * degree[j] else 2 * T[i, j] - 1
                for j in range(n)
            ]
            for i in range(n)
        ]
    )
    return D

a = matrix(
    [
        [0, 0, 1, 0, 0],
        [0, 0, 1, 0, 1],
        [1, 1, 0, 1, 1],
        [0, 0, 1, 0, 0],
        [0, 1, 1, 0, 0],
    ]
)
print(apd(a, 5))

Returns this.

[[0 2 1 2 2]
 [2 0 1 2 1]
 [1 1 0 1 1]
 [2 2 1 0 2]
 [2 1 1 2 0]]

Is this the response you expect to see?