def make_intens_all(w1, w2):
    fig = plt.figure(figsize=(6., 6.))
    gs = gridspec.GridSpec(1,1)
    gs.update(left=0.13, right=0.985, bottom = 0.13, top=0.988)
    ax = plt.subplot(gs[0])
    plt.minorticks_on()
    make_contours()
    labels = ["A", "B", "C", "D"]
    for i, field in enumerate(fields):
        os.chdir(os.path.join(data_dir, "combined_{0}".format(field)))
        image = "collapsed_w{0}_{1}.fits".format(w1, w2)
        intens = pf.getdata(image, verify=False)
        extent = calc_extent(image)
        extent = offset_extent(extent, field)
        plt.imshow(intens, cmap="bone", origin="bottom", extent=extent,
                   vmin=-20, vmax=80)
        verts = calc_verts(intens, extent)
        path = Path(verts, [Path.MOVETO, Path.LINETO, Path.LINETO, Path.LINETO,
                    Path.CLOSEPOLY,])
        patch = patches.PathPatch(path, facecolor='none', lw=2, edgecolor="r")
        ax.add_patch(patch)
        xtext, ytext = np.mean(verts[:-1], axis=0)
        plt.text(xtext-8, ytext+8, labels[i], color="r",
                fontsize=35, fontweight='bold', va='top')
        plt.hold(True)
    plt.xlim(26, -38)
    plt.ylim(-32, 32)
    plt.xlabel("X [kpc]")
    plt.ylabel("Y [kpc]")
    # plt.show()
    plt.savefig(os.path.join(plots_dir, "muse_fields.eps"), dpi=60,
                format="eps")
    plt.savefig(os.path.join(plots_dir, "muse_fields.png"), dpi=200)
    return

def plotErrorBars(dict_to_plot, x_lim, y_lim, xlabel, y_label, title, out_file, margin=[0.05, 0.05], loc=2):

    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel(y_label)

    if y_lim is None:
        y_lim = [1 * float("Inf"), -1 * float("Inf")]

    max_val_seen_y = y_lim[1] - margin[1]
    min_val_seen_y = y_lim[0] + margin[1]
    print min_val_seen_y, max_val_seen_y
    max_val_seen_x = x_lim[1] - margin[0]
    min_val_seen_x = x_lim[0] + margin[0]
    handles = []
    for k in dict_to_plot:
        means, stds, x_vals = dict_to_plot[k]

        min_val_seen_y = min(min(np.array(means) - np.array(stds)), min_val_seen_y)
        max_val_seen_y = max(max(np.array(means) + np.array(stds)), max_val_seen_y)

        min_val_seen_x = min(min(x_vals), min_val_seen_x)
        max_val_seen_x = max(max(x_vals), max_val_seen_x)

        handle = plt.errorbar(x_vals, means, yerr=stds)
        handles.append(handle)
        print max_val_seen_y
    plt.xlim([min_val_seen_x - margin[0], max_val_seen_x + margin[0]])
    plt.ylim([min_val_seen_y - margin[1], max_val_seen_y + margin[1]])
    plt.legend(handles, dict_to_plot.keys(), loc=loc)
    plt.savefig(out_file)

def scree_plot(pca_obj, fname=None): 
    '''
    Scree plot for variance & cumulative variance by component from PCA. 

    Arguments: 
        - pca_obj: a fitted sklearn PCA instance
        - fname: path to write plot to file

    Output: 
        - scree plot 
    '''   
    components = pca_obj.n_components_ 
    variance = pca.explained_variance_ratio_
    plt.figure()
    plt.plot(np.arange(1, components + 1), np.cumsum(variance), label='Cumulative Variance')
    plt.plot(np.arange(1, components + 1), variance, label='Variance')
    plt.xlim([0.8, components]); plt.ylim([0.0, 1.01])
    plt.xlabel('No. Components', labelpad=11); plt.ylabel('Variance Explained', labelpad=11)
    plt.legend(loc='best') 
    plt.tight_layout() 
    if fname is not None:
        plt.savefig(fname)
        plt.close() 
    else:
        plt.show() 
    return 

def tuning(x, y, err=None, smooth=None, ylabel=None, pal=None):
    """
    Plot a tuning curve
    """
    if smooth is not None:
        xs, ys = smoothfit(x, y, smooth)
        plt.plot(xs, ys, linewidth=4, color="black", zorder=1)
    else:
        ys = asarray([0])
    if pal is None:
        pal = sns.color_palette("husl", n_colors=len(x) + 6)
        pal = pal[2 : 2 + len(x)][::-1]
    plt.scatter(x, y, s=300, linewidth=0, color=pal, zorder=2)
    if err is not None:
        plt.errorbar(x, y, yerr=err, linestyle="None", ecolor="black", zorder=1)
    plt.xlabel("Wall distance (mm)")
    plt.ylabel(ylabel)
    plt.xlim([-2.5, 32.5])
    errTmp = err
    errTmp[isnan(err)] = 0
    rng = max([nanmax(ys), nanmax(y + errTmp)])
    plt.ylim([0 - rng * 0.1, rng + rng * 0.1])
    plt.yticks(linspace(0, rng, 3))
    plt.xticks(range(0, 40, 10))
    sns.despine()
    return rng

def draw_stat(actual_price, action):
	price_list = []
	x_list = []
	# idx = np.where(actual_price == 0)[0]
	# print idx
	# print actual_price[np.where(actual_price < 2000)]
	# idx = [0] + idx.tolist()
	# print idx
	# for i in range(len(idx)-1):
	# 	price_list.append(actual_price[idx[i]+1:idx[i+1]-1])
	# 	x_list.append(range(idx[i]+i+1, idx[i+1]+i-1))
	# for i in range(len(idx)-1):
	# 	print x_list[i]
	# 	print price_list[i]
	# 	plt.plot(x_list[i], price_list[i], 'r')
	x_list = range(1,50)
	price_list = actual_price[1:50]
	plt.plot(x_list, price_list, 'k')
	for i in range(1, 50):
		style = 'go'
		if action[i] == 1:
			style = 'ro'
		plt.plot(i, actual_price[i], style)
	plt.ylim(2140, 2144.2)
	# plt.show()
	plt.savefig("action.png")

def make_fish(zoom=False):
    plt.close(1)
    plt.figure(1, figsize=(6, 4))
    plt.plot(plot_limits['pitch'], plot_limits['rolldev'], '-g', lw=3)
    plt.plot(plot_limits['pitch'], -plot_limits['rolldev'], '-g', lw=3)
    plt.plot(pitch.midvals, roll.midvals, '.b', ms=1, alpha=0.7)

    p, r = make_ellipse()  # pitch, off nominal roll
    plt.plot(p, r, '-c', lw=2)

    gf = -0.08  # Fudge on pitch value for illustrative purposes
    plt.plot(greta['pitch'] + gf, -greta['roll'], '.r', ms=1, alpha=0.7)
    plt.plot(greta['pitch'][-1] + gf, -greta['roll'][-1], 'xr', ms=10, mew=2)

    if zoom:
        plt.xlim(46.3, 56.1)
        plt.ylim(4.1, 7.3)
    else:
        plt.ylim(-22, 22)
        plt.xlim(40, 180)
    plt.xlabel('Sun pitch angle (deg)')
    plt.ylabel('Sun off-nominal roll angle (deg)')
    plt.title('Mission off-nominal roll vs. pitch (5 minute samples)')
    plt.grid()
    plt.tight_layout()
    plt.savefig('fish{}.png'.format('_zoom' if zoom else ''))

def entries_histogram(turnstile_weather):
    '''
    Before we perform any analysis, it might be useful to take a
    look at the data we're hoping to analyze. More specifically, lets 
    examine the hourly entries in our NYC subway data and determine what
    distribution the data follows. This data is stored in a dataframe
    called turnstile_weather under the ['ENTRIESn_hourly'] column.
    
    Why don't you plot two histograms on the same axes, showing hourly
    entries when raining vs. when not raining. Here's an example on how
    to plot histograms with pandas and matplotlib:
    turnstile_weather['column_to_graph'].hist()
    
    Your histograph may look similar to the following graph:
    http://i.imgur.com/9TrkKal.png
    
    You can read a bit about using matplotlib and pandas to plot
    histograms:
    http://pandas.pydata.org/pandas-docs/stable/visualization.html#histograms
    
    You can look at the information contained within the turnstile weather data at the link below:
    https://www.dropbox.com/s/meyki2wl9xfa7yk/turnstile_data_master_with_weather.csv
    '''
    plt.figure()
    (turnstile_weather[turnstile_weather.rain==0].ENTRIESn_hourly).hist(bins=175) # your code here to plot a historgram for hourly entries when it is not raining
    (turnstile_weather[turnstile_weather.rain==1].ENTRIESn_hourly).hist(bins=175) # your code here to plot a historgram for hourly entries when it is raining
    plt.ylim(ymax = 45000, ymin = 0)
    plt.xlim(xmax = 6000, xmin = 0)
    return plt

def main( args ):

  hash = get_genes_with_features(args['file'])
  for key, featurearray in hash.iteritems():
    cluster, branch = key.split()
    length = int(featurearray[0][0])
    import matplotlib.pyplot as P
    x = [e+1 for e in range(length+1)]
    y1 = [0] * (length+1)
    y2 = [0] * (length+1)
    for feature in featurearray:
      length, pos, aa, prob = feature[0:4]
      if prob > 0.95: y1[pos] = prob
      else: y2[pos] = prob
    
    P.bar(x, y1, color='#000000', edgecolor='#000000')
    P.bar(x, y2, color='#bbbbbb', edgecolor='#bbbbbb')
    P.ylim(ymin=0, ymax=1)
    P.xlim(xmin=0, xmax=length)
    P.xlabel("position in the ungapped alignment [aa]")
    P.ylabel(r'$P (\omega > 1)$')
    P.title(cluster + " (branch " + branch + ")")

    P.axhline(y=.95, xmin=0, xmax=length, linestyle=":", color="k")
    P.savefig(cluster + "." + branch + ".png", format="png")
    P.close()

    def exec_transmissions():
        IP,IP_AP,files=parser_reduce()
        plt.figure("GRAPHE_D'EVOLUTION_DES_TRANSMISSIONS")
        ENS_TEMPS_, TRANSMISSION_ = transmissions(files)
        plt.plot(ENS_TEMPS_, TRANSMISSION_,"r.", label="Transmissions: ")

        lot = map(inet_aton, IP)
        lot.sort()
        iplist1 = map(inet_ntoa, lot)

        for i in iplist1: #ici j'affiche les annotations et vérifie si j'ai des @ip de longueur 9 ou 8 pour connaitre la taille de la fenetre du graphe
                if len(i)==9:
                    maxim_=i[-2:] #Sera utilisé pour la taille de la fenetre du graphe
                    plt.annotate('   Machine: '+ i ,horizontalalignment='left', xy=(1, float(i[-2:])), xytext=(1, float(i[-2:])-0.4),arrowprops=dict(facecolor='black', shrink=0.05),)
                else:
                    maxim_=i[-1:] #Sera utilisé pour la taille de la fenetre du graphe
                    plt.annotate('   Machine: '+ i ,horizontalalignment='left', xy=(1, float(i[7])), xytext=(1, float(i[7])-0.4),arrowprops=dict(facecolor='black', shrink=0.05),)
        for i in IP_AP: #ACCESS POINT ( cas spécial )
            if i[-2:]:
                plt.annotate('   access point: '+ i , xy=(1, i[7]), xytext=(1, float(i[7])-0.4),arrowprops=dict(facecolor='black', shrink=0.05),)

        plt.ylim(0, (float(maxim_))+1) #C'est à ça que sert le tri
        plt.xlim(1, 1.1)
        plt.legend(loc='best',prop={'size':10})
        plt.xlabel('Temps (s)')
        plt.ylabel('IP machines transmettrices')
        plt.grid(True)
        plt.title("GRAPHE_D'EVOLUTION_DES_TRANSMISSIONS")
        plt.legend(loc='best')
        plt.show()

def plot_1D_pmf(coord,title):
    ''' Plot a 1D pmf for a coordinate.'''

    x = get_data(coord)

    path = os.getcwd()
    savedir = path+"/pmfs"
    if os.path.exists(savedir) == False:
        os.mkdir(savedir)

    if coord in ["Rg","rmsd"]:
        skip = 80
    else:
        skip = 4 
    vals = np.unique(list(x))[::skip]

    n,bins = np.histogram(x,bins=vals,density=True)
    np.savetxt(savedir+"/"+coord+"_n.dat",n,delimiter=" ",fmt="%.4f")
    np.savetxt(savedir+"/"+coord+"_bins.dat",bins,delimiter=" ",fmt="%.4f")

    pmf = -np.log(n)
    pmf -= min(pmf)

    plt.figure()
    plt.plot(bins[1:]/max(bins),pmf)
    plt.xlabel(coord,fontsize="xx-large")
    plt.ylabel("F("+coord+") / kT",fontsize="xx-large")
    plt.title("F("+coord+") "+title,fontsize="xx-large")
    plt.ylim(0,6)
    plt.xlim(0,1)
    plt.savefig(savedir+"/"+coord+"_pmf.pdf")

def test_draw_residual_blh_norm():
    np.random.seed(0)
    data = np.random.randn(1000)
    blh = BinnedLH(gaussian, data)
    blh.draw_residual(args=(0., 1.), norm=True)
    plt.ylim(-4., 3.)
    plt.xlim(-4., 3.)

def visualizeEigenvalues(eVal, verboseLevel):
	real = []
	imag = []

	for z in eVal:
		rp = z.real
		im = z.imag
		if not (rp == np.inf or rp == - np.inf) \
				and not (im == np.inf or im == - np.inf):
			real.append(rp)
			imag.append(im)

	if verboseLevel>=1:
		print("length of regular real values=" + str(len(real)))
		print("length of regular imag values=" + str(len(imag)))
		print("minimal real part=" + str(min(real)), "& maximal real part=" + str(max(real)))
		print("minimal imag part=" + str(min(imag)), "& maximal imag part=" + str(max(imag)))
	if verboseLevel==2:
		print("all real values:", str(real))
		print("all imag values:", str(imag))


	# plt.scatter(real[4:],img[4:])
	plt.scatter(real, imag)
	plt.grid(True)
	plt.xlabel("realpart")
	plt.ylabel("imagpart")
	plt.xlim(-10, 10)
	plt.ylim(-10, 10)
	plt.show()

def example(show=True, save=False):

    # Settings:
    t0 = 0.
    dt = .0001
    dv = .0001
    tf = .1
    verbose = True
    update_method = 'approx'
    approx_order = 1
    tol = 1e-14
    
    # Run simulation:
    simulation = get_simulation(dv=dv, verbose=verbose, update_method=update_method, approx_order=approx_order, tol=tol)
    simulation.run(dt=dt, tf=tf, t0=t0)
    
    # Visualize:
    i1 = simulation.population_list[1]
    plt.figure(figsize=(3,3))
    plt.plot(i1.t_record, i1.firing_rate_record)
    plt.xlim([0,tf])
    plt.ylim(ymin=0)
    plt.xlabel('Time (s)')
    plt.ylabel('Firing Rate (Hz)')
    plt.tight_layout()
    if save == True: plt.savefig('./excitatory_inhibitory.png')
    if show == True: plt.show()
    
    return i1.t_record, i1.firing_rate_record

def draw(ord_l, gaps):

    axScatter = plt.subplot(3, 1, 1)

    number_samples=0
    # axScatter.scatter([i['seq'] for i in ord_l[-number_samples:]], [i['a'] for i in ord_l[-number_samples:]], s=2, color='r', label='ch1')
    axScatter.scatter([i['seq'] % 24 for i in ord_l[-number_samples:]], [i['d'] for i in ord_l[-number_samples:]], s=2, color='r', label='ch1')
    # axScatter.scatter(time_l[-number_samples:], b_l[-number_samples:], s=2, color='c', label='ch2')
    # axScatter.scatter(time_l[-number_samples:], c_l[-number_samples:], s=2, color='y', label='ch3')
    # axScatter.scatter(time_l[-number_samples:], d_l[-number_samples:], s=2, color='g', label='ch4')
    plt.ylim(-9000000, 9000000)
    plt.legend()
    axScatter.set_xlabel("Sequence Packet")
    axScatter.set_ylabel("Voltage")
    plt.title("Channels Values")


    # time_plot = plt.subplot(3, 1, 2)
    # time_plot.scatter([i['seq'] for i in ord_l[-number_samples:]], [i['delta'] for i in ord_l[-number_samples:]], s=1, color='r', label='delta')
    # time_plot.set_xlabel("Sequence Packet")
    # time_plot.set_ylabel("Delta to referencial")
    # ax2 = time_plot.twinx()
    # ax2.scatter([i['seq'] for i in ord_l[-number_samples:]], [i['ts'] for i in ord_l[-number_samples:]], s=2, color='g', label='Timestamp')
    # ax2.set_ylabel("Kernel time")
    # plt.title("Timestamp deltas")

    gaps_draw = plt.subplot(3, 1, 3)
    gaps_draw.plot([i[0] for i in gaps[-number_samples:]], [i[1] for i in gaps[-number_samples:]], color='b', marker='.', label='gaps')
    gaps_draw.set_ylim(-0.5, 1.5)

    plt.draw()
    # plt.savefig("res.png")
    plt.show()

def plot_decision_regions(X, y, classifier, test_idx=None, resolution=0.02):
    # setup marker generator and color map
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    # plot the decision surface
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                           np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    # plot class samples
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1],
                    alpha=0.8, c=cmap(idx),
                    marker=markers[idx], label=cl)

    # Highlight test samples
    if test_idx:
        X_test, y_test = X[test_idx, :], y[test_idx]
        plt.scatter(X_test[:, 0],
                    X_test[:, 1],
                    c='',
                    alpha=1.0,
                    linewidths=1,
                    marker='o',
                    s=55, label='test set')

def plt_data():
    t = [[0,1], [1,0], [1, 1], [0, 0]]
    t2 = [1, 1, 1, 0]
    X = np.array(t)
    Y = np.array(t2)

    h = .02  # step size in the mesh

    logreg = linear_model.LogisticRegression(C=1e5)

    # we create an instance of Neighbours Classifier and fit the data.
    logreg.fit(X, Y)
    # Plot the decision boundary. For that, we will assign a color to each
    # point in the mesh [x_min, m_max]x[y_min, y_max].
    x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5
    y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
    Z = logreg.predict(np.c_[xx.ravel(), yy.ravel()])

    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    plt.figure(1, figsize=(4, 3))
    plt.pcolormesh(xx, yy, Z, cmap=plt.cm.Paired)

    # Plot also the training points
    plt.scatter(X[:, 0], X[:, 1], c=Y, edgecolors='k', cmap=plt.cm.Paired)
    plt.xlabel('Sepal length')
    plt.ylabel('Sepal width')

    plt.xlim(xx.min(), xx.max())
    plt.ylim(yy.min(), yy.max())
    plt.xticks(())
    plt.yticks(())

    plt.show()

def disc_norm():
    x = np.linspace(-3,3,100)
    y = st.norm.pdf(x,0,1)
    fig, ax = plt.subplots()
    fig.canvas.draw()
    
    ax.plot(x,y)
    
    fill1_x = np.linspace(-2,-1.5,100)
    fill1_y = st.norm.pdf(fill1_x,0,1)
    fill2_x = np.linspace(-1.5,-1,100)
    fill2_y = st.norm.pdf(fill2_x,0,1)
    ax.fill_between(fill1_x,0,fill1_y,facecolor = 'blue', edgecolor = 'k',alpha = 0.75)
    ax.fill_between(fill2_x,0,fill2_y,facecolor = 'blue', edgecolor = 'k',alpha = 0.75)
    for label in ax.get_yticklabels():
        label.set_visible(False)
    for tick in ax.get_xticklines():
        tick.set_visible(False)
    for tick in ax.get_yticklines():
        tick.set_visible(False)
    
    plt.rc("font", size = 16)
    plt.xticks([-2,-1.5,-1])
    labels = [item.get_text() for item in ax.get_xticklabels()]
    labels[0] = r"$v_k$"
    labels[1] = r"$\varepsilon_k$"
    labels[2] = r"$v_{k+1}$"
    ax.set_xticklabels(labels)
    plt.ylim([0, .45])

    
    plt.savefig('discnorm.pdf')
    plt.clf()

def plot_twoscales(name, dict_array, xlabel='', ylabel='', title='', linetypes=['b','r','g','k'], labels=[], xlog=None, ylim=None):
  plt.clf()
  if len(xlabel) > 0:
  	plt.xlabel(xlabel)
  if len(ylabel) > 0:
  	plt.ylabel(ylabel)
  if len(title) > 0:
  	plt.title(title)
  if xlog:
    plt.xscale('log', basex=xlog)
  if ylim:
    plt.ylim(ylim)

  ax1 = plt.figure().add_subplot(111)  
  dicty1 = zip(*sorted(dict_array[0].iteritems()))
  dicty2 = zip(*sorted(dict_array[1].iteritems()))
  ax1.plot(dicty1[0], dicty1[1], linetypes[0], label=labels[0])
  for tl in ax1.get_yticklabels():
    tl.set_color(linetypes[0])
  ax1.set_ylabel(labels[0], color=linetypes[0])
  ax2 = ax1.twinx() 
  ax2.plot(dicty2[0], dicty2[1], linetypes[1], label=labels[1])
  for tl in ax2.get_yticklabels():
    tl.set_color(linetypes[1])
  ax2.set_ylabel(labels[1], color=linetypes[1])
  plt.savefig('%s.eps' % name)

def make_overview_plot(filename, title, noip_arrs, ip_arrs):
    plt.title("Inner parallelism - " + title)

    
    plt.ylabel('Time (ms)', fontsize=12)

    x = 0
    barwidth = 0.5
    bargroupspacing = 1.5

    for z in zip(noip_arrs, ip_arrs):
        noip,ip = z
        noip_mean,noip_conf = conf_stats(noip)
        ip_mean,ip_conf = conf_stats(ip)

        b_noip = plt.bar(x, noip_mean, barwidth, color='r', yerr=noip_conf, ecolor='black', alpha=0.7)
        x += barwidth

        b_ip = plt.bar(x, ip_mean, barwidth, color='b', yerr=ip_conf, ecolor='black', alpha=0.7)
        x += bargroupspacing

    plt.xticks([0.5, 2.5, 4.5], ['50k', '100k', '200k'], rotation='horizontal')

    fontP = FontProperties()
    fontP.set_size('small')

    plt.legend([b_noip, b_ip], \
        ('no inner parallelism', 'inner parallelism'), \
        prop=fontP, loc='upper center', bbox_to_anchor=(0.5, -0.05), fancybox=True, shadow=True, ncol=2)
   
    plt.ylim([0,62000])
    plt.savefig(output_file(filename))
    plt.clf()

def plot_obs_expc_new(obs, expc, expc_upper, expc_lower, analysis, log, ax = None):
    """Modified version of obs-expc plot suggested by R2. The points are separated by whether their CIs are above, below, 
    
    or overlapping the empirical value
    Input: 
    obs - list of observed values
    expc_mean - list of mean simulated values for the corresponding observed values
    expc_upper - list of the 97.5% quantile of the simulated vlaues
    expc_lower - list of the 2.5% quantile of the simulated values
    analysis - whether it is patitions or compositions
    log - whether the y axis is to be transformed. If True, expc/obs is plotted. If Flase, expc - obs is plotted.
    ax - whether the plot is generated on a given figure, or a new plot object is to be created
    
    """
    obs, expc, expc_upper, expc_lower = list(obs), list(expc), list(expc_upper), list(expc_lower)
    if not ax:
        fig = plt.figure(figsize = (3.5, 3.5))
        ax = plt.subplot(111)
    
    ind_above = [i for i in range(len(obs)) if expc_lower[i] > obs[i]]
    ind_below = [i for i in range(len(obs)) if expc_upper[i] < obs[i]]
    ind_overlap = [i for i in range(len(obs)) if expc_lower[i] <= obs[i] <= expc_upper[i]]
    
    if log:
        expc_standardize = [expc[i] / obs[i] for i in range(len(obs))]
        expc_upper_standardize = [expc_upper[i] / obs[i] for i in range(len(obs))]
        expc_lower_standardize = [expc_lower[i] / obs[i] for i in range(len(obs))]
        axis_min = 0.9 * min([expc_lower_standardize[i] for i in range(len(expc_lower_standardize)) if expc_lower_standardize[i] != 0])
        axis_max = 1.5 * max(expc_upper_standardize)
    else:
        expc_standardize = [expc[i] - obs[i] for i in range(len(obs))]
        expc_upper_standardize = [expc_upper[i] - obs[i] for i in range(len(obs))]
        expc_lower_standardize = [expc_lower[i] - obs[i] for i in range(len(obs))]
        axis_min = 1.1 * min(expc_lower_standardize)
        axis_max = 1.1 * max(expc_upper_standardize)
   
    if analysis == 'partition': col = '#228B22'
    else: col = '#CD69C9'
    ind_full = [] 
    for index in [ind_below, ind_overlap, ind_above]:
        expc_standardize_ind = [expc_standardize[i] for i in index]
        sort_ind_ind = sorted(range(len(expc_standardize_ind)), key = lambda i: expc_standardize_ind[i])
        sorted_index = [index[i] for i in sort_ind_ind]
        ind_full.extend(sorted_index)

    xaxis_max = len(ind_full)
    for i, ind in enumerate(ind_full):
        plt.plot([i, i],[expc_lower_standardize[ind], expc_upper_standardize[ind]], '-', c = col, linewidth = 0.4)
    plt.scatter(range(len(ind_full)), [expc_standardize[i] for i in ind_full], c = col,  edgecolors='none', s = 8)    
    if log: 
        plt.plot([0, xaxis_max + 1], [1, 1], 'k-', linewidth = 1.5)
        ax.set_yscale('log')
    else: plt.plot([0, xaxis_max + 1], [0, 0], 'k-', linewidth = 1.5)
    plt.plot([len(ind_below) - 0.5, len(ind_below) - 0.5], [axis_min, axis_max], 'k--')
    plt.plot([len(ind_below) + len(ind_overlap) - 0.5, len(ind_below) + len(ind_overlap) - 0.5], [axis_min, axis_max], 'k--')
    plt.xlim(0, xaxis_max)
    plt.ylim(axis_min, axis_max)
    plt.tick_params(axis = 'y', which = 'major', labelsize = 8, labelleft = 'on')
    plt.tick_params(axis = 'x', which = 'major', top = 'off', bottom = 'off', labelbottom = 'off')
    return ax