code: def test_joint_pdf_gradients_sparse(rng):
      h = 1e-4
  
      # Make sure dictionary entries are processed in the same order regardless
      # of the platform. Otherwise any random numbers drawn within the loop
      # would make the test non-deterministic even if we fix the seed before
      # the loop.Right now, this test does not draw any samples, but we still
      # sort the entries to prevent future related failures.
  
      for ttype in sorted(factors):
          dim = ttype[1]
          if dim == 2:
              nslices = 1
              interp_method = interpolate_scalar_2d
          else:
              nslices = 45
              interp_method = interpolate_scalar_3d
  
          transform = regtransforms[ttype]
          factor = factors[ttype]
          theta = transform.get_identity_parameters()
  
          (
              static,
              moving,
              static_g2w,
              moving_g2w,
              smask,
              mmask,
              M,
          ) = setup_random_transform(
              transform, factor, nslices, 5.0, rng=rng
          )
          parzen_hist = ParzenJointHistogram(32)
          parzen_hist.setup(static, moving, smask, mmask)
  
          # Sample the fixed-image domain
          k = 3
          sigma = 0.25
          shape = np.array(static.shape, dtype=np.int32)
          samples = sample_domain_regular(
              k, shape, static_g2w, sigma, rng
          )
          samples = np.array(samples)
          samples = np.hstack(
              (samples, np.ones(samples.shape[0])[:, None])
          )
          sp_to_static = np.linalg.inv(static_g2w)
          samples_static_grid = sp_to_static.dot(samples.T).T[..., :dim]
          intensities_static, inside = interp_method(
              static.astype(np.float32), samples_static_grid
          )
          # The routines in vector_fields operate, mostly, with float32 because
          # they were thought to be used for non-linear registration. We may need
          # to write some float64 counterparts for affine registration, where
          # memory is not so big issue
          intensities_static = np.array(
              intensities_static, dtype=np.float64
          )
  
          # Compute the gradient at theta with the implementation under test
          M = transform.param_to_matrix(theta)
          sp_to_moving = np.linalg.inv(moving_g2w).dot(M)
          samples_moving_grid = sp_to_moving.dot(samples.T).T[..., :dim]
          intensities_moving, inside = interp_method(
              moving.astype(np.float32), samples_moving_grid
          )
          intensities_moving = np.array(
              intensities_moving, dtype=np.float64
          )
          parzen_hist.update_pdfs_sparse(
              intensities_static, intensities_moving
          )
          # Get the joint distribution evaluated at theta
          J0 = np.copy(parzen_hist.joint)
  
          spacing = np.ones(dim + 1, dtype=np.float64)
          mgrad, inside = vf.sparse_gradient(
              moving.astype(np.float32), sp_to_moving, spacing, samples
          )
          parzen_hist.update_gradient_sparse(
              theta,
              transform,
              intensities_static,
              intensities_moving,
              samples[..., :dim],
              mgrad,
          )
          # Get the gradient of the joint distribution w.r.t. the transform
          # parameters
          actual = np.copy(parzen_hist.joint_grad)
  
          # Compute the gradient using finite differences
          n = transform.get_number_of_parameters()
          expected = np.empty_like(actual)
          for i in range(n):
              dtheta = theta.copy()
              dtheta[i] += h
              # Update the joint distribution with the transformed moving image
              M = transform.param_to_matrix(dtheta)
              sp_to_moving = np.linalg.inv(moving_g2w).dot(M)
              samples_moving_grid = sp_to_moving.dot(samples.T).T
              intensities_moving, inside = interp_method(
                  moving.astype(np.float32), samples_moving_grid
              )
              intensities_moving = np.array(
                  intensities_moving, dtype=np.float64
              )
              parzen_hist.update_pdfs_sparse(
                  intensities_static, intensities_moving
              )
              J1 = np.copy(parzen_hist.joint)
              expected[..., i] = (J1 - J0) / h
  
          # Dot product and norms of gradients of each joint histogram cell
          # i.e. the derivatives of each cell w.r.t. all parameters
          P = (expected * actual).sum(2)
          enorms = np.sqrt((expected**2).sum(2))
          anorms = np.sqrt((actual**2).sum(2))
          prodnorms = enorms * anorms
          # Cosine of angle between the expected and actual gradients.
          # Exclude very small gradients
          P[prodnorms > 1e-6] /= prodnorms[prodnorms > 1e-6]
          P[prodnorms <= 1e-6] = 0
          # Verify that a large proportion of the gradients point almost in
          # the same direction. Disregard very small gradients
          mean_cosine = P[P != 0].mean()
          std_cosine = P[P != 0].std()
          assert mean_cosine > 0.98
          assert std_cosine < 0.16
  
imports: import numpy
  import scipy
  from functools import reduce
  from operator import mul
  from dipy.core.ndindex import ndindex
  from dipy.core.interpolation import interpolate_scalar_2d, interpolate_scalar_3d
  from dipy.data import get_fnames
  from dipy.align import vector_fields
  from dipy.align.transforms import regtransforms
  from dipy.align.parzenhist import ParzenJointHistogram, cubic_spline, cubic_spline_derivative, sample_domain_regular
  from numpy.testing import assert_array_equal, assert_array_almost_equal, assert_almost_equal, assert_equal, assert_raises
  from dipy.testing.decorators import set_random_number_generator
label_desc: Use image node_id 8 for calling a global function with example usage: test_joint_pdf_gradients_sparse(rng) without return types