cifti analysis: broadcasting error

[mszinte]

Hi,

I have an error when running latest pybest with cifti files. Here is the command I use :

pybest /scratch/mszinte/data/PredictEye/deriv_data/fmriprep_new/fmriprep/ 
/scratch/mszinte/data/PredictEye/bids_data/
/scratch/mszinte/data/PredictEye/deriv_data/pybest_new/indices/cifti_indices.hdf5 
'Left_indices' 'Right_indices' 'Subcortex_indices' 
--out-dir /scratch/mszinte/data/PredictEye/deriv_data/pybest_new/ 
--subject '01' --space 'fsLR_den-170k' --iscifti 'y' --mode 'all' 
--noise-source fmriprep --skip-signalproc  --verbose 'DEBUG'

the code goes through the first steps but next reach a broadcasting error :

/home/mszinte/anaconda3/envs/idp/lib/python3.6/site-packages/nilearn/datasets/__init__.py:95: FutureWarning: Fetchers from the nilearn.datasets module will be updated in version 0.9 to return python strings instead of bytes and Pandas dataframes instead of Numpy arrays.
  "Numpy arrays.", FutureWarning)
2021-08-30 18:02 [INFO   ]  Using BIDS directory /scratch/mszinte/data/PredictEye/bids_data/
2021-08-30 18:02 [INFO   ]  No RETROICOR directory, so assuming no physio data.
2021-08-30 18:02 [WARNING]  No session identifier given; assuming a single session.
2021-08-30 18:02 [WARNING]  No single-trial-id found! Treating trials as conditions
2021-08-30 18:02 [INFO   ]  Found 3 session(s) for sub-01 ['01', '02', '03']
2021-08-30 18:02 [INFO   ]  Found 4 task(s) for sub-01 and ses-01 ['PurLoc', 'pMF', 'rest', 'SacLoc']
2021-08-30 18:02 [INFO   ]  Found 1 task(s) for sub-01 and ses-02 ['pRF']
2021-08-30 18:02 [INFO   ]  Found 2 task(s) for sub-01 and ses-03 ['PurVELoc', 'SacVELoc']
2021-08-30 18:02 [INFO   ]  Starting process for sub-01, ses-01, task-PurLoc
2021-08-30 18:02 [INFO   ]  Found 2 complete runs for task PurLoc
2021-08-30 18:02 [INFO   ]  Starting preprocessing of functional data ... 
2021-08-30 18:02 [INFO   ]  Preprocessing funcs:  ███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████  2/2
Traceback (most recent call last):
  File "/home/mszinte/anaconda3/envs/idp/bin/pybest", line 33, in <module>
    sys.exit(load_entry_point('pybest', 'console_scripts', 'pybest')())
  File "/home/mszinte/anaconda3/envs/idp/lib/python3.6/site-packages/click/core.py", line 1137, in __call__
    return self.main(*args, **kwargs)
  File "/home/mszinte/anaconda3/envs/idp/lib/python3.6/site-packages/click/core.py", line 1062, in main
    rv = self.invoke(ctx)
  File "/home/mszinte/anaconda3/envs/idp/lib/python3.6/site-packages/click/core.py", line 1404, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/home/mszinte/anaconda3/envs/idp/lib/python3.6/site-packages/click/core.py", line 763, in invoke
    return __callback(*args, **kwargs)
  File "/home/mszinte/softwares/pybest/pybest/cli.py", line 132, in main
    ddict = preprocess_funcs(ddict, cfg, logger)
  File "/home/mszinte/softwares/pybest/pybest/preproc.py", line 57, in preprocess_funcs
    save_data(data, cfg, ddict, par_dir='preproc', run=None, desc='preproc', dtype='bold')
  File "/home/mszinte/softwares/pybest/pybest/utils.py", line 283, in save_data
    subc_orig[pos] = subc_data
ValueError: shape mismatch: value array of shape (62053,416) could not be broadcast to indexing result of shape (62053,208)

indeed the TR indice is incorrect as it is multiplied by the number of run (208 x 2 = 416), any thought how to fix that ?

[lukassnoek]

Hey Martin, I haven't implemented the Cifti stuff. Perhaps @ecasimiro can take a look at this?

[ecasimiro]

Hey Martin, could you try adding something to the source code for me? Starting on line 272 of utils.py (below if cfg['cifti'] ...:)

runs = len(ddict['funcs'])
data_splitted = [data[:,x:x + runs] for x in range(0, data.shape[1], runs)]
full_data = dict()
for idx, run_data in enumerate(data_splitted):
    subc_data = cfg['subc_original']
    pos = cfg['pos']
    subc_data[pos] = run_data.T
    full_data[idx] = subc_data
    data_to_save = np.concatenate([v for k, v in full_data.items()], 3)
    np.save(f_out + '_subc.npy', data_to_save)

[mszinte]

Ok it was still incorrect, i'm working on a solution now, i will share it later... it might take some time as after fixing this error i received some other ones. Thanks.

[ecasimiro]

I was expecting more errors to pop up. If you need any help let me know, but then I'll need to look for a similar cifti file with multiple runs to test.

[mszinte]

Hi @ecasimiro you have written:

but after loading, 'tr' is in seconds, so it should be *=1000, isn't it ? Althought after doing that there is a warning showing up :

/home/mszinte/anaconda3/envs/idp/lib/python3.6/site-packages/nilearn/glm/first_level/design_matrix.py:108: UserWarning: High-pass filter will span all accessible frequencies and saturate the design matrix. You may want to reduce the high_pass value.The provided value is 0.01 Hz
  'The provided value is {0} Hz'.format(high_pass))

What do you think shoudl be the correct value ?

[ecasimiro]

Hey @mszinte, I think you're right that the cifti one does not need the /= 1000 and should just be the TR. I added that line because I was following the load_gifti as an example (that part was already there) to set up my code. The load_gifti one I think returns TRs as msec, so does need it.