code: def _D_constraint(self, neg_pairs, w):
      
      diff = neg_pairs[:, 0, :] - neg_pairs[:, 1, :]
      diff_sq = diff * diff
      dist = np.sqrt(diff_sq.dot(w))
      sum_deri1 = np.einsum('ij,i', diff_sq, 0.5 / np.maximum(dist, 1e-6))
      sum_deri2 = np.einsum(
          'ij,ik->jk',
          diff_sq,
          diff_sq / (-4 * np.maximum(1e-6, dist**3))[:,None]
      )
      sum_dist = dist.sum()
      return (
        np.log(sum_dist),
        sum_deri1 / sum_dist,
        sum_deri2 / sum_dist - np.outer(sum_deri1, sum_deri1) / (sum_dist * sum_dist)
      )
docs: Compute the value, 1st derivative, second derivative (Hessian) of
      a dissimilarity constraint function gF(sum_ij distance(d_ij A d_ij))
      where A is a diagonal matrix (in the form of a column vector 'w').