randomise -i OneSamp4D -o OneSampT -1 -T
import matplotlib
import matplotlib.pylab as plt
import numpy as np
%matplotlib inline
import nibabel as nib
from nilearn.plotting import plot_stat_map, plot_roi
from scipy.spatial.distance import pdist
from scipy.stats import norm, t
import seaborn as sns

all_copes_file = "/home/vanessa/Desktop/tfMRI_LANGUAGE_STORY.nii_tstat1.nii.gz"
all_copes = nib.load(all_copes_file)
plot_stat_map(all_copes)
print "Here is our map created with randomise for all 486 subjects, for the story contrast"

# Function to flag outliers
def plot_outliers(image,n_std=6):
  mr = nib.load(image)
  data = mr.get_data()
  mean = data.mean()
  std = data.std()
  six_dev_up = mean + n_std * std
  six_dev_down = mean - n_std*std
  empty_brain = np.zeros(data.shape)
  empty_brain[data>=six_dev_up] = 1
  empty_brain[data<=six_dev_down] = 1 
  outlier_nii = nib.nifti1.Nifti1Image(empty_brain,affine=mr.get_affine(),header=mr.get_header())
  plot_roi(outlier_nii)

plot_outliers(all_copes_file)

data = all_copes.get_data()
data = data[data!=0]
sns.distplot(data.flatten(), label="Original T-Stat Data")
plt.legend()
print "Here is our map created with randomise for all 486 subjects, for the story contrast"

dof=486 - 2 
data = all_copes.get_data()
p_values = t.sf(data, df = dof)
p_values[p_values==1] = 0.99999999999999
sns.distplot(p_values.flatten(), label="P-Values from T-Stat Data")
plt.legend()
print "Here are the p-values created from the t-stat map, including all zeros in the map when we calculate"

z_values = norm.isf(p_values)
sns.distplot(z_values.flatten(), label="Z-Values from T-Stat Data")
plt.legend()
print "Here are the z-values created from the t-stat map, including all zeros in the map when we calculate"


# Need to make sure we look at the same slices :)
def plot_outliers(image,cut_coords,n_std=6):
  mr = nib.load(image)
  data = mr.get_data()
  mean = data.mean()
  std = data.std()
  six_dev_up = mean + n_std * std
  six_dev_down = mean - n_std*std
  empty_brain = np.zeros(data.shape)
  empty_brain[data>=six_dev_up] = 1
  empty_brain[data<=six_dev_down] = 1 
  outlier_nii = nib.nifti1.Nifti1Image(empty_brain,affine=mr.get_affine(),header=mr.get_header())
  plot_roi(outlier_nii,cut_coords=cut_coords)

Z_nii = nib.nifti1.Nifti1Image(z_values,affine=all_copes.get_affine(),header=all_copes.get_header())
nib.save(Z_nii,"/home/vanessa/Desktop/Zimage.nii")
plot_outliers("/home/vanessa/Desktop/Zimage.nii",cut_coords=(7,0,13))

data = all_copes.get_data()

# Let's select just the nonzero voxels
nonzero = data[data!=0]

# We will store our results here
Z = np.zeros(len(nonzero))

# Select values less than or == 0, and greater than zero
c  = np.zeros(len(nonzero))
k1 = (nonzero <= c)
k2 = (nonzero > c)

# Subset the data into two sets
t1 = nonzero[k1]
t2 = nonzero[k2]

# Calculate p values for <=0
p_values_t1 = t.cdf(t1, df = dof)
z_values_t1 = norm.ppf(p_values_t1)

# Calculate p values for > 0
p_values_t2 = t.cdf(-t2, df = dof)
z_values_t2 = -norm.ppf(p_values_t2)
Z[k1] = z_values_t1
Z[k2] = z_values_t2
sns.distplot(Z, label="Z-Values from T-Stat Data")
plt.legend()

empty_nii = np.zeros(all_copes.shape)
empty_nii[all_copes.get_data()!=0] = Z
Z_nii_fixed = nib.nifti1.Nifti1Image(empty_nii,affine=all_copes.get_affine(),header=all_copes.get_header())
nib.save(Z_nii_fixed,"/home/vanessa/Desktop/Zfixed.nii")
plot_stat_map(Z_nii_fixed)