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, all noise types
        ending with sigma (e.g. motion sigma) are mixed together in
        _generate_temporal_noise. The sigma values describe the proportion of
        mixing of these elements. However critically, SFNR is the
        parameter that describes how much noise these components contribute
        to the brain. If you set the noise dict to matched then it will
        fit the parameters to match the participant as best as possible.

    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]

    # 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
def compute_signal_change ( signal_function , noise_function , noise_dict , magnitude , method = 'PSC' , ) : assert type ( magnitude ) is list , '"magnitude" should be a list of floats' if len ( magnitude ) == 1 : magnitude *= signal_function . shape [ 1 ] signal_function /= np . max ( np . abs ( signal_function ) ) signal_function_scaled = np . zeros ( signal_function . shape ) for voxel_counter in range ( signal_function . shape [ 1 ] ) : sig_voxel = signal_function [ : , voxel_counter ] noise_voxel = noise_function [ : , voxel_counter ] magnitude_voxel = magnitude [ voxel_counter ] max_amp = np . max ( np . abs ( sig_voxel ) ) if method == 'SFNR' : temporal_var = noise_voxel . mean ( ) / noise_dict [ 'sfnr' ] new_sig = sig_voxel * ( temporal_var * magnitude_voxel ) elif method == 'CNR_Amp/Noise-SD' : noise_std = np . std ( noise_voxel ) new_sig = sig_voxel * ( magnitude_voxel * noise_std ) elif method == 'CNR_Amp2/Noise-Var_dB' : noise_std = np . std ( noise_voxel ) scale = ( 10 ** ( magnitude_voxel / 20 ) ) * noise_std / max_amp new_sig = sig_voxel * scale elif method == 'CNR_Signal-SD/Noise-SD' : sig_std = np . std ( sig_voxel ) noise_std = np . std ( noise_voxel ) new_sig = sig_voxel * ( ( magnitude_voxel / max_amp ) * noise_std / sig_std ) elif method == 'CNR_Signal-Var/Noise-Var_dB' : sig_std = np . std ( sig_voxel ) noise_std = np . std ( noise_voxel ) scale = ( 10 ** ( magnitude_voxel / 20 ) ) * noise_std / ( max_amp * sig_std ) new_sig = sig_voxel * scale elif method == 'PSC' : scale = ( ( noise_voxel . mean ( ) / 100 ) * magnitude_voxel ) new_sig = sig_voxel * scale signal_function_scaled [ : , voxel_counter ] = new_sig return signal_function_scaled
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, all noise types
        ending with sigma (e.g. motion sigma) are mixed together in
        _generate_temporal_noise. The sigma values describe the proportion of
        mixing of these elements. However critically, SFNR is the
        parameter that describes how much noise these components contribute
        to the brain. If you set the noise dict to matched then it will
        fit the parameters to match the participant as best as possible.

    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