def stat_function_tst(conn, prefix='', OUTPUT_PATH=None, threshold=0.05):
    fc = conn.hurst

    tst = Parallel(n_jobs=3, verbose=5)(delayed(ttest_group)(group, threshold, fc)
                                    for group in groups)
    
    if OUTPUT_PATH is None:
        font = {'family' : 'normal',
            'size'   : 20}
        changefont('font', **font)
        gr = ['v', 'av', 'avn']
        for i in range(3):
            title = prefix + '_'.join(groups[i])
            try:
                img = conn.masker.inverse_transform(tst[i])
                print title
                plot_stat_map(img, cut_coords=(3, -63, 36))
                plt.show()

            except ValueError:
                print "problem with tst " + title
        changefont.func_defaults
            
    else:
        for i in range(3):
            title = prefix + '_'.join(groups[i])
            output_file = os.path.join(OUTPUT_PATH, title)
            try:
                img = conn.masker.inverse_transform(tst[i])
                plot_stat_map(img, cut_coords=(3, -63, 36), output_file=output_file + '.pdf')
            except ValueError:
                print "problem with tst " + title

def createTFCEfMRIOverlayImages(folder,suffix,title='',vmax=8,display_mode='z',slices=range(-20,50,10),threshold=0.94999,plotToAxis=False,f=[],axes=[],colorbar=True,tight_layout=False,draw_cross=False,annotate=False):


    TFCEposImg,posImg,TFCEnegImg,negImg=getFileNamesfromFolder(folder,suffix)

    bg_img='./Templates/MNI152_.5mm_masked_edged.nii.gz'
    # threshold=0.949
    pos=image.math_img("np.multiply(img1,img2)",
                         img1=image.threshold_img(TFCEposImg,threshold=threshold),img2=posImg)
    neg=image.math_img("np.multiply(img1,img2)",
                         img1=image.threshold_img(TFCEnegImg,threshold=threshold),img2=negImg)
    fw=image.math_img("img1-img2",img1=pos,img2=neg)

    if plotToAxis:
        display=plotting.plot_stat_map(fw,display_mode=display_mode,threshold=0,
                                       cut_coords=slices,vmax=vmax,colorbar=colorbar,
                                       bg_img=bg_img,black_bg=False,title=title,dim=0,
                                       figure=f,axes=axes,draw_cross=draw_cross,
                                       annotate=annotate)
    else:
        display=plotting.plot_stat_map(fw,display_mode=display_mode,threshold=0,
        cut_coords=slices,vmax=vmax,colorbar=colorbar,bg_img=bg_img,
        black_bg=False,title=title,dim=0,annotate=annotate)

    if tight_layout:
        display.tight_layout()

    return display

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 draw_brain_map(self):
     cmap = plt.get_cmap('Accent')
     self.fig = plt.figure('brain_map')
     plot_stat_map(self.cluster_img, cut_coords=(0, 0, 0), output_file=None,
                   display_mode='ortho', colorbar=False, figure=self.fig,
                   axes=None, title=None, threshold=0.1, annotate=True,
                   draw_cross=False, black_bg='auto', symmetric_cbar="auto",
                   dim=True, vmax=None, cmap=cmap)

def qc_image_data(dataset, images, plot_dir='qc'):
    # Get ready
    masker = GreyMatterNiftiMasker(memory=Memory(cachedir='nilearn_cache')).fit()
    if op.exists(plot_dir):  # Delete old plots.
        shutil.rmtree(plot_dir)

    # Dataframe to contain summary metadata for neurovault images
    if dataset == 'neurovault':
        fetch_summary = pd.DataFrame(
            columns=('Figure #', 'col_id', 'image_id', 'name',
                     'modality', 'map_type', 'analysis_level',
                     'is_thresholded', 'not_mni', 'brain_coverage',
                     'perc_bad_voxels', 'perc_voxels_outside'))

    for ii, image in enumerate(images):
        im_path = image['absolute_path']
        if im_path is None:
            continue

        ri = ii % 4  # row i
        ci = (ii / 4) % 4  # column i
        pi = ii % 16 + 1  # plot i
        fi = ii / 16  # figure i

        if ri == 0 and ci == 0:
            fh = plt.figure(figsize=(16, 10))
            print('Plot %03d of %d' % (fi + 1, np.ceil(len(images) / 16.)))
        ax = fh.add_subplot(4, 4, pi)
        title = "%s%s" % (
            '(X) ' if image['rejected'] else '', op.basename(im_path))

        if dataset == 'neurovault':
            fetch_summary.loc[ii] = [
                'fig%03d' % (fi + 1), image.get('collection_id'),
                image.get('id'), title, image.get('modality'),
                image.get('map_type'), image.get('analysis_level'),
                image.get('is_thresholded'), image.get('not_mni'),
                image.get('brain_coverage'), image.get('perc_bad_voxels'),
                image.get('perc_voxels_outside')]

        # Images may fail to be transformed, and are of different shapes,
        # so we need to trasnform one-by-one and keep track of failures.
        img = cast_img(im_path, dtype=np.float32)
        img = clean_img(img)
        try:
            img = masker.inverse_transform(masker.transform(img))
        except Exception as e:
            print("Failed to mask/reshape image %s: %s" % (title, e))

        plot_stat_map(img, axes=ax, black_bg=True, title=title, colorbar=False)

        if (ri == 3 and ci == 3) or ii == len(images) - 1:
            out_path = op.join(plot_dir, 'fig%03d.png' % (fi + 1))
            save_and_close(out_path)

    # Save fetch_summary
    if dataset == 'neurovault':
        fetch_summary.to_csv(op.join(plot_dir, 'fetch_summary.csv'))

def plot_stat_map2(**kwargs):
    cut_coords = kwargs['cut_coords']
    row_l = kwargs['row_l']
    lines_nb = int(len(cut_coords) / row_l)
    for line in xrange(lines_nb):
        opt = dict(kwargs)
        opt.pop('row_l')
        opt['cut_coords'] = cut_coords[line * row_l: (line +1) *row_l]
        plotting.plot_stat_map(**opt)

def compute_hurst_and_stat(metric='dfa', regu='off', OUTPUT_PATH = '/volatile/hubert/beamer/test_hurst/', plot=False):
    conn = Hurst_Estimator(metric=metric, mask=dataset.mask,smoothing_fwhm=0, regu=regu, n_jobs=5)
    os.write(1,'fit\n')
    fc = conn.fit(dataset.func1)
    #conn.load_map(INPUT_PATH)
    os.write(1,'save\n')
    #stat_function_tst(conn, metric+' '+regu+' ', OUTPUT_PATH)
    conn.save(save_path=OUTPUT_PATH)
    if plot:
        os.write(1,'plot\n')
        a = Parallel(n_jobs=3, verbose=5)(delayed(classify_group)(group, fc)
                                        for group in groups)

        tst = Parallel(n_jobs=3, verbose=5)(delayed(ttest_group)(group, .05, fc)
                                        for group in groups)

        ost = Parallel(n_jobs=3, verbose=5)(delayed(ttest_onesample)(group, 0.05, fc)
                                            for group in ['v', 'av', 'avn'])

        mht = Parallel(n_jobs=3, verbose=5)(delayed(ttest_onesample_Hmean)(group, 0.05, fc)
                                            for group in ['v', 'av', 'avn'])

        mpt = Parallel(n_jobs=3, verbose=5)(delayed(mne_permutation_ttest)(group,0.05, fc, 1)
                                            for group in ['v', 'av', 'avn'])
        
        
        cot = Parallel(n_jobs=3, verbose=5)(delayed(ttest_onesample_coef)(np.reshape(coef['coef'], (coef['coef'].shape[0], coef['coef'].shape[-1])),
                                            0.05, fc)
                                            for coef in a)

        gr = ['v', 'av', 'avn']
        if regu=='off':
            OUTPUT_PATH = os.path.join(OUTPUT_PATH, metric)
        else:
            OUTPUT_PATH = os.path.join(OUTPUT_PATH, metric, regu)

        for i in range(3):
            title = '_'.join(groups[i])
            output_file = os.path.join(OUTPUT_PATH, title)
            img = conn.masker.inverse_transform(tst[i])
            plot_stat_map(img, cut_coords=(3, -63, 36), title=title, output_file=output_file + '.pdf')
            img = conn.masker.inverse_transform(cot[i])
            plot_stat_map(img, title='coef_map ' + title, output_file=output_file + 'coef_map.pdf')

            title = gr[i]
            output_file = os.path.join(OUTPUT_PATH, title)
            img = conn.masker.inverse_transform(ost[i])
            plot_stat_map(img, title='t-test H0 : H = 0.5 pvalue in -log10 scale groupe : ' + title, output_file= output_file + '.pdf')
            img = conn.masker.inverse_transform(mht[i])
            plot_stat_map(img, title='t-test H0 : H = 0.5 pvalue in -log10 scale groupe : ' + title, output_file= output_file + 'meanH.pdf')
            img = conn.masker.inverse_transform(mpt[i])
            plot_stat_map(img, title='t-test H0 : H = 0.5 pvalue in -log10 scale groupe : ' + title, output_file= output_file + 'mnepermutH.pdf')


        plt.figure()
        plt.boxplot(map(lambda x: x['accuracy'], a))
        plt.savefig(os.path.join(OUTPUT_PATH, 'boxplot.pdf'))

def montage(img, thr=0, mode='coronal', rows=5, cloumns=6, fsz=(10, 20)):
    """
    Make a montage using nilearn for the background
    The output figure will be 5 slices wide and 6
    slices deep

    :param img: nilearn image containing the data
    :param thr: threshold for the image
    :param mode: view mode. saggital, coronal, axial
    :param rows: number of rows in the figure
    :param cloumns: number of columns in the figure
    :param fsz: size of the figure
    :return fig: figure handle for saving or whatnot
    """
    # Hardwired view range
    sag_rng = [-65, 65]
    cor_rng = [-100, 65]
    axi_rng = [-71, 85]

    # Get the number of slices
    n_slices = rows * cloumns

    if mode == 'coronal':
        # Get the slice indices
        view_range = np.floor(np.linspace(cor_rng[0], cor_rng[1], n_slices))
        view_mode = 'y'
    if mode == 'axial':
        # Get the slice indices
        view_range = np.floor(np.linspace(axi_rng[0], axi_rng[1], n_slices))
        view_mode = 'z'
    if mode == 'saggital':
        # Get the slice indices
        view_range = np.floor(np.linspace(sag_rng[0], sag_rng[1], n_slices))
        view_mode = 'x'

    # Prepare the figure
    fig = plt.figure(figsize=fsz)
    gs = gridspec.GridSpec(cloumns, 1, hspace=0, wspace=0)
    # Loop through the rows of the image
    for row_id in range(cloumns):
        # Create the axis to show
        ax = fig.add_subplot(gs[row_id, 0])
        # Get the slices in the column direction
        row_range = view_range[row_id*rows:(row_id+1)*rows]
        # Display the thing
        nlp.plot_stat_map(img, cut_coords=row_range,
                          display_mode=view_mode, threshold=thr,
                          axes=ax, black_bg=True)

    return fig

def diff_computed_hurst(metric='wavelet', regu='off', INPUT_PATH = '/volatile/hubert/beamer/test_hurst/', OUTPUT_PATH=''):
    conn = Hurst_Estimator(metric=metric, mask=dataset.mask, regu=regu, n_jobs=5)
    os.write(1,'load\n')
    conn.load_map(INPUT_PATH)
    fc = conn.hurst
    os.write(1,'stat\n')

    tst = ttest_group(['av', 'v'], .05, fc)
    vmean_avmean = np.mean([fc[i] for i in dataset.group_indices['v']], axis=0) - np.mean([fc[i] for i in dataset.group_indices['av']], axis=0)
    vmean_avmean[tst == 0] = 0
    
    img = conn.masker.inverse_transform(vmean_avmean)
    plot_stat_map(img)
    plt.show()

def plot_contrast(first_level_model):
    """ Given a first model, specify, enstimate and plot the main contrasts"""
    design_matrix = first_level_model.design_matrices_[0]
    # Call the contrast specification within the function
    contrasts = make_localizer_contrasts(design_matrix)
    fig = plt.figure(figsize=(11, 3))
    # compute the per-contrast z-map
    for index, (contrast_id, contrast_val) in enumerate(contrasts.items()):
        ax = plt.subplot(1, len(contrasts), 1 + index)
        z_map = first_level_model.compute_contrast(
            contrast_val, output_type='z_score')
        plotting.plot_stat_map(
            z_map, display_mode='z', threshold=3.0, title=contrast_id, axes=ax,
            cut_coords=1)

def run(idx, reduction, alpha, mask, raw, n_components, init, func_filenames):
    output_dir = join(trace_folder, 'experiment_%i' % idx)
    try:
        os.makedirs(output_dir)
    except OSError:
        pass
    dict_fact = SpcaFmri(mask=mask,
                         smoothing_fwhm=3,
                         batch_size=40,
                         shelve=not raw,
                         n_components=n_components,
                         replacement=False,
                         dict_init=fetch_atlas_smith_2009().rsn70 if
                         init else None,
                         reduction=reduction,
                         alpha=alpha,
                         random_state=0,
                         n_epochs=2,
                         l1_ratio=0.5,
                         backend='c',
                         memory=expanduser("~/nilearn_cache"), memory_level=2,
                         verbose=5,
                         n_jobs=1,
                         trace_folder=output_dir
                         )

    print('[Example] Learning maps')
    t0 = time.time()
    dict_fact.fit(func_filenames, raw=raw)
    t1 = time.time() - t0
    print('[Example] Dumping results')
    # Decomposition estimator embeds their own masker
    masker = dict_fact.masker_
    components_img = masker.inverse_transform(dict_fact.components_)
    components_img.to_filename(join(output_dir, 'components_final.nii.gz'))
    print('[Example] Run in %.2f s' % t1)
    # Show components from both methods using 4D plotting tools
    import matplotlib.pyplot as plt
    from nilearn.plotting import plot_prob_atlas, show

    print('[Example] Displaying')
    fig, axes = plt.subplots(2, 1)
    plot_prob_atlas(components_img, view_type="filled_contours",
                    axes=axes[0])
    plot_stat_map(index_img(components_img, 0),
                  axes=axes[1],
                  colorbar=False,
                  threshold=0)
    plt.savefig(join(output_dir, 'components.pdf'))
    show()

def make_thresholded_slices(regions, colors, display_mode='z', overplot=True, binarize=True, **kwargs):
    """ Plots on axial slices numerous images
    regions: Nibabel images
    colors: List of colors (rgb tuples)
    overplot: Overlay images?
    binarize: Binarize images or plot full stat maps
    """             

    from matplotlib.colors import LinearSegmentedColormap
    from nilearn import plotting as niplt
    
    if binarize:
        for reg in regions:
             reg.get_data()[reg.get_data().nonzero()] = 1
                                   
    for i, reg in enumerate(regions):
        reg_color = LinearSegmentedColormap.from_list('reg1', [colors[i], colors[i]])
        if i == 0:
            plot = niplt.plot_stat_map(reg, draw_cross=False,  display_mode=display_mode, cmap = reg_color, alpha=0.9, colorbar=False, **kwargs)
        else:
            if overplot:
                plot.add_overlay(reg, cmap = reg_color, alpha=.72)
            else:
                plt.plot_stat_map(reg, draw_cross=False,  display_mode=display_mode, cmap = reg_color, colorbar=False, **kwargs)
    
    return plot

def make_stat_image(nifti_file,png_img_file=None):
    """Make statmap image"""
    nifti_file = str(nifti_file)
    brain = plot_stat_map(nifti_file)
    if png_img_file:    
        brain.savefig(png_img_file)
    plt.close('all')
    return brain

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 plot_stat_overlay(stat_img, contour_img, bg_img, **kwargs):
    """Plot over bg_img a stat_img and the countour."""
    import nilearn.plotting as niplot

    if bg_img is not None:
        kwargs['bg_img'] = bg_img

    display = niplot.plot_stat_map(stat_img, **kwargs)
    display.add_contours(contour_img, filled=True, alpha=0.6, levels=[0.5], colors='g')
    return display

def dump_comps(masker, compressor, components, threshold=2, fwhm=None,
               perc=None):
    from scipy.stats import zscore
    from nilearn.plotting import plot_stat_map
    from nilearn.image import smooth_img
    from scipy.stats import scoreatpercentile

    if isinstance(compressor, basestring):
        comp_name = compressor
    else:
        comp_name = compressor.__str__().split('(')[0]

    for i_c, comp in enumerate(components):
        path_mask = op.join(WRITE_DIR, '%s_%i-%i' % (comp_name,
                                                     n_comp, i_c + 1))
        nii_raw = masker.inverse_transform(comp)
        nii_raw.to_filename(path_mask + '.nii.gz')
        
        comp_z = zscore(comp)
        
        if perc is not None:
            cur_thresh = scoreatpercentile(np.abs(comp_z), per=perc)
            path_mask += '_perc%i' % perc
            print('Applying percentile %.2f (threshold: %.2f)' % (perc, cur_thresh))
        else:
            cur_thresh = threshold
            path_mask += '_thr%.2f' % cur_thresh
            print('Applying threshold: %.2f' % cur_thresh)

        nii_z = masker.inverse_transform(comp_z)
        gz_path = path_mask + '_zmap.nii.gz'
        nii_z.to_filename(gz_path)
        plot_stat_map(gz_path, bg_img='colin.nii', threshold=cur_thresh,
                      cut_coords=(0, -2, 0), draw_cross=False,
                      output_file=path_mask + 'zmap.png')
                      
        # optional: do smoothing
        if fwhm is not None:
            nii_z_fwhm = smooth_img(nii_z, fwhm=fwhm)
            plot_stat_map(nii_z_fwhm, bg_img='colin.nii', threshold=cur_thresh,
                          cut_coords=(0, -2, 0), draw_cross=False,
                          output_file=path_mask +
                          ('zmap_%imm.png' % fwhm))

def plot_components(ica_image, hemi='', out_dir=None,
                    bg_img=datasets.load_mni152_template()):
    print("Plotting %s components..." % hemi)

    # Determine threshoold and vmax for all the plots
    # get nonzero part of the image for proper thresholding of
    # r- or l- only component
    nonzero_img = ica_image.get_data()[np.nonzero(ica_image.get_data())]
    thr = stats.scoreatpercentile(np.abs(nonzero_img), 90)
    vmax = stats.scoreatpercentile(np.abs(nonzero_img), 99.99)
    for ci, ic_img in enumerate(iter_img(ica_image)):

        title = _title_from_terms(terms=ica_image.terms, ic_idx=ci, label=hemi)
        fh = plt.figure(figsize=(14, 6))
        plot_stat_map(ic_img, axes=fh.gca(), threshold=thr, vmax=vmax,
                      colorbar=True, title=title, black_bg=True, bg_img=bg_img)

        # Save images instead of displaying
        if out_dir is not None:
            save_and_close(out_path=op.join(
                out_dir, '%s_component_%i.png' % (hemi, ci)))

def dump_comps(masker, compressor, components, threshold=2):
    from scipy.stats import zscore
    from nilearn.plotting import plot_stat_map

    if isinstance(compressor, basestring):
        comp_name = compressor
    else:
        comp_name = compressor.__str__().split('(')[0]

    for i_c, comp in enumerate(components):
        path_mask = op.join(WRITE_DIR, '%s_%i-%i' % (comp_name,
                                                     n_comp, i_c + 1))
        nii_raw = masker.inverse_transform(comp)
        nii_raw.to_filename(path_mask + '.nii.gz')

        nii_z = masker.inverse_transform(zscore(comp))
        gz_path = path_mask + '_zmap.nii.gz'
        nii_z.to_filename(gz_path)
        plot_stat_map(gz_path, bg_img='colin.nii', threshold=threshold,
                      cut_coords=(0, -2, 0), draw_cross=False,
                      output_file=path_mask + 'zmap.png')

def main(stat_map="nii/tstat.nii.gz", template="~/NIdata/templates/medres_QBI_chr.nii.gz", black_bg=False, cut_coords=(-50,8,45)):
	template = path.expanduser(template)

	colors_plus = plt.cm.autumn(np.linspace(0., 1, 128))
	colors_minus = plt.cm.winter(np.linspace(0, 1, 128))
	colors = np.vstack((colors_minus, colors_plus[::-1]))

	mymap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors)

	for i in ['none', 'nearest', 'bilinear', 'bicubic', 'spline16','spline36', 'hanning', 'hamming', 'hermite', 'kaiser', 'quadric','catrom', 'gaussian', 'bessel', 'mitchell', 'sinc', 'lanczos']:
		display = plotting.plot_stat_map(stat_map, bg_img=template,threshold=2.5, vmax=40, cmap=mymap, black_bg=black_bg, cut_coords=cut_coords, annotate=True, title=i+" interpolation", draw_cross=False, interpolation=i)
		plt.savefig(i, dpi=700)

def generate_images(components_img, n_components, images_dir, glass=False):
    # Remove existing images
    if os.path.exists(images_dir):
        shutil.rmtree(images_dir)
    os.makedirs(images_dir)
    output_filenames = [osp.join(images_dir, 'IC_{}.png'.format(i))
                        for i in range(n_components)]

    for i, output_file in enumerate(output_filenames):
        plot_stat_map(nibabel.Nifti1Image(components_img.get_data()[..., i],
                                          components_img.get_affine()),
                      display_mode="z", title="IC %d" % i, cut_coords=7,
                      colorbar=False, output_file=output_file)
    if glass:
        output_filenames = [osp.join(images_dir, 'glass_IC_{}.png'.format(i))
                            for i in range(n_components)]
        for i, output_file in enumerate(output_filenames):
            plot_glass_brain(nibabel.Nifti1Image(
                components_img.get_data()[..., i],
                components_img.get_affine()),
                display_mode="ortho", title="IC %d" % i,
                             output_file=output_file)