code: def compute_signal_change(
      signal_function,
      noise_function,
      noise_dict,
      magnitude,
      method="PSC",
  ):
      """Rescale the signal to be a given magnitude, based on a specified
      metric (e.g. percent signal change). Metrics are heavily inspired by
      Welvaert & Rosseel (2013). The rescaling is based on the maximal
      activity in the timecourse. Importantly, all values within the
      signal_function are scaled to have a min of -1 or max of 1, meaning that
      the voxel value will be the same as the magnitude.
  
      Parameters
      ----------
  
  
      signal_function : timepoint by voxel array
          The signal time course to be altered. This can have
          multiple time courses specified as different columns in this
          array. Conceivably you could use the output of
          generate_stimfunction as the input but the temporal variance
          will be incorrect. Critically, different values across voxels are
          considered relative to each other, not independently. E.g., if the
          voxel has a peak signal twice as high as another voxel's, then this
          means that the signal after these transformations will still be
          twice as high (according to the metric) in the first voxel relative
          to the second
  
      noise_function : timepoint by voxel numpy array
          The time course of noise (a voxel created from generate_noise)
          for each voxel specified in signal_function. This is necessary
          for computing the mean evoked activity and the noise variability
  
      noise_dict : dict
          A dictionary specifying the types of noise in this experiment. The
          noise types interact in important ways. First, autoregressive,
          physiological and task-based noise types are mixed together in
          _generate_temporal_noise. The parameter values for 'auto_reg_sigma',
          'physiological_sigma' and 'task_sigma' describe the proportion of
          mixing of these elements, respectively. However critically, 'SFNR' is
          the parameter that controls how much noise these components contribute
          to the brain. 'auto_reg_rho' and 'ma_rho' set parameters for the
          autoregressive noise being simulated. Second, drift noise is added to
          this, according to the size of 'drift_sigma'. Thirdly, system noise is
          added based on the 'SNR' parameter. Finally, 'fwhm' is used to estimate
          the smoothness of the noise being inserted. If 'matched' is set to
          true, then it will fit the parameters to match the participant as best
          as possible.
  
          Variables defined as follows:
  
          snr [float]: Ratio of MR signal to the spatial noise
          sfnr [float]: Ratio of the MR signal to the temporal noise. This is the
          total variability that the following sigmas 'sum' to:
  
          task_sigma [float]: Size of the variance of task specific noise
          auto_reg_sigma [float]: Size of the variance of autoregressive
          noise. This is an ARMA process where the AR and MA components can be
          separately specified
          physiological_sigma [float]: Size of the variance of physiological
          noise
  
          drift_sigma [float]: Size of the variance of drift noise
  
          auto_reg_rho [list]: The coefficients of the autoregressive
          components you are modeling
          ma_rho [list]:The coefficients of the moving average components you
          are modeling
          max_activity [float]: The max value of the averaged brain in order
          to reference the template
          voxel_size [list]: The mm size of the voxels
          fwhm [float]: The gaussian smoothing kernel size (mm)
          matched [bool]: Specify whether you are fitting the noise parameters
  
          The volumes of brain noise that are generated have smoothness
          specified by 'fwhm'
  
      magnitude : list of floats
          This specifies the size, in terms of the metric choosen below,
          of the signal being generated. This can be a single number,
          and thus apply to all signal timecourses, or it can be array and
          thus different for each voxel.
  
      method : str
          Select the procedure used to calculate the signal magnitude,
          some of which are based on the definitions outlined in Welvaert &
          Rosseel (2013):
          - 'SFNR': Change proportional to the temporal variability,
          as represented by the (desired) SFNR
          - 'CNR_Amp/Noise-SD': Signal magnitude relative to the temporal
          noise
          - 'CNR_Amp2/Noise-Var_dB': Same as above but converted to decibels
          - 'CNR_Signal-SD/Noise-SD': Standard deviation in signal
          relative to standard deviation in noise
          - 'CNR_Signal-Var/Noise-Var_dB': Same as above but converted to
          decibels
          - 'PSC': Calculate the percent signal change based on the
          average activity of the noise (mean / 100 * magnitude)
  
  
      Returns
      ----------
      signal_function_scaled : 4d numpy array
          The new signal volume with the appropriately set signal change
  
      """
  
      # If you have only one magnitude value, duplicate the magnitude for each
      #  timecourse you have
      assert (
          type(magnitude) is list
      ), '"magnitude" should be a list of floats'
      if len(magnitude) == 1:
          magnitude *= signal_function.shape[1]
  
      # Check that the signal_function and noise_function are the same size
      if signal_function.shape != noise_function.shape:
          msg = "noise_function is not the same size as signal_function"
          raise ValueError(msg)
  
      # Scale all signals that to have a range of -1 to 1. This is
      # so that any values less than this will be scaled appropriately
      signal_function /= np.max(np.abs(signal_function))
  
      # Iterate through the timecourses and calculate the metric
      signal_function_scaled = np.zeros(signal_function.shape)
      for voxel_counter in range(signal_function.shape[1]):
          # Pull out the values for this voxel
          sig_voxel = signal_function[:, voxel_counter]
          noise_voxel = noise_function[:, voxel_counter]
          magnitude_voxel = magnitude[voxel_counter]
  
          # Calculate the maximum signal amplitude (likely to be 1,
          # but not necessarily)
          max_amp = np.max(np.abs(sig_voxel))
  
          # Calculate the scaled time course using the specified method
          if method == "SFNR":
              # How much temporal variation is there, relative to the mean
              # activity
              temporal_var = noise_voxel.mean() / noise_dict["sfnr"]
  
              # Multiply the timecourse by the variability metric
              new_sig = sig_voxel * (temporal_var * magnitude_voxel)
  
          elif method == "CNR_Amp/Noise-SD":
              # What is the standard deviation of the noise
              noise_std = np.std(noise_voxel)
  
              # Multiply the signal timecourse by the the CNR and noise (
              # rearranging eq.)
              new_sig = sig_voxel * (magnitude_voxel * noise_std)
  
          elif method == "CNR_Amp2/Noise-Var_dB":
              # What is the standard deviation of the noise
              noise_std = np.std(noise_voxel)
  
              # Rearrange the equation to compute the size of signal change in
              #  decibels
              scale = (
                  (10 ** (magnitude_voxel / 20)) * noise_std / max_amp
              )
  
              new_sig = sig_voxel * scale
  
          elif method == "CNR_Signal-SD/Noise-SD":
              # What is the standard deviation of the signal and noise
              sig_std = np.std(sig_voxel)
              noise_std = np.std(noise_voxel)
  
              # Multiply the signal timecourse by the the CNR and noise (
              # rearranging eq.)
              new_sig = sig_voxel * (
                  (magnitude_voxel / max_amp) * noise_std / sig_std
              )
  
          elif method == "CNR_Signal-Var/Noise-Var_dB":
              # What is the standard deviation of the signal and noise
              sig_std = np.std(sig_voxel)
              noise_std = np.std(noise_voxel)
  
              # Rearrange the equation to compute the size of signal change in
              #  decibels
              scale = (
                  (10 ** (magnitude_voxel / 20))
                  * noise_std
                  / (max_amp * sig_std)
              )
  
              new_sig = sig_voxel * scale
  
          elif method == "PSC":
              # What is the average activity divided by percentage
              scale = (noise_voxel.mean() / 100) * magnitude_voxel
              new_sig = sig_voxel * scale
  
          signal_function_scaled[:, voxel_counter] = new_sig
  
      # Return the scaled time course
      return signal_function_scaled
  
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 27 for calling a global function with example usage: compute_signal_change(signal_function, noise_function, noise_dict, magnitude, method) and returns: signal_function_scaled