CMDS Additive Constant: correction for non-Euclidean distances.
    Procedure taken from R function cmdscale.
    The additive constant is given by the largest eigenvalue (real part)
    of this 2x2 block matrix -

    /-------+-------\
    |       |       |
    |   0   | 2*dbc |    NB: dbc function(m: matrix):
    |       | (d^2) |        double_centre() [see above]
    +-------+-------+
    |       |       |
    |  -I   |-4*dbc |
    |       |  (2d) |
    \-------+-------/

    corrected matrix = matrix + additive constant (diagonal kept as 0)

def _estimate_additive_constant(matrix):
    """
    CMDS Additive Constant: correction for non-Euclidean distances.
    Procedure taken from R function cmdscale.
    The additive constant is given by the largest eigenvalue (real part)
    of this 2x2 block matrix -

    /-------+-------\
    |       |       |
    |   0   | 2*dbc |    NB: dbc function(m: matrix):
    |       | (d^2) |        double_centre() [see above]
    +-------+-------+
    |       |       |
    |  -I   |-4*dbc |
    |       |  (2d) |
    \-------+-------/

    corrected matrix = matrix + additive constant (diagonal kept as 0)
    """
    topleft = np.zeros(matrix.shape)
    topright = 2*double_centre(matrix)
    bottomleft = -np.eye(matrix.shape[0])
    bottomright = -4*double_centre(matrix, square_input=False)
    Z = np.vstack([np.hstack([topleft,topright]),
                   np.hstack([bottomleft,bottomright])])
    return max(np.real(np.linalg.eigvals(Z)))