code: def haversine_grad(x, y):
      # spectral initialization puts many points near the poles
      # currently, adding pi/2 to the latitude avoids problems
      # TODO: reimplement with quaternions to avoid singularity
  
      if x.shape[0] != 2:
          raise ValueError(
              "haversine is only defined for 2 dimensional data"
          )
      sin_lat = np.sin(0.5 * (x[0] - y[0]))
      cos_lat = np.cos(0.5 * (x[0] - y[0]))
      sin_long = np.sin(0.5 * (x[1] - y[1]))
      cos_long = np.cos(0.5 * (x[1] - y[1]))
  
      a_0 = (
          np.cos(x[0] + np.pi / 2)
          * np.cos(y[0] + np.pi / 2)
          * sin_long**2
      )
      a_1 = a_0 + sin_lat**2
  
      d = 2.0 * np.arcsin(np.sqrt(min(max(abs(a_1), 0), 1)))
      denom = np.sqrt(abs(a_1 - 1)) * np.sqrt(abs(a_1))
      grad = np.array(
          [
              (
                  sin_lat * cos_lat
                  - np.sin(x[0] + np.pi / 2)
                  * np.cos(y[0] + np.pi / 2)
                  * sin_long**2
              ),
              (
                  np.cos(x[0] + np.pi / 2)
                  * np.cos(y[0] + np.pi / 2)
                  * sin_long
                  * cos_long
              ),
          ]
      ) / (denom + 1e-6)
      return d, grad
  
imports: import numba
  import numpy
  import scipy.stats
  from sklearn.metrics import pairwise_distances
label_desc: Use image node_id 32 for calling a global function with example usage: haversine_grad(x, y) and returns: d,  grad