def skymap (tipo1, tipo2="none"):
	if  tipo2 == "none":
		if tipo1 == "casa":
			pl.figure(2, figsize=(6, 5), facecolor='w', edgecolor='k')
			pl.hist2d(ra_casa, dec_casa, bins=51, norm=LogNorm())
			pl.colorbar()
			pl.xlabel('Ascencion Recta [grados]')
			pl.ylabel('Declinacion [grados]')
			pl.show()
		elif tipo1 == "off":
			pl.figure(2, figsize=(6, 5), facecolor='w', edgecolor='k')
			pl.hist2d(ra_off, dec_off, bins=51, norm=LogNorm())
			pl.colorbar()
			pl.xlabel('Ascencion Recta [grados]')
			pl.ylabel('Declinacion [grados]')
			pl.show()
		else: print "no tenemos esos datos"
	elif tipo1 == "casa" and tipo2 == "off":
		np.seterr(divide='ignore', invalid='ignore')
		pl.figure(3, figsize=(5, 5), facecolor='w', edgecolor='k')
		hist_casa, xedge, yedge = np.histogram2d(ra_casa, dec_casa, bins=15)
		hist_off, xedge, yedge = np.histogram2d(ra_off, dec_off, bins=15)
		hist_excess = np.subtract(hist_casa, hist_off)
		hist_excess = np.divide(hist_excess, hist_off)
		pl.imshow(hist_excess, interpolation='gaussian')
		pl.xlabel('Ascencion Recta [u.a.]')
		pl.ylabel('Declinacion [u.a.]')
		pl.show()
	else: print "no tenemos esos datos"

    def histogram(self, index1, index2, **options):
        data = self.datafile.data()

        # param 1
        parameter1 = self.datafile.parameters[index1]
        data1 = data[:, index1]

        parameter2 = self.datafile.parameters[index2]
        data2 = data[:, index2]

        plot.clf()
        plot.figure(figsize=(10, 10), dpi=80)

        # x axis
        plot.xlabel(parameter1[0])
        xmin = options['xmin'] if options['xmin'] != None else parameter1[1]
        xmax = options['xmax'] if options['xmax'] != None else parameter1[2]
        plot.xlim(xmin, xmax)

        # y axis
        plot.ylabel(parameter2[0])
        ymin = options['ymin'] if options['ymin'] != None else parameter2[1]
        ymax = options['ymax'] if options['ymax'] != None else parameter2[2]
        plot.ylim(ymin, ymax)

        # plot
        plot.hist2d(data1, data2, bins=100)
        plot.tight_layout()

        plot.savefig(self.pdffile)

    def _fill_hist(self, x, y, mapsize, data_coords, what='train'):
        x = np.arange(.5, mapsize[1]+.5, 1)
        y = np.arange(.5, mapsize[0]+.5, 1)
        X, Y = np.meshgrid(x, y)

        if what == 'train':
            a = plt.hist2d(x, y, bins=(mapsize[1], mapsize[0]), alpha=.0,
                           cmap=cm.jet)
            area = a[0].T * 12
            plt.scatter(data_coords[:, 1], mapsize[0] - .5 - data_coords[:, 0],
                        s=area.flatten(), alpha=.9, c='None', marker='o',
                        cmap='jet', linewidths=3, edgecolor='r')

        else:
            a = plt.hist2d(x, y, bins=(mapsize[1], mapsize[0]), alpha=.0,
                           cmap=cm.jet, norm=LogNorm())
            area = a[0].T*50
            plt.scatter(data_coords[:, 1] + .5,
                        mapsize[0] - .5 - data_coords[:, 0],
                        s=area, alpha=0.9, c='None', marker='o', cmap='jet',
                        linewidths=3, edgecolor='r')
            plt.scatter(data_coords[:, 1]+.5, mapsize[0]-.5-data_coords[:, 0],
                        s=area, alpha=0.2, c='b', marker='o', cmap='jet',
                        linewidths=3, edgecolor='r')

        plt.xlim(0, mapsize[1])
        plt.ylim(0, mapsize[0])

    def plot(self):
        """
        Gibt die räumliche Verteilung sowie die Energiespektren der Bereiche
        innerhalb eines 4 cm Radius um den Nullpunkt und außerhalb dessen
        in Konsole und Datei aus.
        """
        plt.figure()
        plt.hist2d(self.particles.coords[:, 1], self.particles.coords[:, 2], bins=100)
        plt.colorbar()
        plt.title("Verteilung auf Detektor")
        plt.xlabel("y-Position")
        plt.ylabel("z-Position")
        plt.savefig("distribution.png")

        self.inner = (np.sqrt(np.sum(self.particles.coords[:, 1::] ** 2, 1)) < 40) * self.survivors
        plt.figure()
        plt.hist(self.particles.energy[self.inner] * 1e3, bins=50)
        plt.title("Spektrum in 4 cm Radius")
        plt.xlabel("E / keV")
        plt.ylabel("Anzahl")
        plt.savefig("inner.png")

        plt.figure()
        plt.hist(self.particles.energy[np.logical_not(self.inner)] * 1e3, bins=50)
        plt.title("Spektrum ausserhalb")
        plt.xlabel("E / keV")
        plt.ylabel("Anzahl")
        plt.savefig("outer.png")

def velocity_density_field(P, Rbox, x, y, z, vx, vy, vz, foutname, focusshell=2):
	print 'lets graph! ',foutname
	xmin = x - Rbox
	ymin = y - Rbox
	zmin = z - Rbox
	if (xmin < 0): xmin =0 
	if (ymin < 0): ymin = 0
	if (zmin < 0): zmin = 0
	xmax = x + Rbox
	ymax = y + Rbox
	zmax = z + Rbox
	
	tempx = np.linspace(xmin, xmax, 50)
	tempy = np.linspace(ymin, ymax, 50)
	
	velx = P['v'][:,0] - vx 
	vely = P['v'][:,1] - vy
	
	X,Y = meshgrid(tempx, tempy)
	U =theplt.hist2d(P['p'][:,0], P['p'][:,1], range=[[xmin,xmax],[ymin,ymax]], bins=50, norm=LogNorm(), vmin=1, weights=velx)
	V = theplt.hist2d(P['p'][:,0], P['p'][:,1], range=[[xmin,xmax],[ymin,ymax]], bins=50, norm=LogNorm(), vmin=1, weights=vely)
	M = theplt.hist2d(P['p'][:,0], P['p'][:,1], range=[[xmin,xmax],[ymin,ymax]], bins=50, norm=LogNorm(), vmin=1)
	Q = theplt.quiver( X[::3, ::3], Y[::3, ::3], U[::3, ::3], V[::3, ::3], M[::3, ::3], units='x', pivot='tip', linewidths=(2,), edgecolors=('k'), headaxislength=5, cmap=cm.get_cmap('hot'))
	theplt.colorbar()
	theplt.savefig(foutname)
	theplt.clf()

def run_diagnostics(samples, function=None, plots=True):
    if plots:
        xlim = (-0.5, 1.5)
        ylim = (-1.5, 1.)
        
        # plot the sample distribution
        f = plt.gcf()
        f.set_size_inches(8, 8)
        plt.hist2d(samples[:,1], samples[:,0], bins=50, cmap=reds, zorder=100)
        plt.xlabel("Longitude")
        plt.ylabel("Latitude")
        
        # overlay the true function
        if function:
            plot_true_function(function, xlim, ylim)
        
        plt.show()
        
        plot_diagnostics(samples)
    
    gelman_rubin(samples)

    # gewecke
    #geweke_val = pymc.diagnostics.geweke(samples, intervals=1)[0][0][1]
    Geweke(samples)

def partD():
    ls = np.linspace(-math.pi/2, math.pi/2, DATA_DIM*2)
    rs = np.linspace(0.001, 25, DATA_DIM)
    
    l_list = []
    vel_list = []

    for i in range(len(ls)):
        for j in range(len(rs)):
            R = math.sqrt(rs[j]**2 + 10**2 - 2 * 10 * rs[j] * math.cos(ls[i]))
            if (R < 15):
                if (R > 4 and R < 6):
                    vel_add = 50*math.cos( math.asin( 10*(math.sin(abs(ls[i]))/R )) )
                    #vel_add = 0 #subbing out radial expansion
                    vel = vel4(rs[j], ls[i], R) - vel_add*(ls[i]/abs(ls[i]))
                    l_list.append(ls[i])
                    l_list.append(ls[i])
                    vel_list.append(vel)
                    vel_list.append(vel)
                else:
                    l_list.append(ls[i])
                    vel_list.append(vel4(rs[j], ls[i], R))

    plt.hist2d(l_list, vel_list, PLOT_DIM, cmin=1)
    plt.savefig("partd.png")

def plot_replay_memory_2d_state_histogramm(states):
    x,v = zip(*states)
    plt.hist2d(x, v, bins=40, norm=LogNorm())
    plt.xlabel("position")
    plt.ylabel("velocity")
    plt.colorbar()
    plt.show()

def histo2d(r1, r2):
    nBin = 8
    h1 = r1*nBin/255.0
    h2 = r2*nBin/255.0
    plt.hist2d(h1, h2, bins=nBin+1, norm=LogNorm())
    plt.colorbar()
    plt.show()

def plotMAtreeFINALhist( scale1, mass1, Vmaxarray, plotfunc_count):
	print("Entered Plot Function")
	plot_title="Mass Accretion History Histogram"   #Can code the number in with treemax
	x_axis="scale time"
	y_axis="total mass"
	figure_name=os.path.expanduser('~/Jan14HISTMassAccretionfigure' +'.png')
	#Choose which type of plot you would like: Commented out.
	
	plt.hist2d(scale1, mass1, (100, 100), cmap=plt.cm.jet)

	#plt.plot(scale1, mass1, linestyle="", marker="o", markersize=3, edgecolor='none')
	#plt.plot(scale1, Vmaxarray, linestyle="-", marker="o")
	

	#plt.scatter(scale1, mass1, label="first tree")
	plt.title(plot_title)
	plt.xlabel(x_axis)
	plt.ylabel(y_axis)
	#plt.yscale("log")

	plt.savefig(figure_name)
	print("Saving plot: %s" % figure_name)
	print
	
	#In order to Plot only a single tree on a plot must clear lists before loop. 
	#Comment out to over plot curves.			
	plt.clf()

	clearmass  = []
	clearscale = []
	clearVmax  = []
	plotfunc_count += 1
	#return clearmass, clearscale, clearVmax, plotfunc_count
	return

def PlotHeatMap(Bins,X,Y,ContourLengthNm,
                LabelColorbar="Data points in bin"):
    plt.hist2d(X,Y,bins=Bins,cmap='afmhot')
    plt.axhline(65,color='g',linestyle='--',linewidth=4,alpha=0.3,
                label="65pN")
    plt.axvline(ContourLengthNm,color='c',linestyle='--',linewidth=4,
                alpha=0.3,label=r"L$_0$={:.1f}nm".format(ContourLengthNm)) 

def plot_impact_map(impactX, impactY, telX, telY, telTypes=None, Outfile="ImpactMap.png"):
    """
    Map of the site with telescopes positions and impact points heatmap

    Parameters
    ----------
    impactX: `numpy.ndarray`
    impactY: `numpy.ndarray`
    telX: `numpy.ndarray`
    telY: `numpy.ndarray`
    telTypes: `numpy.ndarray`
    Outfile: string - name of the output file
    """
    plt.figure(figsize=(12, 12))

    plt.hist2d(impactX, impactY, bins=40)
    plt.colorbar()

    assert (len(telX) == len(telY)), "telX and telY should have the same length"
    if telTypes:
        assert (len(telTypes) == len(telX)), "telTypes and telX should have the same length"
        plt.scatter(telX, telY, color=telTypes, s=30)
    else:
        plt.scatter(telX, telY, color='black', s=30)

    plt.axis('equal')
    plt.savefig(Outfile, bbox_inches="tight", format='png', dpi=200)
    plt.close()

def main():

    # Parse commandline arguments
    parser = argparse.ArgumentParser(description="Runs benchmark that has been compiled with contech, generating a task graph and optionally running backend tools.")
    parser.add_argument("inFile", help="Input file, a json file of commRecords")
    parser.add_argument("outFile", help="Output file, a png of inter-thread communication.")
    args = parser.parse_args()
        
    x, y = [], []
    nThreads = 0
    print "Loading {}...".format(args.inFile)
    with open(args.inFile) as file:
        data = json.load(file)
        records = data["records"]
        print "Loaded {} records".format(len(records))
        
        for r in records:
            src = int(r["src"].split(":")[0])
            dst = int(r["dst"].split(":")[0])
            nThreads = max([nThreads, src, dst])
            x.append(dst)
            y.append(src)
            
    print "Plotting for {} threads, {} communication records...".format(nThreads, len(x))
    
#     plt.title(os.path.basename(args.inFile).replace(".json",""))
    plt.xlabel("consumer CPU")    
    plt.ylabel("producer CPU")
    tickPositions = [a + .5 for a in range(nThreads)]
    plt.xticks(tickPositions, range(nThreads))
    plt.yticks(tickPositions, range(nThreads))
    plt.hist2d(x, y, bins=nThreads, cmap=matplotlib.cm.Greys)
    plt.colorbar()
    plt.savefig(args.outFile)        

def boundary_sim(xyini, exyini, a, xy_step, dt, tmax, expmt, ephi, vb):
	"""
	Run the LE simulation from (x0,y0), stopping if x<xmin.
	"""
	me = me0+".boundary_sim: "
		
	## Initialisation
	x0,y0 = xyini
	nstp = int(tmax/dt)
	
	## Simulate eta
	if vb: t0 = time.time()
	exy = np.vstack([sim_eta(exyini[0], expmt, nstp, a, dt), sim_eta(exyini[1], expmt, nstp, a, dt)]).T
	if vb: print me+"Simulation of eta",round(time.time()-t0,2),"seconds for",nstp,"steps"
				
	## Spatial variables
	if vb: t0 = time.time()
	xy = np.zeros([nstp,2]); xy[0] = [x0,y0]
	
	## Calculate trajectory
	for i in xrange(0,nstp-1):
		r = np.sqrt((xy[i]*xy[i]).sum())
		xy[i+1] = xy[i] + xy_step(xy[i],r,exy[i])
		
	if vb: print me+"Simulation of x",round(time.time()-t0,2),"seconds for",nstp,"steps"
	
	rcoord = np.sqrt((xy*xy).sum(axis=1))
	ercoord = np.sqrt((exy*exy).sum(axis=1))
	
	## -----------------===================-----------------
	R = 2.0; S = 1.0; lam = 0.5; nu = 1.0
	## Distribution of spatial steps and eta
	if 0:
		from LE_RunPlot import plot_step_wall, plot_eta_wall, plot_step_bulk
		## In wall regions
		plot_step_wall(xy,rcoord,R,S,a,dt,vb)
		plot_eta_wall(xy,rcoord,exy,ercoord,R,S,a,dt,vb)
		## In bulk region
		plot_step_bulk(xy,rcoord,ercoord,R,S,a,dt,vb)
		exit()
	## Trajectory plot with force arrows
	if 0:
		from LE_RunPlot import plot_traj	
		plot_traj(xy,rcoord,R,S,lam,nu,force_dnu,a,dt,vb)
		exit()
	## 2D Density plot
	if 0:
		R=5.0
		plt.hist2d(xy[:,0],xy[:,1],1000,cmap="Blues")
		ang = np.linspace(0,2*np.pi,60)
		plt.plot(R*np.cos(ang),R*np.sin(ang),"g--",lw=2)
		plt.show()
		exit()
	## -----------------===================-----------------
			
	if ephi:
		epcoord = np.arctan2(exy[:,1],exy[:,0]) - np.arctan2(xy[:,1],xy[:,0])
		return [rcoord, ercoord, epcoord]
	else:
		return [rcoord, ercoord]

def group_mark_pos(step, per_step):
    i = 0
    while i<=step:
        x = []
        y = []
        filename = "./data/gif_dat/mark_pos_metropolis_"+str(i)+".dat"
        with open(filename,"r") as f:
            line = f.readline().split()
            size = line[1]
            line = f.readline().split()
            T = line[1]
            line = f.readline().split()
            J = line[1]
            line = f.readline().split()
            H = line[1]
            tail = "_"+size+"_T"+T+"_J"+J+"_H"+H
            line = f.readline().split()
            while line:
                line = f.readline().split()
                if len(line)==2:
                    x.append(int(line[0]))
                    y.append(int(line[1]))
        fig = plt.figure(figsize = (5,5))
        plt.title("mark_pos")
        plt.xlabel("x")
        plt.ylabel("y")
        plt.hist2d(x,y,bins = int(size))
        plt.xlim(0,int(size)-1)
        plt.ylim(0,int(size)-1)
        #plt.savefig(filename[:-15]+tail+".png")
        plt.savefig(filename[:-4]+".png")
        plt.close()
        i += per_step
        if (step//i != step//(i+1)):
            print(step//i)

def draw_scatter_2(obj, varname1, varname2, loclist, idlist, stoparray):
    """
    TODO
    """
    
    varlist1 = np.array([])
    varlist2 = np.array([])
    
    
    cnt = 0   # continuous counter for startarray and stoparray
    for i in range(len(loclist)):
        print 'Plotting file', i+1, 'of', len(loclist)
        rootgrp = nc.Dataset(loclist[i], 'r')
        
        var1 = rootgrp.variables[varname1][:, :]
        var2 = rootgrp.variables[varname2][:, :]
        p = rootgrp.variables['P'][:, :]
        
        for j in idlist[i]:
            mask = (np.isfinite(p[:, j])) & (p[:, j] != 0)
            mask &= (np.isfinite(var1[:, j])) & (np.isfinite(var2[:, j]))
            if not stoparray == None:
                mask[:stoparray[cnt]] = False
            varlist1 = np.append(varlist1, var1[:, j][mask])
            varlist2 = np.append(varlist2, var2[:, j][mask])
            cnt += 1
    
    varlist2 *= 1.e6
    #plt.scatter(varlist1, varlist2)
    x_edges = np.arange(0, 1000 + 20, 20)
    y_edges = np.arange(-20, 20 + 0.8, 0.8)
    plt.hist2d(varlist1, varlist2, cmin = 1, 
               norm = clr.LogNorm(1, 10000),
               bins = [x_edges, y_edges])
    plt.colorbar()

def generate_beam():
    plt.clf()
    
    #grab data
    x = data_io.wireX
    y = data_io.wireY
    
    if len(x) == 0 or len(y) == 0:
        return None, None

    #make plot
    ax1 = plt.subplot2grid((3,3),(0,0), colspan=2)
    plt.hist(x, bins=50, range=[-5,5])
    plt.xlim(-5,5)
    ax1.axes.get_xaxis().set_visible(False)
    ax1.axes.get_yaxis().set_visible(False)
    ax2 = plt.subplot2grid((3,3),(1,0), colspan=2, rowspan=2)
    plt.hist2d(x, y, bins=50, range=np.array([(-5,5),(-5,5)]))
    plt.xlim(-5,5)
    plt.ylim(-5,5)
    plt.xlabel('x position')
    plt.ylabel('y position')
    ax3 = plt.subplot2grid((3,3),(1,2), rowspan=2)
    plt.hist(y, bins=50, orientation='horizontal', range=[-5,5])
    plt.ylim(-5,5)
    ax3.axes.get_xaxis().set_visible(False)
    ax3.axes.get_yaxis().set_visible(False)
    
    for temp_file in glob.glob(app.config['UPLOAD_FOLDER'] + '/temp_beam*'):
        os.remove(temp_file)
    filename = unique_filename('temp_beam.png')
    filepath = upload_path(filename)
    plt.savefig(filepath)

    return filename, filepath

def doit(_DATADIR, _HANDLE):

    fullData = pd.DataFrame({'minutes':[],'fly_x':[],'fly_y':[],'FlyID':[],'GROUP':[]})

    for directory in glob.glob(_DATADIR + '*' + _HANDLE + '*' + '*zoom*'):
        FLY_ID, FMF_TIME, GROUP = parse_fmftime(directory)
        try:
            df = pd.read_pickle(directory + '/frame_by_frame_synced.pickle')
            df['minutes'] = df.synced_time.astype('timedelta64[m]')
            tempdf = pd.DataFrame({'minutes':df.minutes, 'fly_x':df.fly_x, 'fly_y':df.fly_y, 'FlyID':FLY_ID, 'group':GROUP})
            fullData = pd.concat([fullData, tempdf], axis=0)
        except:
            pass

    fullData.to_pickle(_DATADIR + '2D_positions_' + _HANDLE + '.pickle')

    plt.figure(figsize=(24,2))

    for x in np.arange(-1,11):
        plt.subplot2grid((1,49),(0,4*(x+1)), colspan=4)
        slyce = fullData[fullData.minutes == x]
        print x, len(slyce), len(slyce.fly_x.notnull())
        plt.hist2d(slyce.fly_x.values, slyce.fly_y.values, bins=50, norm=LogNorm())
        plt.xlim(0,659)
        plt.ylim(0,494)
        plt.axis('off')
        if x == 10:
            axes = plt.subplot2grid((1,49),(0,48), colspan=1)
            plt.colorbar(cax=axes)

    plt.savefig(_DATADIR + 'position_plots_' + _HANDLE + '.svg', bbox_inches='tight')
    return

def survey_overlap(xs, ys, xs1, ys1, c_graph = True, radius = -0.1, nbins = 1000):

	"""
	Takes two sets of points and sees if one is inside the other
	Might need nbins and radius tweaking
	xs, ys are the points to draw the contours off
	xs1, ys1 are the points to check if in contour
	"""

	import numpy as np
	import matplotlib.pyplot as plt

	xwalls = np.linspace( min(xs) - 5.0, max(xs) + 5.0, nbins + 1 )
	ywalls = np.linspace( min(ys) - 5.0, max(ys) + 5.0, nbins + 1 )

	im, xs_bin, ys_bin, ax = plt.hist2d(xs, ys, bins = (xwalls, ywalls) )
	xs_mids = 0.5*(xs_bin[:-1] + xs_bin[1:])
	ys_mids = 0.5*(ys_bin[:-1] + ys_bin[1:])
	plt.close()
	im[im>0] = 1
	conts = plt.contour(xs_mids, ys_mids, im.T, 1)
	
	xwalls1 = np.linspace( min(xs1) - 5.0, max(xs1) + 5.0, nbins + 1)
	ywalls1 = np.linspace( min(ys1) - 5.0, max(ys1) + 5.0, nbins + 1)
	im1, xs_bin1, ys_bin1, ax1 = plt.hist2d(xs1, ys1, bins = (xwalls1, ywalls1))
	xs_mids1 = 0.5*(xs_bin1[:-1] + xs_bin1[1:])
	ys_mids1 = 0.5*(ys_bin1[:-1] + ys_bin1[1:])
	plt.close()
	im1[im1>0] = 1

	conts = plt.contour(xs_mids, ys_mids, im.T, 1)
	conts1 = plt.contour(xs_mids1, ys_mids1, im1.T, 1)
	plt.show()

def hist2d(data, bins):
    data = np.array(data)
    x = data[:,0]
    y = data[:,1]

    plt.hist2d(x, y, bins)
    plt.colorbar()