code: def _estimate_additive_constant(matrix):
      
      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)))
docs: 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)