def test():
	figure(figsize=(6,12))
	t = linspace(0, 1, 1001)
	y = sin(2*pi*t*6) + sin(2*pi*t*10) + sin(2*pi*t*13)
	subplot(311)
	plot(t, y, 'b-')
	xlabel("TIME (sec)")
	ylabel("SIGNAL MAGNITUDE")
	# compute FFT and plot the magnitude spectrum
	F = fft(y)
	N = len(t)             # number of samples
	dt = 0.001             # inter-sample time difference
	w = fftfreq(N, dt)     # gives us a list of frequencies for the FFT
	ipos = where(w>0)
	freqs = w[ipos]        # only look at positive frequencies
	mags = abs(F[ipos])    # magnitude spectrum
	subplot(312)
	plot(freqs, mags, 'b-')
	ylabel("POWER")
	subplot(313)
	plot(freqs, mags, 'b-')
	xlim([0, 50])          # replot but zoom in on freqs 0-50 Hz
	ylabel("POWER")
	xlabel("FREQUENCY (Hz)")
	savefig("signal_3freqs.jpg", dpi=150)

def test_radar_view():
	lats=linspace(-13.0, -11.5, 40)
	lons=linspace(130., 132., 40)
	ber_loc=[-12.4, 130.85] #location of Berrimah radar
	gp_loc=[-12.2492,  131.0444]#location of CPOL at Gunn Point
	h=2.5*1000.0
	t0=systime()
	i_a, j_a, k_a=propigation.unit_vector_grid(lats, lons, h, gp_loc)
	print "unit_vector_compute took ", systime()-t0, " seconds"
	fw=0.1#degrees
	m_bump=5.0
	#b1=simul_winds.speed_bump(lats, lons, [-12.0, 131.0], fw)*m_bump
	#b2=-1.0*simul_winds.speed_bump(lats, lons, [-12.0, 131.1], fw)*m_bump
	u,v=simul_winds.unif_wind(lats, lons, 5.0, 75.0)
	up,vp=array(simul_winds.vortex(lats, lons,  [-12.5, 131.1], fw))*m_bump
	#up=u+(b1-b2)
	#vp=v-(b1-b2)
	w=v*0.0
	v_r=i_a*(up+u)+j_a*(vp+v)+k_a*w
	f=figure()
	mapobj=pres.generate_darwin_plot()
	pres.contour_vr(mapobj,lats, lons, v_r)
	savefig(os.getenv('HOME')+'/bom_mds/output/test_radar_view_gp.png')
	close(f)
	f=figure()
	mapobj=pres.generate_darwin_plot()
	pres.quiver_contour_winds(mapobj, lats, lons, (up+u),(vp+v))
	savefig(os.getenv('HOME')+'/bom_mds/output/test_pert2.png')
	close(f)

def show_dct_fig():
    figure(figsize=(12,7))
    for u in range(8):
        subplot(2, 4, u+1)
        ylim((-1, 1))
        title(str(u))
        plot(arange(0,8,0.1), dct0[u, :])
        plot(dct[u, :],'ro')

def test_gracon():
	#setup
	noise_level=0.0#m/s
	nx=40
	ny=40
	fw=0.1
	m_bump=10.00
	t0=systime()
	lats=linspace(-13.5, -12.0, 40)
	lons=linspace(130.5, 131.5, 40)
	ber_loc=[-12.4, 130.85] #location of Berrimah radar
	gp_loc=[-12.2492,  131.0444]#location of CPOL at Gunn Point
	h=2.5*1000.0
	print 'calculating berimah UV', systime()-t0
	i_ber, j_ber, k_ber=propigation.unit_vector_grid(lats, lons, h, ber_loc)
	print 'calculating gp UV', systime()-t0
	i_gp, j_gp, k_gp=propigation.unit_vector_grid(lats, lons, h, gp_loc)
	#make winds
	u,v=simul_winds.unif_wind(lats, lons, 10.0, 75.0)
	up,vp=array(simul_winds.vortex(lats, lons,  [-12.5, 131.1], fw))*m_bump
	#make V_r measurements
	vr_ber=i_ber*(up+u)+j_ber*(vp+v) + (random.random([nx,ny])-0.5)*(noise_level*2.0)
	vr_gp=i_gp*(up+u)+j_gp*(vp+v)+ (random.random([nx,ny])-0.5)*(noise_level*2.0)
	#try to reconstruct the wind field
	igu, igv= simul_winds.unif_wind(lats, lons, 0.0, 90.0)
	gv_u=zeros(u.shape)
	gv_v=zeros(v.shape)
	f=0.0
	print igu.mean()
	
	angs=array(propigation.make_lobe_grid(ber_loc, gp_loc, lats,lons))
	wts=zeros(angs.shape, dtype=float)+1.0
	#for i in range(angs.shape[0]):
	#	for j in range(angs.shape[1]):
	#		if (angs[i,j] < 150.0) and (angs[i,j] > 30.0): wts[i,j]=1.0
	print 'Into fortran'
	gv_u,gv_v,f,u_array,v_array = gracon_vel2d.gracon_vel2d(gv_u,gv_v,f,igu,igv,i_ber,j_ber,i_gp,j_gp,vr_ber,vr_gp,wts, nx=nx,ny=ny)
	print u_array.mean()
	print f
	bnds=[0.,20.]
	f=figure()
	mapobj=pres.generate_darwin_plot()
	pres.quiver_contour_winds(mapobj, lats, lons, (up+u),(vp+v), bounds=bnds)
	savefig(os.getenv('HOME')+'/bom_mds/output/orig_winds_clean.png')
	close(f)
	f=figure()
	mapobj=pres.generate_darwin_plot()
	pres.quiver_contour_winds(mapobj, lats, lons, (wts*u_array +0.001),(wts*v_array +0.001), bounds=bnds)
	savefig(os.getenv('HOME')+'/bom_mds/output/recon_winds_clean.png')
	close(f)
	f=figure()
	mapobj=pres.generate_darwin_plot()
	pres.quiver_contour_winds(mapobj, lats, lons, (wts*u_array - (up+u)),(wts*v_array -(vp+v)))
	savefig(os.getenv('HOME')+'/bom_mds/output/errors_clean.png')
	close(f)

def plotSolutions(x,states=None): #Default func values is trivial

    plt.figure(figsize=(11,8.5))

    #get the exact values
    f = open('exact_results.txt', 'r')
    x_e = []
    u_e = []
    p_e = []
    rho_e = []
    e_e = []
    for line in f:
        if len(line.split())==1:
            t = line.split()
        else:
            data = line.split()
            x_e.append(float(data[0]))
            u_e.append(float(data[1]))
            p_e.append(float(data[2]))
            rho_e.append(float(data[4]))
            e_e.append(float(data[3]))


    if states==None:
        raise ValueError("Need to pass in states")
    else:
        u = []
        p = []
        rho = []
        e = []
        for i in states:
            u.append(i.u)
            p.append(i.p)
            rho.append(i.rho)
            e.append(i.e)

    #get edge values
    x_cent = [0.5*(x[i]+x[i+1]) for i in range(len(x)-1)]
    
    if u != None:
        plot2D(x_cent,u,"$u$",x_ex=x_e,y_ex=u_e)
    
    if rho != None:
        plot2D(x_cent,rho,r"$\rho$",x_ex=x_e,y_ex=rho_e) 

    if p != None:
        plot2D(x_cent,p,r"$p$",x_ex=x_e,y_ex=p_e)

    if e != None:
        plot2D(x_cent,e,r"$e$",x_ex=x_e,y_ex=e_e)

    plt.show(block=False) #show all plots generated to this point
    raw_input("Press anything to continue...")
    plot2D.fig_num=0

def generate_plot(array, vmin, vmax, figNumber=1):
    plt.figure(figNumber)
    plt.subplot(2,3,i)
    print i
    plt.imshow(array, vmin = vmin, vmax= vmax, interpolation = None) 
    plt.xlabel('Sample')
    plt.ylabel('Line')
    plt.title(row[0])
    cb = plt.colorbar(orientation='hor', spacing='prop',ticks = [vmin,vmax],format = '%.2f')
    cb.set_label('Reflectance / cos({0:.2f})'.format(incAnglerad*180.0/math.pi))
    plt.grid(True)

	def plot_f_score(self, disag_filename):
		plt.figure()
		from nilmtk.metrics import f1_score
		disag = DataSet(disag_filename)
		disag_elec = disag.buildings[building].elec
		f1 = f1_score(disag_elec, test_elec)
		f1.index = disag_elec.get_labels(f1.index)
		f1.plot(kind='barh')
		plt.ylabel('appliance');
		plt.xlabel('f-score');
		plt.title(type(self.model).__name__);

def compress_kmeans(im, k=4):
    height, width, depth = im.shape

    data = im.reshape((height * width, depth))
    labels, centers = kmeans(data, k, 1e-2)
    rep = closest(data, centers)
    data_compressed = centers[rep]

    im_compressed = data_compressed.reshape((height, width, depth))
    plt.figure()
    plt.imshow(im_compressed)
    plt.show()

 def plot(self, output):
     plt.figure(figsize=output.fsize, dpi=output.dpi)
     for ii in range(0, len(self.v)):
         imsize = [self.t[0], self.t[-1], self.x[ii][-1], self.x[ii][0]]
         lim = amax(absolute(self.v[ii])) / output.scale_sat
         plt.imshow(self.v[ii], extent=imsize, vmin=-lim, vmax=lim, cmap=cm.gray, origin='upper', aspect='auto')
         plt.title("%s-Velocity for Trace #%i" % (self.comp.upper(), ii))
         plt.xlabel('Time (s)')
         plt.ylabel('Offset (km)')
         #plt.colorbar()
         plt.savefig("Trace_%i_v%s.pdf" % (ii, self.comp))
         plt.clf()

def plotCentroidFitDiagnostic(img, hdr, ccdMod, ccdOut, res, prfObj):
    """Some diagnostic plots showing the performance of fitPrfCentroid()

    Inputs:
    -------------
    img
        (np 2d array) Image of star to be fit. Image is in the
        format img[row, col]. img should not contain Nans

    hdr
        (Fits header object) header associated with the TPF file the
        image was drawn from

    ccdMod, ccdOut
        (int) CCD module and output of image. Needed to
        create the correct PRF model

    prfObj
        An object of the class prf.KeplerPrf()


    Returns:
    -------------
    **None**

    Output:
    ----------
    A three panel subplot is created
    """
    mp.figure(1)
    mp.clf()
    mp.subplot(131)
    plotTpf.plotCadence(img, hdr)
    mp.colorbar()
    mp.title("Input Image")

    mp.subplot(132)
    c,r = res.x[0], res.x[1]
    bbox = getBoundingBoxForImage(img, hdr)
    model = prfObj.getPrfForBbox(ccdMod, ccdOut, c, r, bbox)
    model *= res.x[2]
    plotTpf.plotCadence(model, hdr)
    mp.colorbar()
    mp.title("Best fit model")

    mp.subplot(133)
    diff = img-model
    plotTpf.plotCadence(diff, hdr)
    mp.colorbar()
    mp.title("Residuals")

    print "Performance %.3f" %(np.max(np.abs(diff))/np.max(img))

def plot_confusion_matrix(cm, labels, title='Confusion matrix', cmap=plt.cm.Blues, save=False):
    plt.figure()
    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    plt.title(title)
    plt.colorbar()
    tick_marks = np.arange(len(labels))
    plt.xticks(tick_marks, labels, rotation=45)
    plt.yticks(tick_marks, labels)
    plt.tight_layout()
    plt.ylabel('True label')
    plt.xlabel('Predicted label')
    plt.show()
    if save:
        plt.savefig(save)

def get_trajectories_for_DF(DF):
    GetXYS=ct.get_xys(DF)
    xys_s=ct.get_xys_s(GetXYS['xys'],GetXYS['Nmin'])
    plt.figure(figsize=(5, 5),frameon=False)
    for m in list(range(9)):
        plt.plot()
        plt.subplot(3,3,m+1)
        xys_s_x_n=xys_s[m]['X']-min(xys_s[m]['X'])
        xys_s_y_n=xys_s[m]['Y']-min(xys_s[m]['Y'])
        plt.plot(xys_s_x_n,xys_s_y_n)
        plt.axis('off')
        axes = plt.gca()
        axes.set_ylim([0,125])
        axes.set_xlim([0,125])

def get_trajectories_singleAxis_for_DF(DF):
    from random import randint
    sns.set_palette(sns.color_palette("Paired"))
    fig = plt.figure(figsize=(5, 5),frameon=False)
    ax=fig.add_subplot(1, 1, 1)
    ax.spines['left'].set_position('zero')
    ax.spines['right'].set_color('none')
    ax.spines['bottom'].set_position('zero')
    ax.spines['top'].set_color('none')
    #ax.spines['left'].set_smart_bounds(True)
    #ax.spines['bottom'].set_smart_bounds(True)
    ax.xaxis.set_ticks_position('bottom')
    ax.yaxis.set_ticks_position('left')
    ax.set_ylim([-300,300])
    ax.set_xlim([-300,300])
    ticklab = ax.xaxis.get_ticklabels()[0]
    ax.xaxis.set_label_coords(300, -40,transform=ticklab.get_transform())
    ax.set_xlabel('x($\mu$m)',fontsize=14)
    ticklab = ax.yaxis.get_ticklabels()[0]
    ax.yaxis.set_label_coords(90, 280,transform=ticklab.get_transform())
    ax.set_ylabel('y($\mu$m)',rotation=0,fontsize=14)
    
    
    GetXYS=ct.get_xys(DF)
    xys_s=ct.get_xys_s(GetXYS['xys'],GetXYS['Nmin'])
    for n in list(range(12)):
        m=randint(0,len(xys_s))-1
        xys_s_x_n=xys_s[m]['X']-(xys_s[m]['X'][xys_s[m]['X'].index[0]])
        xys_s_y_n=xys_s[m]['Y']-(xys_s[m]['Y'][xys_s[m]['X'].index[0]])
        xys_s_x_n=[x*(100/(60.)) for x in xys_s_x_n]
        xys_s_y_n=[x*(100/(60.)) for x in xys_s_y_n]
        ax.plot(xys_s_x_n,xys_s_y_n)
    return fig

def plotRetinaSpikes(retina=None, label=""):
    
    assert retina is not None, "Network is not initialised! Visualising failed."
    import matplotlib.pyplot as plt
    from matplotlib import animation
    
    print "Visualising {0} Spikes...".format(label) 

    spikes = [x.getSpikes() for x in retina]
#     print spikes
    
    sortedSpikes = sortSpikes(spikes)
#     print sortedSpikes
    
    framesOfSpikes = generateFrames(sortedSpikes)
#     print framesOfSpikes
    
    x = range(0, dimensionRetinaX)
    y = range(0, dimensionRetinaY)
    from numpy import meshgrid
    rows, pixels = meshgrid(x,y)
    
    fig = plt.figure()
    
    initialData = createInitialisingData()
    
    imNet = plt.imshow(initialData, cmap='green', interpolation='none', origin='upper')
    
    plt.xticks(range(0, dimensionRetinaX)) 
    plt.yticks(range(0, dimensionRetinaY))
    args = (framesOfSpikes, imNet)
    anim = animation.FuncAnimation(fig, animate, fargs=args, frames=int(simulationTime)*10, interval=30)
          
    plt.show()

def plotColorCodedNetworkSpikes(network):
    assert network is not None, "Network is not initialised! Visualising failed."
    import matplotlib as plt
    from NetworkBuilder import sameDisparityInd
    
    cellsOutSortedByDisp = []
    spikes = []
    for disp in range(0, maxDisparity+1):
        cellsOutSortedByDisp.append([network[x][2] for x in sameDisparityInd[disp]])
        spikes.append([x.getSpikes() for x in cellsOutSortedByDisp[disp]])
    
    sortedSpikes = sortSpikesByColor(spikes)
    print sortedSpikes
    framesOfSpikes = generateColoredFrames(sortedSpikes)
    print framesOfSpikes
    
    fig = plt.figure()
    
    initialData = createInitialisingDataColoredPlot()
    
    imNet = plt.imshow(initialData[0], c=initialData[1], cmap=plt.cm.coolwarm, interpolation='none', origin='upper')
    
    plt.xticks(range(0, dimensionRetinaX)) 
    plt.yticks(range(0, dimensionRetinaY))
    plt.title("Disparity Map {0}".format(disparity))
    args = (framesOfSpikes, imNet)
    anim = animation.FuncAnimation(fig, animate, fargs=args, frames=int(simulationTime)*10, interval=30)
          
    plt.show()

def plot_images(header):
    ''' function to plot images from header.

    It plots images, return nothing
    Parameters
    ----------
        header : databroker header object
            header pulled out from central file system
    '''
    # prepare header
    if type(list(headers)[1]) == str:
        header_list = list()
        header_list.append(headers)
    else:
        header_list = headers

    for header in header_list:
        uid = header.start.uid
        img_field = _identify_image_field(header)
        imgs = np.array(get_images(header, img_field))
        print('Plotting your data now...')
        for i in range(imgs.shape[0]):
            img = imgs[i]
            plot_title = '_'.join(uid, str(i))
            # just display user uid and index of this image
            try:
                fig = plt.figure(plot_title)
                plt.imshow(img)
                plt.show()
            except:
                pass # allow matplotlib to crash without stopping other function

def showScatterPlot(data, labels, idx1, idx2): 
    import matplotlib.pyplot as plt 
    fig = plt.figure() 
    ax = fig.add_subplot(111) 
    # X-axis data, Y-axis data, Size for each sample, Color for each sample 
    ax.scatter(data[:,idx1], data[:,idx2], 100.0*(1 + np.array(labels)), 100.0*(1 + np.array(labels))) 
    plt.show()

def showScatterPlot(data, idx1, idx2): 
    import matplotlib.pyplot as plt 
    fig = plt.figure() 
    ax = fig.add_subplot(111) 
    
    # X-axis data, Y-axis data, Size for each sample, Color for each sample 
    ax.scatter(data[:,idx1], data[:,idx2]) 
    plt.show()

def simple_reconstruction_3d_pytest(tim, lvl_str, use_guess):
	lvl=int(lvl_str)
	ber, gp=netcdf_utis.load_cube('/bm/gdata/scollis/cube_data/20060122_'+tim+'_ver_hr_big.nc')
	print gp['levs'][lvl]
	Re=6371.0*1000.0
	rad_at_radar=Re*sin(pi/2.0 -abs(gp['zero_loc'][0]*pi/180.0))#ax_radius(float(lat_cpol), units='degrees')
	lons=gp['zero_loc'][1]+360.0*gp['xar']/(rad_at_radar*2.0*pi)
	lats=gp['zero_loc'][0] + 360.0*gp['yar']/(Re*2.0*pi)	
	ber_loc=[-12.457, 130.925]
	gp_loc=	 [-12.2492,  131.0444]
	if use_guess=='none':
		igu=ones(ber['CZ'].shape, dtype=float)*0.0
		igv=ones(ber['CZ'].shape, dtype=float)*0.0
	else:
		ber_ig, gp_ig=netcdf_utis.load_cube(use_guess)
		print gp_ig.keys()
		igu=gp_ig['u_array']
		igv=gp_ig['v_array']
	mywts=ones(ber['CZ'].shape, dtype=float)
	angs=array(propigation.make_lobe_grid(ber_loc, gp_loc, lats,lons))
	wts_ang=zeros(gp['CZ'][:,:,0].shape, dtype=float)
	for i in range(angs.shape[0]):
			for j in range(angs.shape[1]):
				if (angs[i,j] < 150.0) and (angs[i,j] > 30.0): wts_ang[i,j]=1.0
	for lvl_num in range(len(gp['levs'])):
		#create a weighting grid
		mask_reflect=10.0#dBZ	
		mask=(gp['CZ'][:,:,lvl_num]/mask_reflect).round().clip(min=0., max=1.0) 
		mask_vel_ber=(ber['VR'][:,:,lvl_num]+100.).clip(min=0., max=1.)
		mywts[:,:,lvl_num]=mask*mask_vel_ber*wts_ang	
	f=0.0
	gv_u=zeros(ber['CZ'].shape, dtype=float)
	gv_v=zeros(ber['CZ'].shape, dtype=float)
	wts=mask*mask_vel_ber*wts_ang
	gu,gv,f= grad_conj_solver_plus_plus.meas_cost(gv_u, gv_v, f, igu, igv, ber['i_comp'], ber['j_comp'], gp['i_comp'], gp['j_comp'],  ber['VR'], gp['VR'], mywts)
	print "Mean U gradient", gu.mean(), "gv mean", gv.mean(), "F ", f
	for i in range(len(gp['levs'])):
		print "U,V ", (igu[:,:,i]).sum()/mywts[:,:,i].sum(), (igv[:,:,i]).sum()/mywts[:,:,i].sum()
		#gv_u,gv_v,cost = vel_2d_cost(gv_u,gv_v,cost,u_array,v_array,i_cmpt_r1,j_cmpt_r1,i_cmpt_r2,j_cmpt_r2,vr1,vr2,weights,nx=shape(gv_u,0),ny=shape(gv_u,1))
      		print gracon_vel2d.vel_2d_cost(gv_u[:,:,i]*0.0,gv_v[:,:,i]*0.0,0.0,igu[:,:,i],igv[:,:,i],ber['i_comp'][:,:,i], ber['j_comp'][:,:,i], gp['i_comp'][:,:,i], gp['j_comp'][:,:,i],  ber['VR'][:,:,i], gp['VR'][:,:,i],mywts[:,:,i])[2]
	gv_u,gv_v,f,u_array,v_array = grad_conj_solver_plus_plus.gracon_vel2d_3d( gv_u, gv_v, f, igu, igv, ber['i_comp'], ber['j_comp'], gp['i_comp'], gp['j_comp'], ber['VR'], gp['VR'], mywts)#, nx=nx, ny=ny, nz=nz)
	gp.update({'u_array': u_array, 'v_array':v_array})
	netcdf_utis.save_data_cube(ber, gp, '/bm/gdata/scollis/cube_data/20060122_'+tim+'_winds_ver1.nc', gp['zero_loc'])
	plotit=True
	if plotit:
		for lvl in range(len(gp['levs'])):
			print lvl
			f=figure()
			mapobj=pres.generate_darwin_plot(box=[130.8, 131.2, -12.4, -12.0])
			diff=gp['VR']-(u_array*gp['i_comp']+ v_array*gp['j_comp'])
			gp.update({'diff':diff})
			pres.reconstruction_plot(mapobj, lats, lons, gp, lvl, 'diff',u_array[:,:,lvl],v_array[:,:,lvl], angs, mywts[:,:,lvl])
			#pres.quiver_contour_winds(mapobj, lats, lons, (wts*u_array).clip(min=-50, max=50),(wts*v_array).clip(min=-50, max=50))
			t1='Gunn Point CAPPI (dBZ) and reconstructed winds (m/s) at %(lev)05dm \n 22/01/06 ' %{'lev':gp['levs'][lvl]}
			title(t1+tim) 
			ff=os.getenv('HOME')+'/bom_mds/output/recons_22012006/real_%(lev)05d_' %{'lev':gp['levs'][lvl]}
			savefig(ff+tim+'_2d_3d.png')
			close(f)	

def test_uniform_winds():
	f=figure()
	lats=linspace(-13.0, -11.5, 20)
	lons=linspace(130., 132., 20)
	u,v=simul_winds.unif_wind(lats, lons, 10.0, 45.0)
	mapobj=pres.generate_darwin_plot()
	pres.quiver_winds(mapobj, lats, lons, u,v)
	savefig(os.getenv('HOME')+'/bom_mds/output/test_uniform.png')
	close(f)