Double-centres the input matrix: From each element: Subtract the row
    mean Subtract the column mean Add the grand mean Divide by -2 Method
    from: Torgerson, W S (1952). Multidimensional scaling: I. Theory and
    method. Alternatively M = -0.5 * (I - 1/n)D[^2](I - 1/n)

def double_centre(matrix, square_input=True):
    """ Double-centres the input matrix: From each element: Subtract the row
    mean Subtract the column mean Add the grand mean Divide by -2 Method
    from: Torgerson, W S (1952). Multidimensional scaling: I. Theory and
    method. Alternatively M = -0.5 * (I - 1/n)D[^2](I - 1/n) """

    m = matrix.copy()
    if square_input:
        m **= 2

    (rows, cols) = m.shape
    cm = np.mean(m, axis=0)  # column means
    rm = np.mean(m, axis=1).reshape((rows, 1))  # row means
    gm = np.mean(cm)  # grand mean
    m -= rm + cm - gm
    m /= -2
    return m