def setup_axes(fig, imx_c, imy_c, h):

    grid = axes_grid.Grid(fig, "111",
                          nrows_ncols=(4, 3), #ngrids=11,
                          direction='row', axes_pad=0.02,
                          add_all=True,
                          share_all=False, #False,
                          share_x=False, share_y=False,
                          label_mode='L',
                          axes_class=(pywcsgrid2.Axes, {"header":h}))
    grid.set_aspect(True)



    # colorbar axes
    from mpl_toolkits.axes_grid.inset_locator import inset_axes
    axins = inset_axes(grid[0],
                       width="5%",
                       height="50%",
                       loc=4,
                       )

    axins.axis[:].toggle(all=False)
    axins.axis["left"].toggle(all=True)
    axins.axis["left"].label.set_size(10)

    return grid, axins

def plot_obs_pred_sad(datasets, data_dir='./data/', radius=2):
    """Multiple obs-predicted plotter"""
    fig = plt.figure()
    num_datasets = len(datasets)
    rows = (3 if num_datasets in (5, 6) else 2)
    for i, dataset in enumerate(datasets):
        obs_pred_data = import_obs_pred_data(data_dir + dataset + '_obs_pred.csv') 
        site = ((obs_pred_data["site"]))
        obs = ((obs_pred_data["obs"]))
        pred = ((obs_pred_data["pred"]))
        
        axis_min = 0.5 * min(obs)
        axis_max = 2 * max(obs)
        ax = fig.add_subplot(rows,2,i+1)
        macroecotools.plot_color_by_pt_dens(pred, obs, radius, loglog=1, 
                                            plot_obj=plt.subplot(rows,2,i+1))      
        plt.plot([axis_min, axis_max],[axis_min, axis_max], 'k-')
        plt.xlim(axis_min, axis_max)
        plt.ylim(axis_min, axis_max)
        plt.subplots_adjust(left=0.2, bottom=0.12, right=0.8, top=0.92, 
                            wspace=0.29, hspace=0.21)  
        
        # Create inset for histogram of site level r^2 values
        axins = inset_axes(ax, width="30%", height="30%", loc=4)
        hist_mete_r2(site, np.log10(obs), np.log10(pred))
        plt.setp(axins, xticks=[], yticks=[])
        
    plt.savefig('obs_pred_plots.png', dpi=400, bbox_inches = 'tight', pad_inches=0)

def plot_eigenvalue_gaps(
        gaps, 
        xlabel=r"$k$", 
        ylabel=r"$\big|\lambda_{k+1}-\lambda_k\big|$", 
        #ylabel=r"$\lambda_{k}$", 
        output="eigen_gap.pdf",
        colors=['b', 'r', 'g', 'c'],
        legend=[r'$\rho_{1}$', r'$\widehat{\rho}_{2}$'],
        marker=['o', '^', 'v']):
    fig = plt.figure(figsize=(fig_width, fig_height))
    ax = fig.add_subplot(111)
    for i, g in enumerate(gaps):
        xs = range(1, len(g)+1)
        ax.plot(xs, g, color=colors[i], linestyle='-', marker=marker[i],
                linewidth=1.5, markersize=5, label=legend[i], alpha=0.7)
    ax.set_ylabel(ylabel)
    ax.set_xlabel(xlabel)
    ax.set_xlim([1, len(xs)])
    #ax.legend(loc=(0.18,0.6), framealpha=.5, ncol=1)
    ax.legend(loc=0, framealpha=.5, ncol=1)

    #with sns.axes_style("whitegrid"):
    axins = inset_axes(ax, width="50%", height="50%", loc=7, 
                        borderpad=1)
    for i, g in enumerate(gaps):
        axins.plot(xs[2:10], g[2:10], color=colors[i], linestyle='-', 
                       marker=marker[i], linewidth=1.5, markersize=5, 
                       alpha=0.7)
    axins.set_xlim([4,8])
    axins.set_xticks([4,5,6,7,8])
    #axins.set_yticks([])
    #axins.set_ylim([0,0.006])

    fig.savefig(output, bbox_inches='tight')

def plot_gap(infile="experiments_data2/energy_synapse_gap.csv", 
             output="gap.pdf", 
             xlabel="$k$", 
             ylabel1=r"$g_k-\left( g_{k+1} - \sigma_{k+1} \right)$", 
             ylabel2=r"$J_k$"):
    df = pd.read_csv(infile, dtype=float)
    fig = plt.figure(figsize=(fig_width, fig_height))
    # plot gaps
    ax = fig.add_subplot(111)
    xs = range(1,len(df["gap"].values)+1)
    #ax.errorbar(xs, df["gap"].values, yerr=df["var"].values, color="b",
    #            linestyle='-', marker="o", markersize=5, elinewidth=.5,
    #            capthick=0.5, linewidth=1.5, barsabove=False)
    ax.plot(xs, df["gap2"].values, color="b", linestyle="-", linewidth=1.5,
                marker="o", markersize=5)
    ax.plot(xs, [0]*len(xs), linestyle='--', color='k')
    ax.set_xlim([1, len(xs)])
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel1)
    #with sns.axes_style("whitegrid"):
    axins = inset_axes(ax, width="50%", height="50%", loc=5, 
                        borderpad=1)
    axins.plot(xs[4:10], df["gap2"].values[4:10], 
                   color='b', linestyle='-', 
                   marker='o', linewidth=1.5, markersize=5, 
                   alpha=0.7)
    axins.set_xticks([5,6,7,8,9])
    axins.set_yticks([])
    fig.savefig(output, bbox_inches='tight')

def plot_obs_pred(obs, pred, radius, loglog, ax = None, inset = False, sites = None):
    """Generic function to generate an observed vs predicted figure with 1:1 line"""
    if not ax:
        fig = plt.figure(figsize = (3.5, 3.5))
        ax = plt.subplot(111)

    axis_min = 0.9 * min(list(obs[obs > 0]) + list(pred[pred > 0]))
    if loglog:
        axis_max = 3 * max(list(obs)+list(pred))
    else:
        axis_max = 1.1 * max(list(obs)+list(pred))
    macroecotools.plot_color_by_pt_dens(np.array(pred), np.array(obs), radius, loglog=loglog, plot_obj = ax)      
    plt.plot([axis_min, axis_max],[axis_min, axis_max], 'k-')
    plt.xlim(axis_min, axis_max)
    plt.ylim(axis_min, axis_max)
    ax.tick_params(axis = 'both', which = 'major', labelsize = 6)
    if loglog:
        plt.annotate(r'$R^2$ = %0.2f' %macroecotools.obs_pred_rsquare(np.log10(obs[(obs != 0) * (pred != 0)]), np.log10(pred[(obs != 0) * (pred != 0)])),
                     xy = (0.05, 0.85), xycoords = 'axes fraction', fontsize = 7)
    else:
        plt.annotate(r'$R^2$ = %0.2f' %macroecotools.obs_pred_rsquare(obs, pred),
                     xy = (0.05, 0.85), xycoords = 'axes fraction', fontsize = 7)
    if inset:
        axins = inset_axes(ax, width="30%", height="30%", loc=4)
        if loglog:
            hist_mete_r2(sites[(obs != 0) * (pred != 0)], np.log10(obs[(obs != 0) * (pred != 0)]), 
                         np.log10(pred[(obs != 0) * (pred != 0)]))
        else:
            hist_mete_r2(sites, obs, pred)
        plt.setp(axins, xticks=[], yticks=[])
    return ax

def plot_extracellular(ax, lfp, ele_pos, num_x, num_y, time_pt):
    """Plots the extracellular potentials at a given potentials"""

    lfp *= 1000.
    lfp_max = np.max(np.abs(lfp[:, time_pt]))
    levels = np.linspace(-lfp_max, lfp_max, 16)
    im2 = plt.contourf(ele_pos[:,0].reshape(num_x, num_y), 
                       ele_pos[:,1].reshape(num_x, num_y), 
                       lfp[:,time_pt].reshape(num_x,num_y), 
                       levels=levels, cmap=plt.cm.PRGn)


    # cb = plt.colorbar(im2, extend='both')
    # tick_locator = ticker.MaxNLocator(nbins=9, trim=False, prune=None)
    # #tick_locator.bin_boundaries(-lfp_max, lfp_max)
    # cb.locator = tick_locator
    # #cb.ax.yaxis.set_major_locator(ticker.AutoLocator())
    # cb.update_ticks()
    # cb.ax.set_title('$\mu$V')

    plt.title('Time='+str(time_pt/10.)+' ms')
    plt.xlabel('X ($\mu$m)')
    plt.ylabel('Y ($\mu$m)')
    plt.ylim(ymin=-2150,ymax=550)
    plt.xlim(xmin=-450,xmax=450)
    cbaxes = inset_axes(ax,
                        width="50%",  # width = 10% of parent_bbox width
                        height="2%",  # height : 50%
                        loc=1, borderpad=1.5)
    cbar = plt.colorbar(cax=cbaxes, ticks=[-lfp_max,0.,lfp_max], orientation='horizontal', format='%.2f')
    cbar.ax.set_xticklabels([round(-lfp_max,2),str('0 $\mu V$'),round(lfp_max,2)])
    return ax

def plot_potential(ax, lfp, max_val, title):
    norm = cm.colors.Normalize(vmax=max_val, vmin=-max_val, clip=False)
    im = plt.imshow(
        lfp[::-1],
        aspect='auto',
        norm=norm,
        interpolation='nearest',
        cmap=plt.cm.PRGn)

    # plt.xticks(np.arange(3000, 5000, 1000), np.arange(300, 500, 100))
    # plt.ylabel('Electrode depth ($\mu$m)')
    # plt.xlabel('Time (ms)')
    plt.title(title, fontweight="bold", fontsize=12)
    # plt.xlim(xmin=2500, xmax=4500)
    # plt.colorbar(extend='both')
    plt.gca().set_yticks(np.arange(28))
    plt.gca().set_yticklabels(np.arange(1, 29))
    for label in plt.gca().yaxis.get_ticklabels()[::2]:
        label.set_visible(False)
    cbaxes = inset_axes(ax,
                        width="40%",  # width = 10% of parent_bbox width
                        height="3%",  # height : 50%
                        loc=1, borderpad=1)
    cbar = plt.colorbar(
        cax=cbaxes,
        ticks=[-max_val,
               0.,
               max_val],
        orientation='horizontal',
        format='%.2f')
    cbar.ax.set_xticklabels(
        [round(-max_val, 2), str('0 $\mu V$'), round(max_val, 2)])
    return ax, im

def plot_elbow_kernel(
        values, 
        xlabel=r"$k$", 
        #ylabel=r"$\textnormal{Tr}(Y^\top G \, Y)$", 
        ylabel=r"$\log Q_{k+1} - \log Q_{k}$", 
        output="elbow.pdf",
        colors=['b', 'r', 'g'],
        legend=[r'$\rho_{1}$', r'$\rho_{0.5}$', r'$\rho_{0.25}$'],
        marker=['o', '^', 'v']):
    fig = plt.figure(figsize=(fig_width, fig_height))
    ax = fig.add_subplot(111)
    for i, g in enumerate(values):
        xs = range(1, len(g)+1)
        ax.plot(xs, g, color=colors[i], linestyle='-', marker=marker[i],
                linewidth=1.5, markersize=5, label=legend[i], alpha=0.7)
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    ax.set_xlim([1, 12])
    ax.legend(loc=2, framealpha=.5, ncol=1)
    
    axins = inset_axes(ax, width="60%", height="60%", loc=1)
    for i, g in enumerate(values):
        axins.plot(xs[4:], g[4:], 
                color=colors[i], linestyle='-', marker=marker[i],
                linewidth=1.5, markersize=5, alpha=0.7)
    #axins.set_xlim([5,12])
    #axins.set_ylim([19.9,20])
    #axins.set_xticks([])
    #axins.set_yticks([])
    fig.savefig(output, bbox_inches='tight')

def plot_eigenvalue_gaps(
        gaps, 
        xlabel=r"$k$", 
        ylabel=r"$\lambda_{k+1}-\lambda_k$", 
        output="eigen_gap.pdf",
        colors=['b', 'r', 'g', 'c'],
        legend=[r'$\rho_{1}$', r'$\rho_{0.5}$', r'$\rho_{0.25}$'],
        marker=['o', '^', 'v']):
    fig = plt.figure(figsize=(fig_width, fig_height))
    ax = fig.add_subplot(111)
    for i, g in enumerate(gaps):
        xs = range(1, len(g)+1)
        ax.plot(xs, g, color=colors[i], linestyle='-', marker=marker[i],
                linewidth=1.5, markersize=5, label=legend[i], alpha=0.7)
    ax.set_ylabel(r'$\lambda_{k+1} - \lambda_k$')
    ax.set_xlabel(r'$k$')
    ax.set_xlim([1, len(xs)])
    ax.legend(loc=2, framealpha=.5, ncol=1)

    axins = inset_axes(ax, width="60%", height="60%", loc=1, 
                        borderpad=0.5)
    for i, g in enumerate(gaps):
        axins.plot(xs, g, color=colors[i], linestyle='-', marker=marker[i],
                linewidth=1.5, markersize=5, alpha=0.7)
    axins.set_xlim([5,12])
    axins.set_ylim([0,0.006])

    fig.savefig(output, bbox_inches='tight')

def plot_spectrogram(X, param, ax, colorbar = False, title = 'Spectrogram', dB= True, freqscale = 'log', dBMax = None, scaling = 'density', **kwargs):
    # TODO: correct t axis scala
    """
    plot the spectrogram of a STFT
    """
    sR = param['sR']
    # PSD
    PSD, freq, t_i =  stft_PSD(X, param, scaling = scaling, **kwargs)
    if dB:
        Z = 10*np.log10(PSD) - 20*np.log10(2e-5)
    else:
        Z = PSD
    # tempo e frequenza per questo plot é ai bordi
    df, __ = frequency_resolution(param['N'],sR)
    tR = param['R']/sR
    t = np.hstack([t_i[0]-tR/2, t_i + tR/2])
    freq = np.hstack([ freq[0]-df/2, freq + df/2])
    X , Y = np.meshgrid(t , freq)
    
    # plotting
    ax.set_title(title, fontsize = 10)
    #cmap = brewer2mpl.get_map('RdPu', 'Sequential', 9).mpl_colormap
    colormap=['#fff7f3','#fde0dd','#fcc5c0','#fa9fb5','#f768a1','#dd3497','#ae017e','#7a0177','#49006a']
    cmap=LinearSegmentedColormap.from_list('YeOrRe',colormap)
    
    #cmap=LinearSegmentedColormap('RdPu',cmap)
    if dBMax==None:
        norm = matplotlib.colors.Normalize(vmin = 0)
    else:
        norm = matplotlib.colors.Normalize(vmin = 0, vmax=dBMax)
    # np.round(np.max(ZdB)-60 ,-1), vmax = np.round(np.max(ZdB)+5,-1), clip = False)
    spect = ax.pcolormesh(X, Y, np.transpose(Z), norm=norm, cmap = cmap)
    #legenda
    if colorbar:
        axcolorbar = inset_axes(ax,
                width="2.5%", # width = 10% of parent_bbox width
                height="100%", # height : 50%
                loc='upper left',
                bbox_to_anchor=(1.01, 0., 1, 1),
                bbox_transform=ax.transAxes,
                borderpad=0,
                )
        axcolorbar.tick_params(axis='both', which='both', labelsize=8)
        ax.figure.colorbar(spect, cax = axcolorbar)
    #
    if freqscale =='log':
        ax.set_yscale('log')
    else:
        ax.set_yscale('linear')
    ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
    ax.grid(which= 'both' ,ls="-", linewidth=0.15, color='#aaaaaa', alpha=0.3)
    ax.set_xlim(t.min(),t.max())
    ax.set_ylim(freq.min(),freq.max())
    if not colorbar:
        return(spect)

def plot_divergence_boxplot_as_inset_axes(divergence, cdf, curr_axes, orig, samplers):
    box_plot_loc = 4 if cdf else 1
    bbox_anchor = (-.02, .12, 1, 1) if cdf else (-.02, -.09, 1, 1)
    in_axes = ins_loc.inset_axes(curr_axes, width="40%", # width = 40% of parent_bbox
        height="40%", # height : 40% of parent box
        loc=box_plot_loc, 
        bbox_to_anchor=bbox_anchor, 
        bbox_transform=curr_axes.transAxes, 
        borderpad=0)
    plot_divergence_boxplot(samplers, orig, divergence, ax=in_axes, notch=1, grid=False)
    return in_axes

def plot_basics(data, data_inst, fig, units):
    from powerlaw import plot_pdf, Fit, pdf
    annotate_coord = (-.4, .95)
    ax1 = fig.add_subplot(n_graphs,n_data,data_inst)
    x, y = pdf(data, linear_bins=True)
    ind = y>0
    y = y[ind]
    x = x[:-1]
    x = x[ind]
    ax1.scatter(x, y, color='r', s=.5)
    plot_pdf(data[data>0], ax=ax1, color='b', linewidth=2)
    from pylab import setp
    setp( ax1.get_xticklabels(), visible=False)

    if data_inst==1:
        ax1.annotate("A", annotate_coord, xycoords="axes fraction", fontproperties=panel_label_font)

    
    from mpl_toolkits.axes_grid.inset_locator import inset_axes
    ax1in = inset_axes(ax1, width = "30%", height = "30%", loc=3)
    ax1in.hist(data, normed=True, color='b')
    ax1in.set_xticks([])
    ax1in.set_yticks([])

    
    ax2 = fig.add_subplot(n_graphs,n_data,n_data+data_inst, sharex=ax1)
    plot_pdf(data, ax=ax2, color='b', linewidth=2)
    fit = Fit(data, xmin=1, discrete=True)
    fit.power_law.plot_pdf(ax=ax2, linestyle=':', color='g')
    p = fit.power_law.pdf()

    ax2.set_xlim(ax1.get_xlim())
    
    fit = Fit(data, discrete=True)
    fit.power_law.plot_pdf(ax=ax2, linestyle='--', color='g')
    from pylab import setp
    setp( ax2.get_xticklabels(), visible=False)

    if data_inst==1:
       ax2.annotate("B", annotate_coord, xycoords="axes fraction", fontproperties=panel_label_font)        
       ax2.set_ylabel(u"p(X)")# (10^n)")
        
    ax3 = fig.add_subplot(n_graphs,n_data,n_data*2+data_inst)#, sharex=ax1)#, sharey=ax2)
    fit.power_law.plot_pdf(ax=ax3, linestyle='--', color='g')
    fit.exponential.plot_pdf(ax=ax3, linestyle='--', color='r')
    fit.plot_pdf(ax=ax3, color='b', linewidth=2)
    
    ax3.set_ylim(ax2.get_ylim())
    ax3.set_xlim(ax1.get_xlim())
    
    if data_inst==1:
        ax3.annotate("C", annotate_coord, xycoords="axes fraction", fontproperties=panel_label_font)

    ax3.set_xlabel(units)

def plot_emp_vs_sim(study_id, data_dir = './out_files/', feas_type = 'partition', ax = None, inset = True, legend = False):
    """Plot of empirical and simulated mean-variance relationships for a given data set
    
    to help visually illustrate our results.
    Includes scatter plot of empirical data and its fitted line, scatter plot and fitted line for one
    set of simulated s_ij^2, 95 quantiles of s_ij^2 for each s_i^2 value, and the distribution of b in an inset.
    
    Input: 
    study_id - ID of the data set of interest, in the form listed in Appendix A. 
    """
    if not ax:
        fig = plt.figure(figsize = (3.5, 3.5))
        ax = plt.subplot(111)
    var_dat = get_var_sample_file(data_dir + 'taylor_QN_var_predicted_' + feas_type + '_full.txt')
    var_study = var_dat[var_dat['study'] == study_id]
    sim_var = [var_study[x][5] for x in xrange(len(var_study))] # take the first simulated sequence
    
    b_emp, inter_emp, r, p, std_err = stats.linregress(np.log(var_study['mean']), np.log(var_study['var']))
    b_list = []
    for k in xrange(len(var_study[0]) - 5):
        study_k = [var_study[x][k + 5] for x in xrange(len(var_study))]
        mean_k = [var_study['mean'][p] for p in xrange(len(var_study)) if study_k[p] != 0]
        study_k = [study_k[p] for p in xrange(len(study_k)) if study_k[p] != 0]
        b_k, inter, r, p, std_err = stats.linregress(np.log(mean_k), np.log(study_k))
        if k == 0: b_0, inter_0 = b_k, inter
        b_list.append(b_k)
   
    ax.set_xscale('log')
    ax.set_yscale('log')
    plt.scatter(var_study['mean'], var_study['var'], s = 8, c = 'black', edgecolors='none')
    emp, = plt.plot(var_study['mean'], np.exp(inter_emp) * var_study['mean'] ** b_emp, '-', c = 'black', linewidth=1.5)
    if feas_type == 'partition': plot_col = '#228B22'
    else: plot_col = '#CD69C9'
    plt.scatter(var_study['mean'], sim_var, s = 8, c = plot_col, edgecolors='none')
    sim, = plt.plot(var_study['mean'], np.exp(inter_0) * var_study['mean'] ** b_0, '-', linewidth=1.5, c = plot_col)
    ax.tick_params(axis = 'both', which = 'major', labelsize = 9)
    ax.set_xlabel('Mean', labelpad = 4, size = 10)
    ax.set_ylabel('Variance', labelpad = 4, size = 10)
    if legend:
        plt.legend([emp, sim], ['Empirical', (feas_type.title()) + 's'], loc = 4, prop = {'size': 8}) 
    if inset:
        axins = inset_axes(ax, width="30%", height="30%", loc=2)
        cov_factor = 0.2
        xs = np.linspace(0.9 * min(b_list + [b_emp]), 1.1 * max(b_list + [b_emp]), 200)
        dens_b = comp_dens(b_list, cov_factor)
        b_dens, = plt.plot(xs, dens_b(xs), c = plot_col, linewidth=1.5)
        ymax = 1.1 * max(dens_b(xs))
        plt.plot((b_emp, b_emp), (0, ymax), 'k-', linewidth = 1.5)
        plt.tick_params(axis = 'y', which = 'major', left = 'off', right = 'off', labelleft = 'off')
        plt.tick_params(axis = 'x', which = 'major', top = 'off', bottom = 'off', labelbottom = 'off')
    return ax

def plot_morp_ele(ax1, src_pos, ele_pos, pot, time_pt):
    """Plots the morphology midpoints and the electrode positions"""
    ax = plt.subplot(121, aspect='equal')
    plt.scatter(src_pos[:, 0], src_pos[:, 1],
                marker='.', alpha=0.7, color='k', s=0.6)
    plt.scatter(ele_pos[:, 0], ele_pos[:, 1],
                marker='x', alpha=0.8, color='r', s=0.9)
    # for tx in range(len(ele_pos[:,0])):
    #    plt.text(ele_pos[tx, 0], ele_pos[tx, 1], str(tx))
    ele_1 = 152
    ele_2 = 148
    plt.scatter(ele_pos[ele_1, 0], ele_pos[ele_1, 1],
                marker='s', color='r', s=14.)
    plt.scatter(ele_pos[ele_2, 0], ele_pos[ele_2, 1],
                marker='s', color='b', s=14.)
    plt.xlabel('X ($\mu$m)')
    plt.ylabel('Y ($\mu$m)')
    plt.title('Morphology and electrodes')
    plt.ylim(ymin=-2150, ymax=550)
    plt.xlim(xmin=-450, xmax=450)

    cbaxes = inset_axes(ax,
                        width="50%",  # width = 10% of parent_bbox width
                        height="17%",  # height : 50%
                        loc=4, borderpad=2.2)

    plt.plot(np.arange(6000), pot[ele_1, :], color='r', linewidth=0.5)
    plt.plot(np.arange(6000), pot[ele_2, :], color='b', linewidth=0.5)

    dummy_line = np.arange(-0.5, 0.5, 0.1)
    plt.plot(np.zeros_like(dummy_line) + time_pt,
             dummy_line, color='black', linewidth=1)

    # ax=plt.gca()
    # ax.arrow(time_pt, -0.1, 0., 0.075, head_width=0.05,
    #          head_length=0.05, width=0.1,
    #          length_includes_head=True, fc='k', ec='k')
    plt.xlim((2750, 3500))  # 4250))
    # plt.xticks(np.arange(3000, 5000, 1000), np.arange(300, 500, 100))
    plt.xticks(np.arange(2750, 3750, 250), np.arange(275, 375, 25))
    plt.ylim((-0.2, 0.12))
    plt.yticks(np.arange(-0.2, 0.1, 0.1), np.arange(-0.2, 0.1, 0.1))
    ax = plt.gca()
    ax.get_yaxis().tick_right()  # set_tick_params(direction='in')

    # inset_plt = plt.plot(cax=cbaxes,
    # cbar = plt.colorbar(cax=cbaxes, ticks=[-lfp_max,0.,lfp_max], orientation='horizontal', format='%.2f')
    # cbar.ax.set_xticklabels([round(-lfp_max,2),str('0 $\mu V$'),round(lfp_max,2)])

    return ax1

def doDistribution(players):
	Ytabu = np.load("ipad_tabu_best.npy")
	YRandom = np.load("humain_random.npy")
	Yhumain = [i["score"] for i in players]

	ntabu, binstabu, patches = plt.hist(Ytabu,bins=np.arange(20,80,0.5), cumulative=True, histtype="step",normed=True)
	nhumain, binshumain, patches = plt.hist(Yhumain,bins=np.arange(20,80,0.5), cumulative=True, histtype="step",normed=True)
	nrandom, binsrandom, patches = plt.hist(YRandom,bins=np.arange(20,200,1), cumulative=True, histtype="step",normed=True)

	plt.clf()

	fig = plt.figure(figsize=(10,5))


	plt.hist(Ytabu,bins=5,normed=True, alpha=1.0,facecolor='#FF9600',edgecolor="white",label="Tabou")
	plt.hist(Yhumain,bins=5,normed=True, alpha=1.0,facecolor='#0074C0',edgecolor="white",label="Humain")
	plt.hist(YRandom,bins=20,normed=True, alpha=1.0,facecolor='gray',edgecolor="white",label=r"Al\'eatoire")


	ax = plt.gca()
	ax.set_xlabel("Distance",fontsize=25)
	ax.set_ylabel(r"Distribution",fontsize=25)
	leg =plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
	           ncol=34, mode="expand", borderaxespad=0.,fontsize=20)
	frame = leg.get_frame()
	frame.set_facecolor('white')
	frame.set_edgecolor('none')
	h=0.2
	# ax.set_yticks(np.arange(0.0,1.0+h,h))
	from mpl_toolkits.axes_grid.inset_locator import inset_axes
	inset_axes = inset_axes(ax, 
	                    width="50%", # width = 30% of parent_bbox
	                    height=2.0, # height : 1 inch
	                    loc=1)

	plt.plot(np.arange(20,80-0.5,0.5),ntabu,linewidth=3,color="#FF9600")
	plt.plot(np.arange(20,80-0.5,0.5),nhumain,linewidth=3,color="#0074C0",alpha=1.0)

	sigma = np.std(YRandom)
	mu = np.mean(YRandom)


	print(norm.cdf(np.mean(Yhumain),mu,sigma))
	plt.xlim([30,80])
	plt.ylim([-0.001,1.001])
	plt.yticks(np.arange(0.0,1.1,0.3))
	plt.ylabel("Cumulative",fontsize=20)

	plt.savefig("../Pres_symposium/figures/distribution_human.pdf",bbox_inches='tight')

def _plot_legend(pos, colors, axis, bads, outlines):
    """Helper function to plot color/channel legends for butterfly plots
    with spatial colors"""
    from mpl_toolkits.axes_grid.inset_locator import inset_axes
    bbox = axis.get_window_extent()  # Determine the correct size.
    ratio = bbox.width / bbox.height
    ax = inset_axes(axis, width=str(30 / ratio) + '%', height='30%', loc=2)
    pos_x, pos_y = _prepare_topomap(pos, ax)
    ax.scatter(pos_x, pos_y, color=colors, s=25, marker='.', zorder=1)
    for idx in bads:
        ax.scatter(pos_x[idx], pos_y[idx], s=5, marker='.', color='w',
                   zorder=1)

    if isinstance(outlines, dict):
        _draw_outlines(ax, outlines)

def _plot_legend(pos, colors, axis, bads, outlines):
    """Helper function to plot color/channel legends for butterfly plots
    with spatial colors"""
    from mpl_toolkits.axes_grid.inset_locator import inset_axes
    bbox = axis.get_window_extent()  # Determine the correct size.
    ratio = bbox.width / bbox.height
    ax = inset_axes(axis, width=str(30 / ratio) + '%', height='30%', loc=2)
    pos_x, pos_y = _prepare_topomap(pos, ax)
    ax.scatter(pos_x, pos_y, color=colors, s=25, marker='.', zorder=0)
    for idx in bads:
        ax.scatter(pos_x[idx], pos_y[idx], s=5, marker='.', color='w',
                   zorder=1)

    if isinstance(outlines, dict):
        outlines_ = dict([(k, v) for k, v in outlines.items() if k not in
                          ['patch', 'autoshrink']])
        for k, (x, y) in outlines_.items():
            if 'mask' in k:
                continue
            ax.plot(x, y, color='k', linewidth=1)

def plot_csd(ax, lfp, max_val, title, start_t, end_t):
    norm = cm.colors.Normalize(vmax=max_val, vmin=-max_val, clip=False)
    im = plt.imshow(lfp[::-1], aspect='auto', norm=norm, interpolation='nearest', cmap=plt.cm.bwr_r)
    #plt.xlim((2750, 4250))
    plt.xticks(np.arange(3000-start_t, 5000-start_t, 1000), np.arange(300, 500, 100))
    #plt.ylabel('Electrode depth ($\mu$m)')
    plt.xlabel('Time (ms)')
    plt.gca().yaxis.set_major_locator(plt.NullLocator())
    plt.title(title)#, fontweight="bold", fontsize=12)
    #plt.yticks(np.arange(28))
    # plt.gca().set_yticks(np.arange(28))
    # plt.gca().set_yticklabels(np.arange(1,29))
    for label in plt.gca().yaxis.get_ticklabels()[::2]:
        label.set_visible(False)
    #plt.xlim(xmin=2500, xmax=4500)
    #plt.colorbar(extend='both')
    cbaxes = inset_axes(ax,
                        width="40%",  # width = 10% of parent_bbox width
                        height="3%",  # height : 50%
                        loc=1, borderpad=1 )
    cbar = plt.colorbar(cax=cbaxes, ticks=[-max_val,0.,max_val], orientation='horizontal', format='%.2f')
    cbar.ax.set_xticklabels(['source',str('0'),'sink'])
    return ax, im

	def plot(self, profile, statsResults):
		if len(statsResults.activeData) <= 0:
			self.emptyAxis()			
			return
			
		axesColour = str(self.preferences['Axes colour'].name())
		
		# *** Get data to plot 
		if self.fieldToPlot == 'p-values':
			data = statsResults.getColumn('pValues')
			xLabel = 'p-value'
		elif self.fieldToPlot == 'p-values (corrected)':
			data = statsResults.getColumn('pValuesCorrected')
			xLabel = 'p-value (corrected)'
		
		# *** Set size of figure
		self.fig.clear()
		self.fig.set_size_inches(self.figWidth, self.figHeight) 
		heightBottomLabels = 0.4	# inches
		widthSideLabel = 0.5			# inches 
		padding = 0.2						 # inches
		axesHist = self.fig.add_axes([widthSideLabel/self.figWidth,heightBottomLabels/self.figHeight,\
																		1.0-(widthSideLabel+padding)/self.figWidth,\
																		1.0-(heightBottomLabels+padding)/self.figHeight])
		
		# *** Histogram plot 
		bins = [0.0]
		binWidth = self.binWidth
		binEnd = binWidth
		while binEnd <= 1.0:
			bins.append(binEnd)
			binEnd += binWidth
			

		n, bins, patch = axesHist.hist(data, bins=bins, log=self.yAxisLogScale,color=(0.5,0.5,0.5))	
		axesHist.set_xlabel(xLabel)
		axesHist.set_ylabel('Number of features')
			

		# *** Prettify plot	 
		for a in axesHist.yaxis.majorTicks:
			a.tick1On=True
			a.tick2On=False
				
		for a in axesHist.xaxis.majorTicks:
			a.tick1On=True
			a.tick2On=False
			
		for line in axesHist.yaxis.get_ticklines(): 
			line.set_color(axesColour)
				
		for line in axesHist.xaxis.get_ticklines(): 
			line.set_color(axesColour)
			
		for loc, spine in axesHist.spines.iteritems():
			if loc in ['right','top']:
				spine.set_color('none') 
			else:
				spine.set_color(axesColour)
			
		# *** Plot inset		
		if self.bShowInset:
			bins = [0.0]
			binWidth = self.insetBinWidth
			binEnd = binWidth
			while binEnd <= self.xLimit:
				bins.append(binEnd)
				binEnd += binWidth
				
			widthStr = str(self.insetWidth) + '%'
			heightStr = str(self.insetHeight) + '%'
			axins = inset_axes(axesHist, width=widthStr, height=heightStr, loc=1)
			filteredData = [d for d in data if d <= self.xLimit]
			if filteredData:
				n, bins, patch = axins.hist(filteredData, bins=bins, log=self.insetLogScale, color=(0.5,0.5,0.5))	
			axins.set_xlim(0, self.xLimit)
			
			# *** Prettify inset	 
			for a in axins.yaxis.majorTicks:
				a.tick1On=True
				a.tick2On=False
					
			for a in axins.xaxis.majorTicks:
				a.tick1On=True
				a.tick2On=False
				
			for line in axins.yaxis.get_ticklines(): 
				line.set_color(axesColour)
				
			for line in axins.xaxis.get_ticklines(): 
				line.set_color(axesColour)
				
			for loc, spine in axins.spines.iteritems():
				if loc in ['right','top']:
					spine.set_color('none') 
				else:
						spine.set_color(axesColour)
				
		self.updateGeometry()			 
		self.draw()

#.........这里部分代码省略.........
        rad = eval(d)
        RADs.append(rad)
    DATA.close()
    
    x = []
    y = []
    for rad in RADs:
        rank = 1
        for j in rad:
            x.append(rank)
            y.append(np.log(j))
            rank += 1

    plt.hexbin(x,y,bins='log',gridsize=100,mincnt=1,cmap=cm.Reds) # bins can be 'log' or None
    plt.tick_params(axis='both', which='major', labelsize=fs-3)
    ranks = range(1,301+1)
    plt.scatter(ranks,np.log(obsRAD), color='0.3', marker = '.')
    plt.text(225, 8, '0.28',fontsize=fs+2)
    plt.xlabel("Rank in abundance",fontsize=fs)
    plt.ylabel("ln(abundance)",fontsize=fs)
    #plt.ylim(-0.1,19.5)
    plt.xlim(0.01,310)
    
    
    ax =fig.add_subplot(2,2,2)
    
    SSADs_freqs = get_SSADs(dataset)
    SSADs = hist_to_rank(SSADs_freqs)
    
    Pvals = []
    for i, obsSSAD in enumerate(SSADs):
        
        Q = sum(obsSSAD)
        N = len(obsSSAD)
        
        partitions = []
        PATH = '/home/kenlocey/combined1/' + str(Q) + '-' + str(N) + '.txt'
        if path.exists(PATH) and path.isfile(PATH) and access(PATH, R_OK):
            DATA = open(PATH,'r')
            DATA.readline() # first line has no abundance info, so skip it
            for d in DATA:
                partitions.append(eval(d))
        else: continue
        
        if len(partitions) == 0:
            print Q,N,'something is wrong'
            sys.exit()  
              
        predSSAD = parts.central_tendency(partitions)
        print sum(obsSSAD),len(obsSSAD),' ',sum(predSSAD),len(predSSAD),
        
        pval = stats.ks_2samp(obsSSAD, predSSAD)[1]
        Pvals.append(pval)
        print pval
        
        if pval > 0.05 and Q > 1000 and switch == 0:
            
            SSADs = []
            for part in partitions:
                ssad = rank_to_PrestonPlot(part) 
                ssad2 = lose_trailing_zeros(ssad)
                SSADs.append(ssad2)
            
            #print SSADs[0]
            #sys.exit()
            
            x = []
            y = []
            for ssad in SSADs:
                ab_class = 0
                for j in ssad:
                    x.append(ab_class)
                    y.append(j)
                    ab_class += 1
            
            plt.hexbin(x,y,bins=None,gridsize=25,mincnt=1,cmap=cm.Reds) 
            obsSSAD2 = rank_to_PrestonPlot(obsSSAD)
            obsSSAD3 = lose_trailing_zeros(obsSSAD2)
            bins = range(0,len(obsSSAD2))
            #plt.plot(obsSSAD, color='0.3', lw=2)
            plt.scatter(bins, obsSSAD3, color='0.3', marker = 'o')
            #plt.xlim(-0.5,max(bins))
            #plt.ylim(-0.5,max(max(predSSAD),max(obsSSAD))+10)
            plt.xlabel("Abundance class",fontsize=fs)
            plt.ylabel("frequency",fontsize=fs)
            plt.tick_params(axis='both', which='major', labelsize=fs-3)
            
            switch += 1
            
        if len(Pvals) > 50 and switch == 1: break
    # Create inset for histogram of site level r^2 values
    axins = inset_axes(ax, width="30%", height="30%", loc=1)
    D = get_kdens(Pvals)
    plt.plot(D[0],D[1],color='0.2',lw=2.0)
    plt.setp(axins, xticks=[0.0, 0.5, 1.0])
    plt.tick_params(axis='both', which='major', labelsize=fs-4)
    plt.xlim(0.0,1.0)
    #plt.ylim(0.0, max(D[1])+0.003)
    plt.subplots_adjust(wspace=0.3)
    plt.savefig('/home/kenlocey/partitions/SAD-SSAD.png', dpi=400, bbox_inches = 'tight', pad_inches=0.03)