code: def test_gradient():
      from brainiak.reprsimil.brsa import GBRSA
      import brainiak.utils.utils as utils
      import numpy as np
      import os.path
      import numdifftools as nd
  
      np.random.seed(100)
      file_path = os.path.join(
          os.path.dirname(__file__), "example_design.1D"
      )
      # Load an example design matrix
      design = utils.ReadDesign(fname=file_path)
  
      # concatenate it by 1, 2, and 3 times, mimicking different length
      # of experiments for different participants
      n_run = [1, 2, 1]
      design_mat = [None] * 3
      n_T = [None] * 3
      n_V = [30, 30, 20]
      for i in range(3):
          design_mat[i] = np.tile(
              design.design_task[:, :-1], [n_run[i], 1]
          )
          n_T[i] = n_run[i] * design.n_TR
  
      # start simulating some data
      n_C = np.size(design_mat[0], axis=1)
  
      noise_bot = 0.5
      noise_top = 1.5
      noise_level = [None] * 3
      for i in range(3):
          noise_level[i] = (
              np.random.rand(n_V[i]) * (noise_top - noise_bot)
              + noise_bot
          )
      # noise level is random.
  
      # AR(1) coefficient
      rho1_top = 0.8
      rho1_bot = -0.2
      rho1 = [None] * 3
  
      # generating noise
      noise = [None] * 3
  
      # baseline
      inten = [None] * 3
      for i in range(3):
          rho1[i] = (
              np.random.rand(n_V[i]) * (rho1_top - rho1_bot) + rho1_bot
          )
          noise[i] = np.zeros([n_T[i], n_V[i]])
          noise[i][0, :] = (
              np.random.randn(n_V[i])
              * noise_level[i]
              / np.sqrt(1 - rho1[i] ** 2)
          )
          for i_t in range(1, n_T[i]):
              noise[i][i_t, :] = (
                  noise[i][i_t - 1, :] * rho1[i]
                  + np.random.randn(n_V[i]) * noise_level[i]
              )
          noise[i] = noise[i] + np.dot(
              np.random.randn(n_T[i], 2), np.random.randn(2, n_V[i])
          )
          inten[i] = np.random.rand(n_V[i]) * 20.0
  
      # ideal covariance matrix
      ideal_cov = np.zeros([n_C, n_C])
      ideal_cov = np.eye(n_C) * 0.6
      ideal_cov[0:4, 0:4] = 0.2
      for cond in range(0, 4):
          ideal_cov[cond, cond] = 2
      ideal_cov[5:9, 5:9] = 0.9
      for cond in range(5, 9):
          ideal_cov[cond, cond] = 1
      L_full = np.linalg.cholesky(ideal_cov)
  
      # generating signal
      snr_top = 5.0  # test with high SNR
      snr_bot = 1.0
      # snr = np.random.rand(n_V)*(snr_top-snr_bot)+snr_bot
      # Notice that accurately speaking this is not snr. the magnitude of signal
      # depends not only on beta but also on x.
  
      snr = [None] * 3
      signal = [None] * 3
      betas_simulated = [None] * 3
      scan_onsets = [None] * 3
      Y = [None] * 3
      for i in range(3):
          snr[i] = (
              np.random.rand(n_V[i]) * (snr_top - snr_bot) + snr_bot
          )
          sqrt_v = noise_level[i] * snr[i]
          betas_simulated[i] = (
              np.dot(L_full, np.random.randn(n_C, n_V[i])) * sqrt_v
          )
          signal[i] = np.dot(design_mat[i], betas_simulated[i])
  
          # Adding noise to signal as data
          Y[i] = signal[i] + noise[i] + inten[i]
  
          scan_onsets[i] = np.linspace(0, n_T[i], num=n_run[i] + 1)
  
      # Get some initial fitting.
      SNR_bins = 11
      rho_bins = 20
      gbrsa = GBRSA(
          n_iter=3,
          rank=n_C,
          SNR_bins=SNR_bins,
          rho_bins=rho_bins,
          logS_range=0.5,
      )
  
      n_grid = SNR_bins * rho_bins
      half_log_det_X0TAX0 = [np.random.randn(n_grid) for i in range(3)]
      log_weights = np.random.randn(n_grid)
      log_fixed_terms = [np.random.randn(n_grid) for i in range(3)]
      l_idx = np.tril_indices(n_C)
      L_vec = np.random.randn(int(n_C * (n_C + 1) / 2))
      n_X0 = [2, 2, 2]
      s = np.linspace(1, SNR_bins, n_grid)
      a = np.linspace(0.5, 1, n_grid)
      s2XTAcorrX = [None] * 3
      YTAcorrY_diag = [None] * 3
      sXTAcorrY = [None] * 3
      # The calculations below are quite arbitrary and do not conform
      # to the model. They simply conform to the symmetry property and shape of
      # the matrix indicated by the model
      for i in range(3):
          YTAcorrY_diag[i] = np.sum(Y[i] * Y[i], axis=0) * a[:, None]
          s2XTAcorrX[i] = (
              np.dot(design_mat[i].T, design_mat[i])
              * s[:, None, None] ** 2
              * a[:, None, None]
          )
          sXTAcorrY[i] = (
              np.dot(design_mat[i].T, Y[i])
              * s[:, None, None]
              * a[:, None, None]
          )
  
      # test if the gradients are correct
      # print(log_fixed_terms)
      ll0, deriv0 = gbrsa._sum_loglike_marginalized(
          L_vec,
          s2XTAcorrX,
          YTAcorrY_diag,
          sXTAcorrY,
          half_log_det_X0TAX0,
          log_weights,
          log_fixed_terms,
          l_idx,
          n_C,
          n_T,
          n_V,
          n_X0,
          n_grid,
          rank=None,
      )
      # We test the gradient to the Cholesky factor
      vec = np.random.randn(np.size(L_vec))
      vec = vec / np.linalg.norm(vec)
      dd = nd.directionaldiff(
          lambda x: gbrsa._sum_loglike_marginalized(
              x,
              s2XTAcorrX,
              YTAcorrY_diag,
              sXTAcorrY,
              half_log_det_X0TAX0,
              log_weights,
              log_fixed_terms,
              l_idx,
              n_C,
              n_T,
              n_V,
              n_X0,
              n_grid,
              rank=None,
          )[0],
          L_vec,
          vec,
      )
      assert np.isclose(
          dd, np.dot(deriv0, vec), rtol=1e-5
      ), "gradient incorrect"
  
label_desc: Use image node_id 3 for calling a global function with example usage: test_gradient() without return types