code: def _generate_noise_temporal_drift(
      trs,
      tr_duration,
      basis="cos_power_drop",
      period=150,
  ):
      """Generate the drift noise
  
      Create a trend (either sine or discrete_cos), of a given period and random
      phase, to represent the drift of the signal over time
  
      Parameters
      ----------
  
      trs : int
          How many volumes (aka TRs) are there
  
      tr_duration : float
          How long in seconds is each volume acqusition
  
      basis : str
          What is the basis function for the drift. Could be made of discrete
          cosines (for longer run durations, more basis functions are
          created) that either have equal power ('discrete_cos') or the power
          diminishes such that 99% of the power is below a specified frequency
          ('cos_power_drop'). Alternatively, this drift could simply be a sine
          wave ('sine')
  
      period : int
          When the basis function is 'cos_power_drop' this is the period over
          which no power of the drift exceeds (i.e. the power of the drift
          asymptotes at this period). However for the other basis functions,
          this is simply how many seconds is the period of oscillation of the
          drift
  
      Returns
      ----------
      noise_drift : one dimensional array, float
          The drift timecourse of activity
  
      """
  
      # Calculate drift differently depending on the basis function
      if basis == "discrete_cos":
          # Specify each tr in terms of its phase with the given period
          timepoints = np.linspace(0, trs - 1, trs)
          timepoints = ((timepoints * tr_duration) / period) * 2 * np.pi
  
          # Specify the other timing information
          duration = trs * tr_duration
          basis_funcs = int(
              np.floor(duration / period)
          )  # How bases do you have
  
          if basis_funcs == 0:
              err_msg = (
                  "Too few timepoints (" + str(trs) + ") to accurately "
                  "model drift"
              )
              logger.warning(err_msg)
              basis_funcs = 1
  
          noise_drift = np.zeros((timepoints.shape[0], basis_funcs))
          for basis_counter in list(range(1, basis_funcs + 1)):
              # What steps do you want to take for this basis function
              timepoints_basis = (timepoints / basis_counter) + (
                  np.random.rand() * np.pi * 2
              )
  
              # Store the drift from this basis func
              noise_drift[:, basis_counter - 1] = np.cos(
                  timepoints_basis
              )
  
          # Average the drift
          noise_drift = np.mean(noise_drift, 1)
  
      elif basis == "sine":
          # Calculate the cycles of the drift for a given function.
          cycles = trs * tr_duration / period
  
          # Create a sine wave with a given number of cycles and random phase
          timepoints = np.linspace(0, trs - 1, trs)
          phaseshift = np.pi * 2 * np.random.random()
          phase = (
              timepoints / (trs - 1) * cycles * 2 * np.pi
          ) + phaseshift
          noise_drift = np.sin(phase)
  
      elif basis == "cos_power_drop":
          # Make a vector counting each TR
          timepoints = np.linspace(0, trs - 1, trs) * tr_duration
  
          # Specify the other timing information
          duration = trs * tr_duration
  
          # How bases do you have? This is to adhere to Nyquist
          basis_funcs = int(trs)
  
          noise_drift = np.zeros((timepoints.shape[0], basis_funcs))
          for basis_counter in list(range(1, basis_funcs + 1)):
              # What steps do you want to take for this basis function
              random_phase = np.random.rand() * np.pi * 2
  
              timepoint_phase = (
                  timepoints / duration * np.pi * basis_counter
              )
  
              # In radians, what is the value for each time point
              timepoints_basis = timepoint_phase + random_phase
  
              # Store the drift from this basis func
              noise_drift[:, basis_counter - 1] = np.cos(
                  timepoints_basis
              )
  
          def power_drop(r, L, F, tr_duration):
              # Function to return the drop rate for the power of basis functions
              # In other words, how much should the weight of each basis function
              # reduce in order to make the power you retain of the period's
              # frequency be 99% of the total power of the highest frequency, as
              # defined by the DCT.
              # For an example where there are 20 time points, there will be 20
              # basis functions in the DCT. If the period of the signal you wish
              # to simulate is such that 99% of the power should drop off after
              # the equivalent of 5 of these basis functions, then the way this
              # code works is it finds the rate at which power must drop off for
              # all of the 20 basis functions such that by the 5th one, there is
              # only 1% of the power remaining.
              # r is the power reduction rate which should be between 0 and 1
              # L is the duration of the run in seconds
              # F is period of the cycle in seconds It is assumed that this will
              # be greater than the tr_duration, or else this will not work
              # tr_duration is the duration of each TR in seconds
  
              # Check the TR duration
              if F < tr_duration:
                  msg = "Period %0.0f > TR duration %0.0f" % (
                      (F, tr_duration)
                  )
                  raise ValueError(msg)
  
              percent_retained = (
                  0.99  # What is the percentage of power retained
              )
  
              # Compare the power at the period frequency (in the numerator) with
              # the power at the frequency of the DCT, AKA the highest possible
              # frequency in the data (in the denominator)
              numerator = 1 - r ** (2 * L / F)  # Power of this period
              denominator = 1 - r ** (
                  2 * L / tr_duration
              )  # Power of DCT freq.
  
              # Calculate the retained power
              power_drop = abs(
                  (numerator / denominator) - percent_retained
              )
              return power_drop
  
          # Solve for power reduction rate.
          # This assumes that r is between 0 and 1
          # Takes the duration and period as arguments
          sol = optimize.minimize_scalar(
              power_drop,
              bounds=(0, 1),
              method="Bounded",
              args=(duration, period, tr_duration),
          )
  
          # Pull out the solution
          r = sol.x
  
          # Weight the basis functions based on the power drop off
          basis_weights = r ** np.arange(basis_funcs)
  
          # Weigh the basis functions
          weighted_basis_funcs = np.multiply(noise_drift, basis_weights)
  
          # Average the drift
          noise_drift = np.mean(weighted_basis_funcs, 1)
  
      # Normalize so the sigma is 1
      noise_drift = stats.zscore(noise_drift)
  
      # Return noise
      return noise_drift
  
imports: import logging
  from itertools import product
  from statsmodels.tsa.arima.model import ARIMA
  import math
  import numpy
  from numpy.linalg import LinAlgError
  from pkg_resources import resource_stream
  from scipy import stats
  from scipy import signal
  import scipy.ndimage
  import copy
  from scipy import optimize
label_desc: Use image node_id 17 for calling a global function with example usage: _generate_noise_temporal_drift(trs, tr_duration, basis, period) and returns: noise_drift, power_drop