def make_graph(folder, c, i):
	ls = "dotted"
	if "f100" in folder:
		ls = "dashed"
	if "f1000" in folder or "qlearning" in folder:
		ls = "solid" 
	df = pd.read_csv("{}/{}/{}".format(directory_path, folder, data_path))
	x = np.array(df.loc[:, x_axis].values)
	y = np.array(df.loc[:, y_axis].values)
	#print x
	#print y
	if rolling_mean == 1:
		y = pd.Series(y).rolling(window=rolling_width).mean()
		#print y
	split = folder.split("_")
	index = None
	for i in range(len(split)):
		if re.match("f" + r"[0-9]+", split[i]):
			index = i
	if update_freq_label == True and index:
		plt.plot(x, y, color = c, label = re.search(r'[0-9]+', split[index]).group(), linestyle = ls)
	else:
		plt.plot(x, y, color = c, label="{}".format(folder), linestyle = ls)
	if mode == "test" and variance == True:
		reward_std = np.array(df.loc[:, "reward_std"].values)
		plt.fill_between(x, y+reward_std, y-reward_std, facecolor=c, alpha=0.3)

def extractEffectiveMass(df,cut_low=0,cut_right=None,n_cuts=5,makePlot=False):
    if cut_right==None:
        cut_right=max(df["t"])
    window=(cut_right-cut_low)/n_cuts
    
    slope=np.zeros(n_cuts)
    intercept=np.zeros(n_cuts)
    slopeError=np.zeros(n_cuts)
    interceptError=np.zeros(n_cuts)
    
    for i in range(0,n_cuts):
        # select the right intervals for the linear fit
        cut_low_current=cut_low + i*window
        cut_high_current=cut_low + (i+1)*window
        df1=df[(df["t"]>cut_low_current) & (df["t"]<cut_high_current) ]
        if len(df1["t"]) <= 3:
            raise not_enough_data()
        
        params,covs=curve_fit(linear,df1["t"],df1["W"],sigma=df1["deltaW"],maxfev=100000)
        slope[i]=params[0]
        intercept[i]=params[1]
        slopeError[i]=sqrt(covs[0][0])
        interceptError[i]=sqrt(covs[1][1])
        
        if makePlot==True:
            up=(slope[i]+slopeError[i])+(intercept[i]+interceptError[i])/df1["t"]
            down=(slope[i]-slopeError[i])+(intercept[i]-interceptError[i])/df1["t"]
            
            plt.fill_between(df1["t"],up,down,alpha=0.4)
            plt.errorbar(np.array(df1["t"]),np.array(df1["W"])/np.array(df1["t"]),np.array(df1["deltaW"])/np.array(df1["t"]),fmt="or")
    return np.array([slope.mean(),sqrt( slopeError.mean()**2 + slope.var())])

def plot_corridor(x,y,axis=1,alpha=0.2,**kwargs):
    x = np.array(x,dtype=float)
    from matplotlib.pylab import fill_between, plot
    fill_between(x, y.mean(axis)-y.std(axis), y.mean(axis)+y.std(axis),alpha=alpha,**kwargs)
    plot(x, y.mean(axis)+y.std(axis),'-',alpha=alpha,**kwargs)
    plot(x, y.mean(axis)-y.std(axis),'-',alpha=alpha,**kwargs)
    plot(x,y.mean(axis),**kwargs)

def  showMeanSTDMetric(metric_name, m_table, title, output_fn = None):
    pl.figure()
    ax = pl.subplot(111)
    #pl.set_context("talk", font_scale=1, rc={"lines.linewidth": 0.2})
    siz = table.shape[0]
  #  ax=pl.errorbar(range(siz),m_table[metric_name],yerr=m_table[metric_name], fmt='o')

    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)

    ax.get_xaxis().tick_bottom()
    ax.get_yaxis().tick_left()

    pl.fill_between(range(siz), m_table[metric_name+'_mean'] - m_table[metric_name+'_std'],
                m_table[metric_name+'_mean'] + m_table[metric_name+'_std'], color="#3F5D7D",alpha=0.8, lw=0.001)
    pl.plot(range(siz), m_table[metric_name+'_mean'], color="r", lw=2)

    labels = range(1,siz,30)

    #g = sns.FacetGrid(m_table)
    #g.map(pl.errorbar, "csv_path", metric_name, yerr= m_table[metric_name+'_std'] , marker="o")

    pl.xticks(labels,labels,fontsize=9)
    pl.title(title)
    pl.xlabel('Image ID 1 ~ ' + str(siz))
    if ( output_fn is None):
        output_fn = data_DIR+'/'+metric_name+title+'.mean_std.pdf'
    pl.savefig(output_fn)
    pl.show()

def plot_effective_mass(data_file,settings_file="input.xml",label="",param_file="effective_mass_t.dat",factor=1,fit_lines=True):
    root=ET.parse(settings_file).getroot()
    skip=int(root.find("measures").find("winding_number").attrib["skip"])
    time_step=float(root.find("method").find("delta_tau").text)
    fit_mean=[]
    fit_low=[]
    fit_up=[]
    
    data=tools.read_matrix_from_file(data_file)
    data[1]=data[1]*factor
    data[2]=data[2]*factor
    #plt.ylim(0,1)
    if (os.path.isfile(param_file) and fit_lines==True):
        
        with open(param_file,"r") as f:
            cut_low=int(f.readline())
            cut_heigh=int(f.readline())
            params=f.readline().split();
            errors=f.readline().split();
            a=(time_step*(skip+1))
            for x in (data[0][cut_low:cut_heigh]*a):

                #values.append(effective_mass_exp_curve(float(x),float(params[0]),float(params[1]),float(params[2])))
                fit_mean.append(hyperbole(float(x),float(params[0]),float(params[1])*a))
                fit_low.append(hyperbole(float(x),float(params[0])-float(errors[0]),float(params[1])*a - float(errors[1])*a))
                fit_up.append(hyperbole(float(x),float(params[0])+float(errors[0]),float(params[1])*a + float(errors[1])*a))
            plt.plot(data[0][cut_low:cut_heigh]*a,fit_mean,linewidth=4)
            plt.fill_between(data[0][cut_low:cut_heigh]*a,fit_low,fit_up,alpha=0.4)
            
    return plt.errorbar(data[0]*time_step*(skip+1),data[1]/(data[0]*time_step*(skip+1)), yerr=data[2]/(data[0]*time_step*(skip+1)),label=label,fmt='o',ms=1)

def plot_beta_dist( ctr, trials, success, alphas, betas, turns ):
	"""
	Pass in the ctr, trials and success, alphas, betas returned
	by the `experiment` function and the number of turns 
	and plot the beta distribution for all the arms in that turn
	"""
	subplot_num = len(turns) / 2
	x = np.linspace( 0.001, .999, 200 )
	fig = plt.figure( figsize = ( 14, 7 ) ) 

	for idx, turn in enumerate(turns):

		plt.subplot( subplot_num, 2, idx + 1 )

		for i in range( len(ctr) ):
			y = beta( alphas[i] + success[ turn, i ], 
					  betas[i] + trials[ turn, i ] - success[ turn, i ] ).pdf(x)
			line = plt.plot( x, y, lw = 2, label = "arm {}".format( i + 1 ) )
			color = line[0].get_color()
			plt.fill_between( x, 0, y, alpha = 0.2, color = color )
			plt.axvline( x = ctr[i], color = color, linestyle = "--", lw = 2 )
			plt.title("Posteriors After {} turns".format(turn) )
			plt.legend( loc = "upper right" )

	return fig

def nova_plot():

	erg2mev=624151.

	fig=plot.figure()
	yrange = [1e-6,2e-4]
	xrange = [1e-1,1e5]
	plot.fill_between([0.2,10e3],[yrange[1],yrange[1]],[yrange[0],yrange[0]],facecolor='yellow',interpolate=True,color='yellow',alpha=0.5)
	plot.annotate('AMEGO',xy=(3,9e-5),xycoords='data',fontsize=26,color='black')

	lat=ascii.read("data/NMon2012.LAT.dat",names=['energy','en_low','en_high','flux','flux_err','tmp'])
	plot.scatter(lat['energy'],lat['flux']*erg2mev,color='red')
	plot.errorbar(lat['energy'],lat['flux']*erg2mev,xerr=[lat['en_low'],lat['en_high']],yerr=lat['flux_err']*erg2mev,ecolor='red',capsize=0,fmt='none')
	latul=ascii.read("data/NMon2012.LAT.limits.dat",names=['energy','en_low','en_high','flux','tmp1','tmp2','tmp3','tmp4'])
	plot.errorbar(latul['energy'],latul['flux']*erg2mev,xerr=[latul['en_low'],latul['en_high']],yerr=0.5*latul['flux']*erg2mev,uplims=True,ecolor='red',capsize=0,fmt='none')
	plot.scatter(latul['energy'],latul['flux']*erg2mev,color='red')

	leptonic=ascii.read("data/sp-NMon12-IC-best-fit-1MeV-30GeV.txt",names=['energy','flux'],data_start=1)
	hadronic=ascii.read("data/sp-NMon12-pi0-and-secondaries.txt",names=['energy','flux1','flux2'],data_start=1)	

	plot.plot(leptonic['energy'],leptonic['flux']*erg2mev,'r--',color='black',lw=2,label='Leptonic')
	plot.plot(hadronic['energy'],hadronic['flux2']*erg2mev,color='black',lw=2,label='Hadronic+Secondary Leptons')

	plot.legend(loc='upper right',fontsize='small',frameon=False,framealpha=0.5)
	plot.xscale('log')
	plot.yscale('log')
	plot.ylim(yrange)
	plot.xlim(xrange)
	plot.xlabel(r'Energy (MeV)')
	plot.ylabel(r'Energy$^2 \times $ Flux (Energy) (erg cm$^{-2}$ s$^{-1}$)')
	plot.title('Nova V339 Del 2013')
	plot.savefig('Nova_SED.png', bbox_inches='tight')
	plot.savefig('Nova_SED.eps', bbox_inches='tight')
	plot.show()
	plot.close()

def ocean():
    dmax = 1650.e3
    dp = np.linspace(0., dmax, 3301)
    zp = maketopo.shelf1(dp) 
    plt.clf()
    plt.plot(dp*1e-3,zp,'k-',linewidth=3)
    plt.fill_between(dp*1e-3,zp,0.,where=(zp<0.), color=[0,0.5,1])
    plt.xlim([0.,1650])
    plt.ylim([-4500.,500])
    plt.title("Topography as function of radius",fontsize=18)
    plt.xlabel("kilometers from center",fontsize=15)
    plt.ylabel("meters",fontsize=15)
    plt.xticks(fontsize=15)
    plt.yticks(fontsize=15)
    plt.savefig('topo1.png')
    #plt.savefig('topo1.tif')
    print "Cross-section plot in topo1.png"

    plt.xlim([1500,1650])
    plt.ylim([-200.,20])
    plt.title("Topography of shelf and beach",fontsize=18)
    plt.xlabel("kilometers from center",fontsize=15)
    plt.ylabel("meters",fontsize=15)
    plt.xticks(fontsize=15)
    plt.yticks(fontsize=15)
    plt.savefig('topo2.png')
    #plt.savefig('topo2.tif')
    print "Zoom near shore in topo2.png"

def learn_curve_plot(estimator,title,X,y,cv=None,train_sizes=np.linspace(0.1,1.0,5)):
    '''
    :param estimator: the model/algorithem you choose
    :param title: plot title
    :param x: train data numpy array style
    :param y: target data vector
    :param xlim: axes x lim
    :param ylim: axes y lim
    :param cv:
    :return: the figure
    '''
    plt.figure()

    train_sizes,train_scores,test_scores=\
        learning_curve(estimator,X,y,cv=cv,train_sizes=train_sizes)
    '''this is the key score function'''
    train_scores_mean=np.mean(train_scores,axis=1)
    train_scores_std=np.std(train_scores,axis=1)
    test_scores_mean=np.mean(test_scores,axis=1)
    test_scores_std=np.std(test_scores,axis=1)

    plt.fill_between(train_sizes,train_scores_mean - train_scores_std,
                     train_scores_mean + train_scores_std,alpha=0.1,color='b')
    plt.fill_between(train_sizes,test_scores_mean - test_scores_std,
                     test_scores_mean + test_scores_std,alpha=0.1,color='g')
    plt.plot(train_sizes,train_scores_mean,'o-',color='b',label='training score')
    plt.plot(train_sizes,test_scores_mean,'o-',color='g',label='cross valid score')
    plt.xlabel('training examples')
    plt.ylabel('score')
    plt.legend(loc='best')
    plt.grid('on')
    plt.title(title)
    plt.show()

def estimateErrorBlocking(valuesIn,minBlocks=10,makePlot=False):
    
    if len(valuesIn)<minBlocks:
        raise NotEnoughData()
        
    values=valuesIn
    errors=[]
    Ms=[]
    i=len(values)
    
    while i>=minBlocks:
        values=block(values,i)
        errors.append(np.sqrt(np.mean(values**2)-np.mean(values)**2)/sqrt(i))
        i=i/2

    errors=np.array(errors)
    Ms=2**np.linspace(0,len(errors),num=len(errors))
   
    index2=1./np.sqrt(Ms)
    params,covs=optimize.curve_fit(lambda x,a,b: a*x + b,index2,errors)
    res=stats.linregress(index2,errors)
    errorsFit=[np.sqrt(covs[0][0]),np.sqrt(covs[1][1])]
    
    if makePlot:
        x=np.linspace(np.min(Ms),np.max(Ms),num=1000)
        plt.plot(Ms,errors,"o")
        plt.plot(x,params[1] + params[0]/np.sqrt(x))
        y1=params[1]-errorsFit[1] + (params[0]-errorsFit[0])/np.sqrt(x)
        y2=params[1]+errorsFit[1] + (params[0]+errorsFit[0])/np.sqrt(x)
        plt.fill_between(x, y1, y2, where=y2 >= y1, interpolate=True,alpha=0.2)
        plt.plot(x,params[1]+0*x ,"--")
    
    return [params[1],errorsFit[1],abs(params[0]/(params[1]*sqrt(np.max(Ms))))]

	def transition_related_averaging_run(self, simulation_data, smoothing_kernel_width = 200, sampling_interval = [-50, 150], plot = True ):
		"""docstring for transition_related_averaging"""
		transition_occurrence_times = self.transition_occurrence_times(simulation_data = simulation_data, smoothing_kernel_width = smoothing_kernel_width)
		# make sure only valid transition_occurrence_times survive
		transition_occurrence_times = transition_occurrence_times[(transition_occurrence_times > -sampling_interval[0]) * (transition_occurrence_times < (simulation_data.shape[0] - sampling_interval[1]))]
		
		# separated into on-and off periods:
		transition_occurrence_times_separated = [transition_occurrence_times[::2], transition_occurrence_times[1::2]]
		
		mean_time_course, std_time_course = np.zeros((2, sampling_interval[1] - sampling_interval[0], 5)), np.zeros((2, sampling_interval[1] - sampling_interval[0], 5))
		if transition_occurrence_times_separated[0].shape[0] > 2:
			
			for k in [0,1]:
				averaging_interval_times = np.array([transition_occurrence_times_separated[k] + sampling_interval[0],transition_occurrence_times_separated[k] + sampling_interval[1]]).T
				interval_data = np.array([simulation_data[avit[0]:avit[1]] for avit in averaging_interval_times])
				mean_time_course[k] = interval_data.mean(axis = 0)
				std_time_course[k] = (interval_data.std(axis = 0) / np.sqrt(interval_data.shape[0]))
			
			if plot:
				f = pl.figure(figsize = (10,8))
				for i in range(simulation_data.shape[1]):
					s = f.add_subplot(simulation_data.shape[1], 1, 1 + i)
					for j in [0,1]:
						pl.plot(np.arange(mean_time_course[j].T[i].shape[0]), mean_time_course[j].T[i], ['r--','b--'][j], linewidth = 2.0 )
						pl.fill_between(np.arange(mean_time_course[j].shape[0]), mean_time_course[j].T[i] + std_time_course[j].T[i], mean_time_course[j].T[i] - std_time_course[j].T[i], ['r','b'][j], alpha = 0.2)
					s.set_title(self.variable_names[i])
				pl.draw()
			
		return (mean_time_course, std_time_course)

def my_hist(data, bins, range, field, color, label):
    #plt.hist(data, bins, label=label, alpha=0.2, normed=True,
    #         range=range, color=color)
    y, bin_edges = np.histogram(data, bins=bins, normed=True)
    bin_centers = 0.5 * (bin_edges[1:] + bin_edges[:-1])
    plt.fill_between(bin_centers, y, color=color, alpha=0.2, label=label)
    plt.xlabel('$log_{10}($ ' + field+' $)$')
    plt.ylabel('Density')

def plot_particles(PARTICLE_FILENAME, TRUTH_FILENAME, OUT_FILENAME):
    T_DELTA = 1/30.

    a = np.load(PARTICLE_FILENAME)
    truth = pickle.load(open(TRUTH_FILENAME, 'r'))
    weights = a['weights']
    particles = a['particles']
    truth_state = truth['state'][:len(particles)]

    STATEVARS = ['x', 'y', 'xdot', 'ydot', 'phi', 'theta']
    vals = dict([(x, []) for x in STATEVARS])
    for p in particles:
        for v in STATEVARS:
            vals[v].append([s[v] for s in p])
    for v in STATEVARS:
        vals[v] = np.array(vals[v])

    pylab.figure()
    for vi, v in enumerate(STATEVARS):


        if 'dot' in v:
            v_truth = np.diff(truth_state[v[0]])/T_DELTA
        else:
            v_truth = truth_state[v]

        v_bar = np.average(vals[v], axis=1, weights=weights)
        x = np.arange(0, len(v_bar))
        cred = np.zeros((len(x), 2), dtype=np.float)
        for ci, (p, w) in enumerate(zip(vals[v], weights)):
            cred[ci] = util.credible_interval(p, w)

        pylab.subplot(len(STATEVARS) + 1,1, 1+vi)
        pylab.plot(v_truth, color='g', linestyle='-', 
                   linewidth=1)

        pylab.plot(x, v_bar, color='b')
        pylab.fill_between(x, cred[:, 0],
                           cred[:, 1], facecolor='b', 
                           alpha=0.4)


        pylab.ylim([np.min(v_truth), 
                    np.max(v_truth)])
        
    pylab.subplot(len(STATEVARS) + 1, 1, len(STATEVARS)+1)
    # now plot the # of particles consuming 95% of the prob mass
    real_particle_num = []
    for w in weights:
        w = w / np.sum(w) # make sure they're normalized
        w = np.sort(w)[::-1] # sort, reverse order
        wcs = np.cumsum(w)
        wcsi = np.searchsorted(wcs, 0.95)
        real_particle_num.append(wcsi)

    pylab.plot(real_particle_num)

    pylab.savefig(OUT_FILENAME)

def plot_label(label,offset=0):
    Y = df.groupby(cut)[label].mean()
    YERR = df.groupby(cut)[label].std()
    X = df.groupby(cut)["L-cut"].mean()

    color = current_palette.next()
    plt.plot(X+offset,Y,label=label,lw=2,color=color)
    plt.fill_between(X+offset,Y+YERR, Y-YERR,alpha=0.15,color=color)
    return X,Y

def make_test_graph(directory_path, comment):
	df = pd.read_csv("{}/{}/evaluation/evaluation.csv".format(directory_path, comment))
	total_step = np.array(df.loc[:, "total_step"].values)
	reward_mean = np.array(df.loc[:, "reward_mean"].values)
	reward_std = np.array(df.loc[:, "reward_std"].values)
	#step_mean = np.array(df.loc[:, "step_mean"].values, dtype=np.float)
	#step_std = np.array(df.loc[:, "step_std"].values, dtype=np.float)
	plt.figure()
	plt.plot(total_step, reward_mean, color="red")
	plt.fill_between(total_step, reward_mean+reward_std, reward_mean-reward_std, facecolor='red', alpha=0.3)
	plt.savefig("{}/{}/evaluation/reward.png".format(directory_path, comment))

 def plot_auc(self, auc_score, tpr, fpr, label=None):
     pylab.figure(num = None, figsize=(6, 5))
     pylab.xlim([0.0, 1.0])
     pylab.ylim([0.0, 1.0])
     pylab.xlabel('False Positive Rate')
     pylab.ylabel('True Positive Rate')
     pylab.title('ROC (AUC=%0.2f) / %s' % (auc_score, label))
     pylab.fill_between(fpr, tpr, alpha=0.5)
     pylab.grid(True, linestyle='-', color='0.75')
     pylab.plot(fpr, tpr, lw=1)
     pylab.show()

def plot_pr(auc_score, precision, recall, label=None):
    pylab.figure(num=None, figsize=(6, 5))
    pylab.xlim([0.0, 1.0])
    pylab.ylim([0.0, 1.0])
    pylab.xlabel('Recall')
    pylab.ylabel('Precision')
    pylab.title('P/R (AUC=%0.2f) / %s' % (auc_score, label))
    pylab.fill_between(recall, precision, alpha=0.5)
    pylab.grid(True, linestyle='-', color='0.75')
    pylab.plot(recall, precision, lw=1)
    pylab.show()

def plot_pr(auc_score, name, precision, recall, label=None):
    pylab.figure(num=None, figsize=(6, 5))
    pylab.xlim([0.0, 1.0])
    pylab.ylim([0.0, 1.0])
    pylab.xlabel('Recall')
    pylab.ylabel('Precision')
    pylab.title('P/R (AUC=%0.2f) / %s' % (auc_score, label))
    pylab.fill_between(recall, precision, alpha=0.5)
    pylab.grid(True, linestyle='-', color='0.75')
    pylab.plot(recall, precision, lw=1)
    filename = name.replace(" ", "_")
    pylab.savefig(os.path.join(CHART_DIR, "pr_" + filename + ".png"))

def plot_pr(auc_score, name, phase, precision, recall, label=None):
    pylab.clf()
    pylab.figure(num=None, figsize=(5, 4))
    pylab.grid(True)
    pylab.fill_between(recall, precision, alpha=0.5)
    pylab.plot(recall, precision, lw=1)
    pylab.xlim([0.0, 1.0])
    pylab.ylim([0.0, 1.0])
    pylab.xlabel('Recall')
    pylab.ylabel('Precision')
    pylab.title('P/R curve (AUC=%0.2f) / %s' % (auc_score, label))
    filename = name.replace(" ", "_")
    pylab.savefig(os.path.join(CHART_DIR, "pr_%s_%s.png"%(filename, phase)), bbox_inches="tight")

def plot_pr(auc_score, name, phase, precision, recall, label=None):
    from matplotlib import pylab
    pylab.clf()
    pylab.figure(num=None, figsize=(5, 4))
    pylab.grid(True)
    pylab.fill_between(recall, precision, alpha=0.5)
    pylab.plot(recall, precision, lw=1)
    pylab.xlim([0.0, 1.0])
    pylab.ylim([0.0, 1.0])
    pylab.xlabel('Recall')
    pylab.ylabel('Precision')
    pylab.title('P/R curve (AUC=%0.2f) / %s' % (auc_score, label))
    filename = name.replace(" ", "_")