code: def test_fit():
      from brainiak.reprsimil.brsa import GBRSA
      import brainiak.utils.utils as utils
      import scipy.stats
      import numpy as np
      import os.path
  
      np.random.seed(10)
      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 = [2, 1, 1]
      design_mat = [None] * 3
      n_T = [None] * 3
      n_V = [40, 60, 60]
      for i in range(3):
          design_mat[i] = np.tile(
              design.design_task[:, :-1], [n_run[i], 1]
          )
  
      # 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
  
      # 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):
          design_mat[i] = np.tile(
              design.design_task[:, :-1], [n_run[i], 1]
          )
          n_T[i] = n_run[i] * design.n_TR
          noise_level[i] = (
              np.random.rand(n_V[i]) * (noise_top - noise_bot)
              + noise_bot
          )
          # noise level is random.
          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)
  
      # Test fitting.
      n_nureg = 2
      gbrsa = GBRSA(
          n_iter=15,
          auto_nuisance=True,
          logS_range=0.5,
          SNR_bins=11,
          rho_bins=16,
          n_nureg=n_nureg,
          optimizer="L-BFGS-B",
      )
  
      gbrsa.fit(X=Y, design=design_mat, scan_onsets=scan_onsets)
  
      # Check that result is significantly correlated with the ideal covariance
      # matrix
      u_b = gbrsa.U_
      u_i = ideal_cov
      p = scipy.stats.spearmanr(
          u_b[np.tril_indices_from(u_b)], u_i[np.tril_indices_from(u_i)]
      )[1]
      assert p < 0.01, (
          "Fitted covariance matrix does not correlate with ideal covariance "
          "matrix!"
      )
      # check that the recovered SNRs makes sense
      p = scipy.stats.pearsonr(gbrsa.nSNR_[0], snr[0])[1]
      assert (
          p < 0.01
      ), "Fitted SNR does not correlate with simulated SNR!"
      p = scipy.stats.pearsonr(gbrsa.sigma_[1], noise_level[1])[1]
      assert (
          p < 0.01
      ), "Fitted noise level does not correlate with simulated noise level!"
      p = scipy.stats.pearsonr(gbrsa.rho_[2], rho1[2])[1]
      assert (
          p < 0.01
      ), "Fitted AR(1) coefficient does not correlate with simulated values!"
      assert np.shape(gbrsa.X0_[1]) == (
          n_T[1],
          n_nureg + 1,
      ), "Wrong size of X0"
  
      Y_new = [None] * 3
      noise_new = [None] * 3
      for i in range(3):
          noise_new[i] = np.zeros([n_T[i], n_V[i]])
          noise_new[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_new[i][i_t, :] = (
                  noise_new[i][i_t - 1, :] * rho1[i]
                  + np.random.randn(n_V[i]) * noise_level[i]
              )
  
          Y_new[i] = signal[i] + noise_new[i] + inten[i]
      ts, ts0 = gbrsa.transform(Y_new, scan_onsets=scan_onsets)
      [score, score_null] = gbrsa.score(
          X=Y_new, design=design_mat, scan_onsets=scan_onsets
      )
      [score_noise, score_null_noise] = gbrsa.score(
          X=noise_new, design=design_mat, scan_onsets=scan_onsets
      )
      for i in range(3):
          assert np.shape(ts[i]) == (n_T[i], n_C) and np.shape(
              ts0[i]
          ) == (n_T[i], n_nureg + 1)
          p = scipy.stats.pearsonr(ts[i][:, 0], design_mat[i][:, 0])[1]
          assert (
              p < 0.01
          ), "Recovered time series does not correlate with true time series!"
  
          assert score[i] > score_null[i], (
              "Full model does not win over null model on data containing "
              "signal"
          )
  
          assert (
              score_noise[i] < score_null_noise[i]
          ), "Null model does not win over full model on data without signal"
  
      [score, score_null] = gbrsa.score(
          X=[None] * 3, design=design_mat, scan_onsets=scan_onsets
      )
      assert (
          score == [None] * 3 and score_null == [None] * 3
      ), "score did not return list of None when data is None"
      ts, ts0 = gbrsa.transform(X=[None] * 3, scan_onsets=scan_onsets)
      assert (
          ts == [None] * 3 and ts0 == [None] * 3
      ), "transform did not return list of None when data is None"
  
label_desc: Use image node_id 2 for calling a global function with example usage: test_fit() without return types