def ica_vis(subj_num):
  # Use the mean as a background
  mean_img_1 = image.mean_img(BOLD_file_1)
  mean_img_2 = image.mean_img(BOLD_file_2)
  mean_img_3 = image.mean_img(BOLD_file_3)

  plot_stat_map(image.index_img(component_img_1, 5), mean_img_1, output_file=os.path.join(data_path,'sub'+subj_num+'_BOLD','task001_run001'+'ica_1'+'.jpg'))
  plot_stat_map(image.index_img(component_img_1, 12), mean_img_1, output_file=os.path.join(data_path,'sub'+subj_num+'_BOLD','task001_run001'+'ica_2'+'.jpg'))

  plot_stat_map(image.index_img(component_img_2, 5), mean_img_2, output_file=os.path.join(data_path,'sub'+subj_num+'_BOLD','task002_run001'+'ica_1'+'.jpg'))
  plot_stat_map(image.index_img(component_img_2, 12), mean_img_2, output_file=os.path.join(data_path,'sub'+subj_num+'_BOLD','task002_run001'+'ica_2'+'.jpg'))

    def _run_interface(self, runtime):
        import matplotlib

        matplotlib.use("Agg")
        import pylab as plt

        wra_img = nib.load(self.inputs.wra_img)
        canonical_img = nib.load(self.inputs.canonical_img)
        title = self.inputs.title
        mean_wraimg = image.mean_img(wra_img)

        if title != "":
            filename = title.replace(" ", "_") + ".pdf"
        else:
            filename = "plot.pdf"

        fig = plotting.plot_anat(
            mean_wraimg, title="wrafunc & canonical single subject", cut_coords=range(-40, 40, 10), display_mode="z"
        )
        fig.add_edges(canonical_img)
        fig.savefig(filename)
        fig.close()

        self._plot = filename

        runtime.returncode = 0
        return runtime

def compute_subj_mask(subj):
    num = 0

    print 'Loading BOLD subject #%i...' % subj
    bold = load(path.join(data_dir, 'sub%03d/BOLD/task001_run00%i/bold_dico.nii.gz' % (subj, num+1)))

    print 'Averaging BOLD...'
    avg = image.mean_img(bold)

    print 'Saving average EPI...'
    avg_path = path.join(mask_dir, 'mean_sub%03d_run00%i.nii.gz' % (subj, num+1))
    affine = avg.get_affine()
    save(avg, avg_path)

    print 'Aligning the surface...'
    cortex.align.automatic('sub%03d' % subj, 'auto01', avg_path)

    m_thick = cortex.db.get_mask('sub%03d' % subj, 'auto01', type='thick')
    print 'Thick mask is computed: %i voxels' % m_thick.sum()
    m_thin = cortex.db.get_mask('sub%03d' % subj, 'auto01', type='thin')
    print 'Thin mask is computed: %i voxels' % m_thin.sum()

# .T is important here, because pycortex and nibabel have different
# assumptions about the axis order
    thick = Nifti1Image(np.float32(m_thick.T), affine)
    thin = Nifti1Image(np.float32(m_thin.T), affine)

    save(thick, path.join(mask_dir, 'sub%03d_mask_thick.nii.gz' % subj))
    save(thin, path.join(mask_dir, 'sub%03d_mask_thin.nii.gz' % subj))

	def _run_interface(self, runtime):
		import matplotlib
		matplotlib.use('Agg')
		import nibabel as nib
		from nilearn import plotting, datasets, image
		from nipype.interfaces.base import isdefined
		import numpy as np
		import pylab as plt
		import os

		wra_img = nib.load(self.inputs.wra_img)
		canonical_img = nib.load(self.inputs.canonical_img)
		title = self.inputs.title
		mean_wraimg = image.mean_img(wra_img)

		if title != "":
			filename = title.replace(" ", "_")+".pdf"
		else:
			filename = "plot.pdf"

		fig = plotting.plot_anat(mean_wraimg, title="wrafunc & canonical single subject", cut_coords=range(-40, 40, 10), display_mode='z')
		fig.add_edges(canonical_img)     
		fig.savefig(filename)
		fig.close()

		self._plot = filename

		runtime.returncode=0
		return runtime

def p_map(task, run, p_values_3d, threshold=0.05):
    """
    Generate three thresholded p-value maps.

    Parameters
    ----------
    task: int
        Task number
    run: int
        Run number
    p_value_3d: 3D array of p_value.
    threshold: The cutoff value to determine significant voxels.

    Returns
    -------
    threshold p-value images
    """
    fmri_img = image.smooth_img('../../../data/sub001/BOLD/' + 'task00' +
                                str(task) + '_run00' + str(run) +
                                '/filtered_func_data_mni.nii.gz',
                                fwhm=6)

    mean_img = image.mean_img(fmri_img)

    log_p_values = -np.log10(p_values_3d)
    log_p_values[np.isnan(log_p_values)] = 0.
    log_p_values[log_p_values > 10.] = 10.
    log_p_values[log_p_values < -np.log10(threshold)] = 0
    plot_stat_map(nib.Nifti1Image(log_p_values, fmri_img.get_affine()),
                  mean_img, title="Thresholded p-values",
                  annotate=False, colorbar=True)

def spikes_mask(in_file, in_mask=None, out_file=None):
    """
    Utility function to calculate a mask in which check
    for :abbr:`EM (electromagnetic)` spikes.
    """

    import os.path as op
    import nibabel as nb
    import numpy as np
    from nilearn.image import mean_img
    from nilearn.plotting import plot_roi
    from scipy import ndimage as nd

    if out_file is None:
        fname, ext = op.splitext(op.basename(in_file))
        if ext == '.gz':
            fname, ext2 = op.splitext(fname)
            ext = ext2 + ext
        out_file = op.abspath('{}_spmask{}'.format(fname, ext))
        out_plot = op.abspath('{}_spmask.pdf'.format(fname))

    in_4d_nii = nb.load(in_file)
    orientation = nb.aff2axcodes(in_4d_nii.affine)

    if in_mask:
        mask_data = nb.load(in_mask).get_data()
        a = np.where(mask_data != 0)
        bbox = np.max(a[0]) - np.min(a[0]), np.max(a[1]) - \
            np.min(a[1]), np.max(a[2]) - np.min(a[2])
        longest_axis = np.argmax(bbox)

        # Input here is a binarized and intersected mask data from previous section
        dil_mask = nd.binary_dilation(
            mask_data, iterations=int(mask_data.shape[longest_axis] / 9))

        rep = list(mask_data.shape)
        rep[longest_axis] = -1
        new_mask_2d = dil_mask.max(axis=longest_axis).reshape(rep)

        rep = [1, 1, 1]
        rep[longest_axis] = mask_data.shape[longest_axis]
        new_mask_3d = np.logical_not(np.tile(new_mask_2d, rep))
    else:
        new_mask_3d = np.zeros(in_4d_nii.shape[:3]) == 1

    if orientation[0] in ['L', 'R']:
        new_mask_3d[0:2, :, :] = True
        new_mask_3d[-3:-1, :, :] = True
    else:
        new_mask_3d[:, 0:2, :] = True
        new_mask_3d[:, -3:-1, :] = True

    mask_nii = nb.Nifti1Image(new_mask_3d.astype(np.uint8), in_4d_nii.get_affine(),
                              in_4d_nii.get_header())
    mask_nii.to_filename(out_file)

    plot_roi(mask_nii, mean_img(in_4d_nii), output_file=out_plot)
    return out_file, out_plot

def mean_img(in_file, out_file=None):
    """ Use nilearn.image.mean_img.
    Returns
    -------
    out_file: str
        The absolute path to the output file.
    """
    import nilearn.image as niimg
    return niimg.mean_img(in_file)

 def mean_nodiff_fct(in_file, bval_file, nodiff_b=0):
     import os, numpy as np
     from nilearn import image
     bvals = np.loadtxt(bval_file)
     nodiff_index = bvals <= nodiff_b
     nodiff_img = image.index_img(in_file, nodiff_index)
     mean_nodiff_img = image.mean_img(nodiff_img)
     out_file = os.path.abspath('mean_nodiff.nii.gz')
     mean_nodiff_img.to_filename(out_file)
     return out_file

def plot_segmentation(
        img, gm_filename, wm_filename=None, csf_filename=None,
        output_filename=None, cut_coords=None, display_mode='ortho',
        cmap=None, title='GM + WM + CSF segmentation', close=False):
    """
    Plot a contour mapping of the GM, WM, and CSF of a subject's anatomical.

    Parameters
    ----------
    img_filename: string or image object
                  path of file containing image data, or image object simply

    gm_filename: string
                 path of file containing Grey Matter template

    wm_filename: string (optional)
                 path of file containing White Matter template

    csf_filename: string (optional)
                 path of file containing Cerebro-Spinal Fluid template


    """
    # misc
    if cmap is None:
        cmap = plt.cm.gray
    if cut_coords is None:
        cut_coords = (-10, -28, 17)
    if display_mode in ['x', 'y', 'z']:
        cut_coords = (cut_coords['xyz'.index(display_mode)],)

    # plot img
    img = mean_img(img)
    img = reorder_img(img, resample="continuous")
    _slicer = plot_img(img, cut_coords=cut_coords, display_mode=display_mode,
                       cmap=cmap, black_bg=True)

    # add TPM contours
    gm = nibabel.load(gm_filename)
    _slicer.add_contours(gm, levels=[.51], colors=["r"])
    if not wm_filename is None:
        _slicer.add_contours(wm_filename, levels=[.51], colors=["g"])
    if not csf_filename is None:
        _slicer.add_contours(csf_filename, levels=[.51], colors=['b'])

    # misc
    _slicer.title(title, size=12, color='w', alpha=0)
    if not output_filename is None:
        plt.savefig(output_filename, bbox_inches='tight', dpi=200,
                    facecolor="k",
                    edgecolor="k")
        if close:
            plt.close()

def loader(anat, downsample, target_affine, dataroot, subject, maskpath, nrun,
           niifilename, labels, **kwargs):
    ''' 
    All parameters are submitted as cfg dictionary.
    Given parameters in cfg, return masked and concatenated over runs data 
    
    Input
    anat: MNI template
    downsample: 1 or 0
    target_affine: downsampling matrix
    dataroot: element of path to data
    subject: folder in dataroot with subject data
    maskpath: path to mask
    nrun: number of runs
    niifilename: how is the data file called
    labels: labels from load_labels function
    
    Output
    dict(nii_func=nii_func,nii_mean=nii_mean,masker=masker,nii_mask=nii_mask)
    nii_func: 4D data
    nii_mean: mean over 4th dimension
    masker: masker object from nibabel
    nii_mask: 3D mask
    '''
    nii_func = list()
    for r in range(nrun):
        fname = '{0}/{1}/run{2}/{3}'.format(dataroot, subject, r+1, niifilename) # Assumption about file location
        nii_img = load(fname, mmap=False)
        nii_img.set_sform(anat.get_sform())
        # Get mean over 4D
        nii_mean = mean_img(nii_img)
        # Masking
        nii_mask = load(maskpath)
        nii_mask.set_sform(anat.get_sform())
        # Binarize the mask
        nii_mask = check_binary(nii_mask)
        if downsample:
            nii_img = resample_img(nii_img, target_affine=target_affine)
            nii_mask = resample_img(nii_mask, target_affine=target_affine, interpolation='nearest')
        masker = NiftiMasker(nii_mask, standardize=True)
        nii_img = masker.fit_transform(nii_img)
        # Drop zero timepoints, zscore
        nii_img = drop_labels(nii_img, labels.get('to_drop_zeros')[r])
        nii_func.append(stats.zscore(nii_img, axis=0)) # zscore over time
    # throw data together
    nii_func = np.concatenate(nii_func)
    return dict(nii_func=nii_func, nii_mean=nii_mean, masker=masker, nii_mask=nii_mask)

def multi_session_time_slice_diffs(img_list):
    """ time slice difference on several 4D images

    Parameters
    ----------
    img_list: list of 4D Niimg-like
        Input multi-session images

    returns
    -------
    results : dict
        see time_slice_diffs docstring for details.

    note
    ----
    The results are accumulated across sessions
    """
    results = {}
    for i, img in enumerate(img_list):
        results_ = time_slice_diffs(img)
        if i == 0:
            for key, val in results_.items():
                # special case for 'session_length' to make
                # aggregation easier later on
                results[key] = val if key != 'session_length' else [val]
        else:
            results['volume_mean_diff2'] = np.hstack((
                    results['volume_mean_diff2'],
                    results_['volume_mean_diff2']))
            results['slice_mean_diff2'] = np.vstack((
                    results['slice_mean_diff2'],
                    results_['slice_mean_diff2']))
            results['volume_means'] = np.hstack((
                    results['volume_means'],
                    results_['volume_means']))
            results['diff2_mean_vol'] = mean_img(
                [results['diff2_mean_vol'], results_['diff2_mean_vol']])
            results['slice_diff2_max_vol'] = nib.Nifti1Image(
                np.maximum(results_['slice_diff2_max_vol'].get_data(),
                           results['slice_diff2_max_vol'].get_data()),
                results['slice_diff2_max_vol'].get_affine()
                )
            results['session_length'].append(results_['session_length'])
    return results

"""
Plot Haxby masks
=================

Small script to plot the masks of the Haxby dataset.
"""
from scipy import linalg
import matplotlib.pyplot as plt

from nilearn import datasets
data = datasets.fetch_haxby()

# Build the mean image because we have no anatomic data
from nilearn import image
mean_img = image.mean_img(data.func[0])

z_slice = -24
from nilearn.image.resampling import coord_transform
affine = mean_img.get_affine()
_, _, k_slice = coord_transform(0, 0, z_slice,
                                linalg.inv(affine))
k_slice = round(k_slice)

fig = plt.figure(figsize=(4, 5.4), facecolor='k')

from nilearn.plotting import plot_anat
display = plot_anat(mean_img, display_mode='z', cut_coords=[z_slice],
                    figure=fig)
display.add_contours(data.mask_vt[0], contours=1, antialiased=False,
                     linewidths=4., levels=[0], colors=['red'])
display.add_contours(data.mask_house[0], contours=1, antialiased=False,

As we vary the smoothing FWHM, note how we decrease the amount of noise,
but also loose spatial details. In general, the best amount of smoothing
for a given analysis depends on the spatial extent of the effects that
are expected.

"""

from nilearn import datasets, plotting, image

data = datasets.fetch_development_fmri(n_subjects=1)

# Print basic information on the dataset
print('First subject functional nifti image (4D) are located at: %s' %
      data.func[0])

first_epi_file = data.func[0]

# First the compute the mean image, from the 4D series of image
mean_func = image.mean_img(first_epi_file)

# Then we smooth, with a varying amount of smoothing, from none to 20mm
# by increments of 5mm
for smoothing in range(0, 25, 5):
    smoothed_img = image.smooth_img(mean_func, smoothing)
    plotting.plot_epi(smoothed_img,
                      title="Smoothing %imm" % smoothing)


plotting.show()

                           memory_level=1)
fmri_masked = nifti_masker.fit_transform(fmri_img)

from sklearn.feature_selection import f_classif
f_values, p_values = f_classif(fmri_masked, y)
p_values = -np.log10(p_values)
p_values[p_values > 10] = 10
p_unmasked = nifti_masker.inverse_transform(p_values).get_data()

### Visualization #############################################################
import matplotlib.pyplot as plt

# Use the fmri mean image as a surrogate of anatomical data
from nilearn import image
from nilearn.plotting import plot_stat_map
mean_fmri = image.mean_img(fmri_img)

plot_stat_map(nibabel.Nifti1Image(searchlight.scores_,
                                  mean_fmri.get_affine()), mean_fmri,
              title="Searchlight", display_mode="z", cut_coords=[-16],
              colorbar=False)

### F_score results
p_ma = np.ma.array(p_unmasked, mask=np.logical_not(process_mask))
plot_stat_map(nibabel.Nifti1Image(p_ma,
                                  mean_fmri.get_affine()), mean_fmri,
              title="F-scores", display_mode="z", cut_coords=[-16],
              colorbar=False)

plt.show()

y_train=np.array(y_train.astype('double'))

y_test[y_test=='scissors']=1
y_test[y_test=='scrambledpix']=-1
y_test=np.array(y_test.astype('double'))



masker = NiftiMasker(mask_strategy='epi',standardize=True)
                        
X_train = masker.fit_transform(X_train)
X_test  = masker.transform(X_test)

mask = masker.mask_img_.get_data().astype(np.bool)
mask= _crop_mask(mask)
background_img = mean_img(data_files.func[0])

X_train, y_train, _, y_train_mean, _ = center_data(X_train, y_train, fit_intercept=True, normalize=False,copy=False)
X_test-=X_train.mean(axis=0)
X_test/=np.std(X_train,axis=0)
alpha=1
ratio=0.5
k=200


solver_params = dict(tol=1e-6, max_iter=5000,prox_max_iter=100)

init=None
w,obj,init=tvksp_solver(X_train,y_train,alpha,ratio,k,mask=mask,init=init,loss="logistic",verbose=1,**solver_params)
coef=w[:-1]
intercept=w[-1]    

hrf_model = 'spm + derivative'  # The hemodunamic response finction is the SPM canonical one

#########################################################################
# Resample the images.
#
# This is achieved by the concat_imgs function of Nilearn.
from nilearn.image import concat_imgs, resample_img, mean_img
fmri_img = [concat_imgs(subject_data.func1, auto_resample=True),
            concat_imgs(subject_data.func2, auto_resample=True)]
affine, shape = fmri_img[0].affine, fmri_img[0].shape
print('Resampling the second image (this takes time)...')
fmri_img[1] = resample_img(fmri_img[1], affine, shape[:3])

#########################################################################
# Create mean image for display
mean_image = mean_img(fmri_img)

#########################################################################
# Make design matrices
import numpy as np
import pandas as pd
from nistats.design_matrix import make_first_level_design_matrix
design_matrices = []

#########################################################################
# loop over the two sessions
for idx, img in enumerate(fmri_img, start=1):
    # Build experimental paradigm
    n_scans = img.shape[-1]
    events = pd.read_table(subject_data['events{}'.format(idx)])
    # Define the sampling times for the design matrix

import numpy as np
import pandas as pd

import nibabel as nib
from nilearn.plotting import plot_stat_map, show
from nilearn.image import mean_img

from nistats.design_matrix import make_design_matrix
from nistats.glm import FirstLevelGLM
from nistats.datasets import fetch_spm_auditory

# fetch spm auditory data
subject_data = fetch_spm_auditory()
fmri_img = nib.concat_images(subject_data.func)
# compute bg unto which activation will be projected
mean_img = mean_img(fmri_img)

# construct experimental paradigm
tr = 7.
n_scans = 96
epoch_duration = 6 * tr  # duration in seconds
conditions = ['rest', 'active'] * 8
n_blocks = len(conditions)
duration = epoch_duration * np.ones(n_blocks)
onset = np.linspace(0, (n_blocks - 1) * epoch_duration, n_blocks)
paradigm = pd.DataFrame(
    {'onset': onset, 'duration': duration, 'name': conditions})

# construct design matrix
frame_times = np.linspace(0, (n_scans - 1) * tr, n_scans)
drift_model = 'Cosine'

                           memory_level=1)
fmri_masked = nifti_masker.fit_transform(fmri_img)

from sklearn.feature_selection import f_classif
f_values, p_values = f_classif(fmri_masked, y)
p_values = -np.log10(p_values)
p_values[np.isnan(p_values)] = 0
p_values[p_values > 10] = 10
p_unmasked = nifti_masker.inverse_transform(p_values).get_data()

### Visualization #############################################################
import matplotlib.pyplot as plt

# Use the fmri mean image as a surrogate of anatomical data
from nilearn import image
mean_fmri = image.mean_img(fmri_img).get_data()

### Searchlight results
plt.figure(1)
# searchlight.scores_ contains per voxel cross validation scores
s_scores = np.ma.array(searchlight.scores_, mask=np.logical_not(process_mask))
plt.imshow(np.rot90(mean_fmri[..., picked_slice]), interpolation='nearest',
          cmap=plt.cm.gray)
plt.imshow(np.rot90(s_scores[..., picked_slice]), interpolation='nearest',
          cmap=plt.cm.hot, vmax=1)
plt.axis('off')
plt.title('Searchlight')

### F_score results
plt.figure(2)
p_ma = np.ma.array(p_unmasked, mask=np.logical_not(process_mask))

# Normalize estimated components, for thresholding to make sense
components_masked -= components_masked.mean(axis=0)
components_masked /= components_masked.std(axis=0)
# Threshold
components_masked[components_masked < 0.8] = 0

# Now invert the masking operation, going back to a full 3D
# representation
component_img = masker.inverse_transform(components_masked)
components = component_img.get_data()

# Using a masked array is important to have transparency in the figures
components = np.ma.masked_equal(components, 0, copy=False)

### Visualize the results #####################################################
# Show some interesting components

# Use the mean as a background
import nibabel
import pylab as plt
from nilearn import image
from nilearn.plotting import plot_stat_map

mean_img = image.mean_img(dataset.func[0])

plot_stat_map(nibabel.Nifti1Image(component_img.get_data()[:, :, :, 5], component_img.get_affine()), mean_img)

plot_stat_map(nibabel.Nifti1Image(component_img.get_data()[:, :, :, 12], component_img.get_affine()), mean_img)

plt.show()

###############################################################################
# From already masked data
from nilearn.input_data import NiftiMasker
import nilearn.image as image
from nilearn.plotting.img_plotting import plot_roi

# Load Miyawaki dataset
miyawaki_dataset = datasets.fetch_miyawaki2008()

# print basic information on the dataset
print('First functional nifti image (4D) is located at: %s' %
      miyawaki_dataset.func[0])  # 4D data

miyawaki_filename = miyawaki_dataset.func[0]
miyawaki_mean_img = image.mean_img(miyawaki_filename)

# This time, we can use the NiftiMasker without changing the default mask
# strategy, as the data has already been masked, and thus lies on a
# homogeneous background

masker = NiftiMasker()
masker.fit(miyawaki_filename)

plot_roi(masker.mask_img_, miyawaki_mean_img,
         title="Mask from already masked data")


###############################################################################
# From raw EPI data