In [3]:
from __future__ import division
import pandas as pd
import numpy as np
import os
import re
import copy
from pprint import pprint
from glob import glob
In [4]:
def create_localizer_block_single_effector(n_trials, response_modality='eye',
block_number=None, pseudorandomize=True,
add_timing=True, null_trials=0, TR=2, feedback=False):
# Only two trial types here: 'left' is correct, or 'right' is correct
trial_types = np.repeat([0, 1], repeats=n_trials/2) # left/right
# trial_types = np.hstack((np.repeat([0, 1], repeats=n_trials / 4), # Pro-saccadic cue, left/right corr
# np.repeat([2, 3], repeats=n_trials / 4))) # Anti-saccadic cue, left/right corr
# Initialize arrays
cue_by_trial = np.zeros(n_trials, dtype='<U5')
correct_answers = np.zeros(n_trials, dtype=np.int8)
# Define the cues for every trial
cue_by_trial[(trial_types == 0)] = 'LEFT'
cue_by_trial[(trial_types == 1)] = 'RIGHT'
# cue_by_trial[(trial_types == 2) | (trial_types == 3)] = 'anti'
# Define the responses ('correct answers')/directions for every trial
correct_answers[trial_types == 0] = 0
correct_answers[trial_types == 1] = 1
# correct_answers[(trial_types == 0) | (trial_types == 2)] = 0 # 0 = respond LEFT
# correct_answers[(trial_types == 1) | (trial_types == 3)] = 1 # 1 = respond RIGHT
# Create dataframe for easier handling
trial_data = pd.DataFrame({'correct_answer': correct_answers,
'cue': cue_by_trial,
'trial_type': trial_types})
# Should we pseudorandomize?
if pseudorandomize:
trial_data = Pseudorandomizer(trial_data, max_identical_iters={'correct_answer': 2,
'cue': 2}).run()
trial_data['null_trial'] = False
if null_trials > 0:
trial_data = add_pseudorandom_null_trials(trial_data, n_null_trials=null_trials)
# Add block number for completeness
if block_number is not None:
trial_data['block'] = block_number
trial_data['block_type'] = 'localizer'
# Usually, we also want to add the duration of all the 'trial phases'
if add_timing:
# Set phase_3 and phase_0 to 0 (no post-cue fixcross, no feedback)
if feedback:
phase_6 = 1
else:
phase_6 = None
trial_data = get_localizer_timing(trial_data, phase_6=phase_6, TR=TR)
trial_data['response_modality'] = response_modality.lower()
return trial_data
In [5]:
def create_cognitive_block(n_trials, block_number=None, response_modality=None,
add_timing=True, pseudorandomize=True, n_null_trials=0, TR=2, cue=None, **kwargs):
"""
Creates a block of SAT-trials; mixing speed and accuracy trials
"""
trial_types = np.hstack((np.repeat([0, 1], repeats=n_trials / 4), # SPEED cue, left/right corr
np.repeat([2, 3], repeats=n_trials / 4))) # ACCURACY cue, left/right corr
if trial_types.shape[0] != n_trials:
raise(ValueError('The provided n_trials (%d) could not be split into the correct number of trial types. '
'Closest option is %d trials' % (n_trials, trial_types.shape[0])))
cue_by_trial = np.zeros(n_trials, dtype='<U5')
correct_answers = np.zeros(n_trials, dtype=np.int8)
cue_by_trial[(trial_types == 0) | (trial_types == 1)] = 'SPD'
cue_by_trial[(trial_types == 2) | (trial_types == 3)] = 'ACC'
correct_answers[(trial_types == 0) | (trial_types == 2)] = 0 # 0 = left is correct
correct_answers[(trial_types == 1) | (trial_types == 3)] = 1 # 1 = right is correct
# Create dataframe for easier handling
trial_data = pd.DataFrame({'correct_answer': correct_answers,
'cue': cue_by_trial,
'trial_type': trial_types})
if pseudorandomize:
trial_data = Pseudorandomizer(trial_data,
max_identical_iters={'cue': 4, 'correct_answer': 4}).run()
if n_null_trials > 0:
trial_data['null_trial'] = False
trial_data = add_pseudorandom_null_trials(trial_data,
n_null_trials=n_null_trials,
null_column_name='null_trial')
if block_number is not None:
trial_data['block'] = block_number
if response_modality is not None:
trial_data['response_modality'] = response_modality
trial_data['block_type'] = 'cognitive_%s' % response_modality
if add_timing:
trial_data = get_block_timing(trial_data, TR=TR, **kwargs) # Add default timing
if cue is not None:
trial_data['cue'] = cue
return trial_data
In [34]:
def create_limbic_block(n_trials, subject_number=1, block_number=None,
response_modality=None, add_timing=True, pseudorandomize=True,
n_null_trials=0, TR=2, **kwargs):
trial_types = np.hstack((np.repeat([0, 1], repeats=n_trials/4), # Neutral cue, left/right corr
np.repeat([2, 3], repeats=n_trials/8), # Left cue, left/right corr
np.repeat([4, 5], repeats=n_trials/8))) # Right cue, left/right corr
if 'random_sample_trials' in kwargs:
# Randomly sample some trials from the created array.
trial_types = trial_types[np.random.choice(a=[0, 2, 5, 4], size=kwargs['random_sample_trials'])]
n_trials = trial_types.shape[0]
kwargs.pop('random_sample_trials')
else:
if trial_types.shape[0] != n_trials:
raise(ValueError('The provided n_trials (%d) could not be split into the correct number of trial types. '
'Closest option is %d trials' % (n_trials, trial_types.shape[0])))
cue_by_trial = np.zeros(n_trials, dtype='<U5')
correct_answers = np.zeros(n_trials, dtype=np.int8)
cue_by_trial[(trial_types == 0) | (trial_types == 1)] = 'NEU'
cue_by_trial[(trial_types == 2) | (trial_types == 3)] = 'LEFT'
cue_by_trial[(trial_types == 4) | (trial_types == 5)] = 'RIGHT'
correct_answers[(trial_types == 0) |
(trial_types == 2) |
(trial_types == 4)] = 0 # 0 = left is correct
correct_answers[(trial_types == 1) |
(trial_types == 3) |
(trial_types == 5)] = 1 # 1 = right is correct
# Create dataframe for easier handling
trial_data = pd.DataFrame({'correct_answer': correct_answers,
'cue': cue_by_trial,
'trial_type': trial_types})
if pseudorandomize:
trial_data = Pseudorandomizer(trial_data,
max_identical_iters={'cue': 4, 'correct_answer': 4}).run()
if n_null_trials > 0:
trial_data['null_trial'] = False
trial_data = add_pseudorandom_null_trials(trial_data,
n_null_trials=n_null_trials,
null_column_name='null_trial')
if block_number is not None:
trial_data['block'] = block_number
if response_modality is not None:
trial_data['response_modality'] = response_modality
trial_data['block_type'] = 'limbic_%s' % response_modality
if add_timing:
trial_data = get_block_timing(trial_data, TR=TR, **kwargs) # Add default timing
return trial_data
In [35]:
def get_localizer_timing(trial_data, phase_0=None, phase_1=None, phase_2=None, phase_3=None, phase_4=None, phase_5=None, phase_6=None, TR=2):
"""
Each localizer trial consists of 7 phases.
In phase_0, we wait for the scanner pulse. Note that phase_0 of trial n is the ITI after trial n-1. Set this timing always to 0: it is the `minimum` time to wait for the pulse
In phase_1, we show the pre-cue fixation cross. By default, timing is jittered (0s, 0.5s, 1s, 1.5s if TR=2 -- 0, .75, 1.5, 2.25 if TR=3).
In phase_2, we show the cue. Defaults to 1.5 second.
In phase_3, we show the post-cue fixation cross. Defaults to 0s.
In phase_4, we assume the participant responds, and wait a bit until we show the fix cross. Defaults to 0.6s
In phase_5 and phase_6, we do nothing (exist for compatibility with the experimental blocks)
Phase_7 is ITI
"""
if TR == 2:
trial_data['phase_0'] = 0 if phase_0 is None else phase_0
trial_data['phase_1'] = np.random.choice([0.2, .7, 1.2, 1.7], size=trial_data.shape[0]) if phase_1 is None else phase_1
trial_data['phase_2'] = 0.8 if phase_2 is None else phase_2
trial_data['phase_3'] = 0 if phase_3 is None else phase_3
trial_data['phase_4'] = 0.6 if phase_4 is None else phase_4
trial_data['phase_5'] = 0 if phase_5 is None else phase_5
trial_data['phase_6'] = 0 if phase_6 is None else phase_6
elif TR == 3:
trial_data['phase_0'] = 0 if phase_0 is None else phase_0
trial_data['phase_1'] = np.random.choice([0, .750, 1.500, 2.250], size=trial_data.shape[0]) if phase_1 is None else phase_1
trial_data['phase_2'] = np.round(np.random.exponential(scale=1/6, size=trial_data.shape[0])+.7, 3) if phase_2 is None else phase_2
trial_data['phase_3'] = 0 if phase_3 is None else phase_3
trial_data['phase_4'] = 0.8 if phase_4 is None else phase_4
trial_data['phase_5'] = 0 if phase_5 is None else phase_5
trial_data['phase_6'] = 0 if phase_6 is None else phase_6
# Calculate duration of trial (depends on random, jittered durations of the fix cross)
trial_data['trial_duration'] = trial_data[['phase_' + str(x) for x in range(7)]].sum(axis=1)
# Because of TR = 2s, some trials can last 8 seconds, but most will last 10. Find trials with total time < 8 seconds
# We calculate the ITI as the difference between the minimum number of pulses necessary for all phases to show.
# [same applies to TR = 3]
# min_TRs = np.ceil(trial_data['trial_duration'].values / TR)
# trial_data['phase_7'] = min_TRs*TR - trial_data['trial_duration'].values
trial_data['phase_7'] = 6 - trial_data['trial_duration'].values
# Recalculate trial duration so it includes the ITI
trial_data['trial_duration'] = trial_data[['phase_' + str(x) for x in range(8)]].sum(axis=1)
# Add trial start times relative to start of block
trial_data['trial_start_time_block'] = trial_data['trial_duration'].shift(1).cumsum()
trial_data.loc[0, 'trial_start_time_block'] = 0
# Add cue onset times relative to start of block
trial_data['cue_onset_time_block'] = trial_data['trial_start_time_block'] + \
trial_data['phase_1']
# Add stimulus onset times relative to start of block
trial_data['stimulus_onset_time_block'] = trial_data['trial_start_time_block'] + \
trial_data['phase_1'] + \
trial_data['phase_2'] + \
trial_data['phase_3']
return trial_data
In [36]:
def get_block_timing(trial_data, phase_0=None, phase_1=None, phase_2=None, phase_3=None, phase_4=None, phase_5=None, phase_6=None, TR=2):
"""
Each trial consists of 7 phases.
In phase_0, we wait for the scanner pulse. Note that phase_0 of trial n is the ITI after trial n-1. Set this timing always to 0: it is the `minimum` time to wait for the pulse
In phase_1, we show the pre-cue fixation cross. By default, timing is jittered (0s, 0.5s, 1s, 1.5s)
In phase_2, we show the cue. In decision-making trials, this is 4.8 seconds.
In phase_3, we show the post-cue fixation cross. Timing is jittered (0s, 0.5s, 1s, 1.5s)
In phase_4, we show the stimulus. Default is 1.5s.
Phase 5 is defined as the period of stimulus presentation, after the participant made a response. The duration is determined by the participant RT, so not set here.
In phase_6, we show feedback. Default is 0.35s.
"""
if TR == 2:
trial_data['phase_0'] = 0 if phase_0 is None else phase_0
trial_data['phase_1'] = np.random.choice([0, .5, 1, 1.5], size=trial_data.shape[0]) if phase_1 is None else phase_1
trial_data['phase_2'] = 4.8 if phase_2 is None else phase_2
trial_data['phase_3'] = np.random.choice([0, .5, 1, 1.5], size=trial_data.shape[0]) if phase_3 is None else phase_3
trial_data['phase_4'] = 1.5 if phase_4 is None else phase_4
trial_data['phase_5'] = 0 if phase_5 is None else phase_5
trial_data['phase_6'] = 0.35 if phase_6 is None else phase_6
elif TR == 3:
trial_data['phase_0'] = 0 if phase_0 is None else phase_0
trial_data['phase_1'] = np.random.choice([0, .750, 1.500, 2.250], size=trial_data.shape[0]) if phase_1 is None else phase_1
trial_data['phase_2'] = 1 if phase_2 is None else phase_2
trial_data['phase_3'] = np.random.choice([.75, 1.50, 2.250, 3.000], size=trial_data.shape[0]) if phase_3 is None else phase_3
trial_data['phase_4'] = 2. if phase_4 is None else phase_4
trial_data['phase_5'] = 0 if phase_5 is None else phase_5
trial_data['phase_6'] = 0.5 if phase_6 is None else phase_6
elif TR == 0:
trial_data['phase_0'] = 0 if phase_0 is None else phase_0
trial_data['phase_1'] = np.random.choice([.750, 1.500, 2.250], size=trial_data.shape[0]) if phase_1 is None else phase_1
trial_data['phase_2'] = 1 if phase_2 is None else phase_2
trial_data['phase_3'] = np.random.choice([.750, 1.500, 2.250, 3.000], size=trial_data.shape[0]) if phase_3 is None else phase_3
trial_data['phase_4'] = 2. if phase_4 is None else phase_4
trial_data['phase_5'] = 0 if phase_5 is None else phase_5
trial_data['phase_6'] = 0.5 if phase_6 is None else phase_6
# Calculate duration of trial (depends on random, jittered durations of the fix cross)
trial_data['trial_duration'] = trial_data[['phase_' + str(x) for x in range(7)]].sum(axis=1)
if TR == 2:
# Because of TR = 2s, some trials can last 8 seconds, but most will last 10. Find trials with total time < 8 seconds
# We calculate the ITI as the difference between the minimum number of pulses necessary for all phases to show.
min_TRs = np.ceil(trial_data['trial_duration'].values / TR)
trial_data['phase_7'] = min_TRs*TR - trial_data['trial_duration'].values
elif TR == 3:
# In this case, fill all trials until 9s have passed.
trial_data['phase_7'] = 9 - trial_data['trial_duration'].values
elif TR == 0:
trial_data['phase_7'] = 0
# Recalculate trial duration so it includes the ITI
trial_data['trial_duration'] = trial_data[['phase_' + str(x) for x in range(8)]].sum(axis=1)
# Add trial start times relative to start of block
trial_data['trial_start_time_block'] = trial_data['trial_duration'].shift(1).cumsum()
trial_data.loc[0, 'trial_start_time_block'] = 0
# Add cue onset times relative to start of block
trial_data['cue_onset_time_block'] = trial_data['trial_start_time_block'] + \
trial_data['phase_1']
# Add stimulus onset times relative to start of block
trial_data['stimulus_onset_time_block'] = trial_data['trial_start_time_block'] + \
trial_data['phase_1'] + \
trial_data['phase_2'] + \
trial_data['phase_3']
return trial_data
In [37]:
class Pseudorandomizer(object):
def __init__(self, data, max_identical_iters={'cue': 4, 'correct_answer': 4}):
self.data = data
self.max_identical_iters = {x: y+1 for x, y in max_identical_iters.items()}
# add 1: if 4 rows is allowed, only give an error after 5 identical rows
def check_trial_rows(self, data, row_n):
"""
Returns True if any of the conditions for pseudorandomization are violated for the given rows,
False if they are fine.
"""
# First, check for the maximum iterations
for column, max_iter in self.max_identical_iters.items():
if row_n - max_iter < 0:
continue
# Select rows [max_iter-1 - row_n] we're going to check. Never select any row with index < 0
row_selection = [x for x in np.arange(row_n, row_n-max_iter, -1)]
# Next, we check if the selected rows only contain *1* trial type.
# If so, this means we have max_iter rows of the same trials, and we need to change something.
if data.iloc[row_selection][column].nunique() == 1:
return True
return False
def run(self, debug=False):
"""
Pseudorandomizes: makes sure that it is not possible to have more than x iterations for every type of column, specified in columns.
"""
# Start by copying from original data, and shuffle
self.data = self.data.sample(frac=1,
random_state=np.random.randint(0, 1e7, dtype='int')).reset_index(drop=True)
if debug:
outer_while_i = 0
debug_print_after_i = 100
good_set = False
while not good_set:
if debug:
outer_while_i += 1
reshuffle = False # Assume the dataset does not need reshuffling.
for row_n in range(0, self.data.shape[0]):
# Select rows [max_iter-1 - row_n] we're going to check
# Check if the current row, and the (max_iters-1) rows before, are the same value (number of unique values = 1).
# If so, then move the current row number to the bottom of the dataframe. However, we need to re-check the same four rows again
# after moving a row to the bottom: therefore, a while loop is necessary.
checked_row = False
n_attempts_at_moving = 0
if debug:
inner_while_i = 0
while not checked_row:
if debug:
inner_while_i += 1
if inner_while_i > debug_print_after_i:
print('New inner loop started for current row')
if self.check_trial_rows(self.data, row_n):
if debug and inner_while_i > debug_print_after_i:
print('Found too many consecutively identical rows.')
# If there are too many consecutively identical rows at the bottom of the dataframe,
# break and start over/shuffle
if row_n >= (self.data.shape[0] - self.max_identical_iters[self.max_identical_iters.keys()[0]]):
if debug and inner_while_i > debug_print_after_i:
print('These occurred at row_n %d, which is at the bottom of the DF.' % row_n)
checked_row = True
reshuffle = True
# Too many consecutive identical rows? Move row_n to the bottom, and check again with the new row_n.
else:
if debug and inner_while_i > debug_print_after_i:
print('These occurred at row_n %d. Checking the remainder of the DF.' % row_n)
# Check if moving to the bottom even makes sense: if all remaining values are identical, it doesn't.
if (self.data.iloc[row_n:][self.max_identical_iters.keys()].nunique().values < 2).any():
if debug and inner_while_i > debug_print_after_i:
print('All remaining values are identical. I should stop the for-loop, and start over.')
checked_row = True
reshuffle = True
else:
if n_attempts_at_moving < 50:
n_attempts_at_moving += 1
if debug and inner_while_i > debug_print_after_i:
print('Not all remaining values are identical. I should move the final part to the bottom.')
# If not, move the current row to the bottom
row_to_move = self.data.iloc[row_n,:]
# Delete row from df
self.data.drop(row_n, axis=0, inplace=True)
# Append original row to end. Make sure to reset index
self.data = self.data.append(row_to_move).reset_index(drop=True)
# If we already tried moving the current row to the bottom for 50 times, let's forget about it and restart
else:
checked_row = True
reshuffle = True
else:
if debug and inner_while_i > debug_print_after_i:
print('Checked row, but the row is fine. Next row.')
checked_row = True
if reshuffle:
good_set = False
break # out of the for loop
# Reached the bottom of the dataframe, but no reshuffle call? Then we're set.
if row_n == self.data.shape[0]-1:
good_set = True
if reshuffle:
# Shuffle, reset index to ensure trial_data.drop(row_n) rows
self.data = self.data.sample(frac=1, random_state=np.random.randint(0, 1e7, dtype='int')).reset_index(drop=True)
return self.data
def add_pseudorandom_null_trials(data, min_row=4, max_row=4, min_n_rows_separate=7,
n_null_trials=10, null_column_name=''):
"""
Adds null trials interspersed at pseudorandom locations. You can determine the minimum
number of trials at the start before a null trial, the minimum number of trials at the end in which no
nulls are shown, and the minimum number of trials that the null trials have to be separated
"""
good_idx = False
while not good_idx:
indx = np.random.choice(np.arange(min_row, data.shape[0]-max_row),
replace=False, size=n_null_trials)
diffs = np.diff(np.sort(indx))
if (diffs >= min_n_rows_separate).all():
good_idx = True
data.index = np.setdiff1d(np.arange(data.shape[0] + n_null_trials), indx)
new_rows = pd.DataFrame({null_column_name: [True]*n_null_trials}, columns=data.columns, index=indx)
data = data.append(new_rows).sort_index()
return data
In [38]:
os.chdir('../designs')
n_trials_per_localizer_block = 6
n_trials_limbic = 8
n_trials_cognitive = 8
TR = 3
# Empty DataFrame
practice_design = pd.DataFrame()
# Start with localizer: 5 trials manual first?
loc_block = create_localizer_block_single_effector(n_trials=n_trials_per_localizer_block,
response_modality='hand',
block_number=0,
pseudorandomize=True, TR=TR, feedback=True)
practice_design = practice_design.append(loc_block)
# Then 5 trials without feedback
loc_block = create_localizer_block_single_effector(n_trials=n_trials_per_localizer_block,
response_modality='hand',
block_number=1,
pseudorandomize=True, TR=TR)
practice_design = practice_design.append(loc_block)
# Then 5 trials with eyes
loc_block = create_localizer_block_single_effector(n_trials=n_trials_per_localizer_block,
response_modality='eye',
block_number=2,
pseudorandomize=True, TR=TR)
practice_design = practice_design.append(loc_block)
# Localizer done?
# Flashtask: start with SAT, hand. In this first block, ignore cue; and speed-up everything
block_data = create_cognitive_block(n_trials=8,
block_number=3,
response_modality='hand',
n_null_trials=0,
TR=0, phase_2=0, phase_3=0, cue='ACC')
practice_design = practice_design.append(block_data)
# Now, introduce SAT-cue
block_data = create_cognitive_block(n_trials=n_trials_cognitive,
block_number=4,
response_modality='hand',
n_null_trials=0,
TR=TR)
practice_design = practice_design.append(block_data)
# Then use eye. Respond with mouse in practice session. No feedback
block_data = create_cognitive_block(n_trials=4,
block_number=5,
response_modality='eye',
n_null_trials=0,
TR=TR, phase_6=0)
practice_design = practice_design.append(block_data)
# Next, limbic block?
block_data = create_limbic_block(n_trials=n_trials_limbic,
block_number=6,
response_modality='hand',
n_null_trials=0,
TR=TR)
practice_design = practice_design.append(block_data)
# Finally, limbic with eye. No feedback
block_data = create_limbic_block(n_trials=8,
block_number=7,
response_modality='eye',
n_null_trials=0,
TR=TR, phase_6=0,
random_sample_trials=4)
practice_design = practice_design.append(block_data)
# Set indices
practice_design.index.name = 'block_trial_ID'
practice_design.reset_index(inplace=True)
practice_design.index.name = 'trial_ID'
# Add trial start times (relative to start of experiment)
practice_design['trial_start_time'] = practice_design['trial_duration'].shift(1).cumsum()
practice_design.loc[0, 'trial_start_time'] = 0
# Add cue onset times (relative to start of experiment)
practice_design['cue_onset_time'] = practice_design['trial_start_time'] + \
practice_design['phase_1']
# Add stimulus onset times (relative to start of experiment)
practice_design['stimulus_onset_time'] = practice_design['trial_start_time'] + \
practice_design['phase_1'] + \
practice_design['phase_2'] + \
practice_design['phase_3']
# Re-order column order for nicety
practice_design = practice_design[['block_trial_ID', 'block', 'block_type', 'null_trial', 'correct_answer', 'cue', 'response_modality', 'trial_type',
'phase_0', 'phase_1', 'phase_2', 'phase_3', 'phase_4', 'phase_5', 'phase_6', 'phase_7', 'trial_duration',
'trial_start_time', 'cue_onset_time', 'stimulus_onset_time',
'trial_start_time_block', 'cue_onset_time_block', 'stimulus_onset_time_block']]
practice_design['null_trial'] = False
# Save full data
if not os.path.exists(os.path.join('practice', 'all_blocks')):
os.makedirs(os.path.join('practice', 'all_blocks'))
practice_design.to_csv(os.path.join('practice', 'all_blocks', 'trials.csv'), index=True)
# Save individual blocks
for block_num in practice_design['block'].unique():
block_type = practice_design.loc[practice_design['block'] == block_num, 'block_type'].values[0]
block = practice_design.loc[practice_design['block'] == block_num]
if not os.path.exists(os.path.join('practice', 'block_%d_type_%s' % (block_num, block_type))):
os.makedirs(os.path.join('practice', 'block_%d_type_%s' % (block_num, block_type)))
block.to_csv(os.path.join('practice', 'block_%d_type_%s' % (block_num, block_type), 'trials.csv'), index=True)
In [ ]: