# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

import nibabel as nib
import glob
import sys
sys.path.append('gradient_data/src/')
import os
import scipy.stats
import numpy as np
import statsmodels.sandbox.stats.multicomp
import nibabel as nib

# +
### Comparison of cerebello-cerebral connectivity between cerebellar Gradient 1 and 2 peaks at each area of
### motor and nonmotor representation was calculated as follows:
### 1) Using workbench view, save dscalar map of connectivity from peak of each Gradient (e.g. "L2midcog_fconn1003subjects.dscalar.nii")
### 2) Using workbench view, annotate correlation values between seeds:
###                           Fisher_z
### R1mot R2mot               0.231264
### L1mot L2mot               0.23966
### R12highcog R3highcog      0.313801
### L12highcog L3highcog      0.370525
### R1midcog R2midcog         0.450554
### L1midcog L2midcog         0.330115
### R1midcog R3midcog         0.330242
### L1midcog L3midcog         0.222103
### R1midcog R3midcog_alt     0.245531
### L1midcog L3midcog_alt     0.336851
### R2midcog R3midcog         0.276653
### L2midcog L3midcog         0.154236
### R2midcog R3midcog_alt     0.188767
### L2midcog L3midcog_alt     0.218249
### "Alt" refers to peak of Gradient 2 in lobule IX if peak was in lobule X, and vice versa. These values were
### not included in the final analysis
### 3) Contrast r values between seeds using the equations from Meng et al., 1992:
### wb_command -cifti-math '((var_zr1)-(var_zr2))*((sqrt(1003-3))/(2*(1-0.301)*((1-(((1-0.301)/(2*(1-(((tanh(var_zr1)^2)+(tanh(var_zr2)^2))/2)))+((1-0.301)/(2*(1-((( tanh(var_zr1)^2)+( tanh(var_zr2)^2))/2)))>1)*(1-(1-0.301)/(2*(1-(((tanh(var_zr1)^2)+(tanh(var_zr2)^2))/2))))))*((( tanh(var_zr1)^2)+( tanh(var_zr2)^2))/2))/(1-(((tanh(var_zr1)^2)+(tanh(var_zr2)^2))/2)))))' comparison_L2midcog_vs_L1midcog_1003subjects.dscalar.nii -var var_zr1 L2midcog_fconn1003subjects.dscalar.nii -var var_zr2 L1midcog_fconn1003subjects.dscalar.nii
### Substitute 0.301 with each particular correlation value between seeds. Note that z values have to be converted to r values (e.g. L1midcog L2midcog z=0.330115 corresponds to r=0.301)
### 4) Correct comparison maps with FDR as follows:
# -

directory = 'xxx' ### write directory with comparison maps here
files = os.listdir(directory)

### Check that files are loaded
for x in files:
    if x.endswith('motor.dscalar.nii'):
        print(x)

for x in files:
    if x.endswith('dscalar.nii'):
        print('calculating:    ' + x)
        zvalues = nib.load(x).get_data()
        pvals = scipy.stats.norm.sf(zvalues)
        pvals = pvals.T
        pvals_onlycortex = pvals[0:59411] ### This selects values only from the cortex
        pvals_onlycortex = pvals_onlycortex.T
        pvals_onlycortex_FDR = statsmodels.sandbox.stats.multicomp.multipletests(pvals_onlycortex[0], alpha=0.05, method='fdr_bh', is_sorted=False, returnsorted=False)
        ### Put cortical FDR corrected p values back to the original pvals matrix
        pvals_onlycortex_FDR[1].shape = (59411, 1)
        pvals[0:59411] = pvals_onlycortex_FDR[1]
        np.save(x + '_pvalues_onetailedcortexonly_FDR.npy', pvals.T[0])
        print('finished calculating:    ' + x)

### Transform to dscalar format
files = os.listdir(directory) #This is to update list of files in directory
for y in files:
    if y.endswith('.npy'):
        res = nib.load(directory2 + '/hcp.tmp.lh.dscalar.nii').get_data()
        cortL = np.squeeze(np.array(np.where(res != 0)[0], dtype=np.int32))
        res = nib.load(directory2 + '/hcp.tmp.rh.dscalar.nii').get_data()
        cortR = np.squeeze(np.array(np.where(res != 0)[0], dtype=np.int32))
        cortLen = len(cortL) + len(cortR)
        del res
        emb = np.load(y)
        emb.shape
        emb.shape = (91282, 1)
        tmp = nib.load('comparison_R1midcog_vs_R2midcog_peakfconn1003.dscalar.nii') #Just to take the shape; has to be dscalar with one map, and brain only
        tmp_cifti = nib.cifti2.load('comparison_R1midcog_vs_R2midcog_peakfconn1003.dscalar.nii')
        data = tmp_cifti.get_data() * 0
        mim = tmp.header.matrix[1]
        for idx, bm in enumerate(mim.brain_models):
            print ((idx, bm.index_offset, bm.brain_structure))
        img = nib.cifti2.Cifti2Image(emb.T, nib.cifti2.Cifti2Header(tmp.header.matrix))
        img.to_filename(y + '.dscalar.nii')

# +
### 5) Open corrected maps with wb_view and use a threshold of 0.05

StarcoderdataJupyterScriptsDedupFiltered