def plot_example_spectrograms(example,rate):
    """
    This function creates a figure with spectrogram sublpots to of the four
    sleep examples. (Recall row 0 is REM, rows 1-3 are NREM stages 1,
    2 and 3/4)
    """
    sleep_stages = ['REM sleep', 'Stage 1 NREM sleep', 'Stage 2 NREM sleep', 'Stage 3 and 4 NREM sleep'];    
    
    plt.figure()
    
    ###YOUR CODE HERE
    for i in range( len(example[:,0]) ):
       
       # plot every sleep stage in a separate plot
       plt.subplot(2,2,i+1)
       
       # plot spectogram
       plt.specgram(example[i, :],NFFT=512,Fs=rate)   
       
       # add legend
       plt.xlabel('Time (Seconds)')
       plt.ylabel('Frequency (Hz)')
       plt.title( 'Spectogram ' + sleep_stages[i] )
       
       plt.ylim(0,60)
       plt.xlim(0,290)
    return

    def plot(self, bit_stream):
        if self.previous_bit_stream != bit_stream.to_list():
            self.previous_bit_stream = bit_stream

            x = []
            y = []
            bit = None

            for bit_time in bit_stream.to_list():
                if bit is None:
                    x.append(bit_time)
                    y.append(0)
                    bit = 0
                elif bit == 0:
                    x.extend([bit_time, bit_time])
                    y.extend([0, 1])
                    bit = 1
                elif bit == 1:
                    x.extend([bit_time, bit_time])
                    y.extend([1, 0])
                    bit = 0

            plt.clf()
            plt.plot(x, y)
            plt.xlim([0, 10000])
            plt.ylim([-0.1, 1.1])
            plt.show()
            plt.pause(0.005)

    def test1():
        x = [0.5]*3
        xbounds = [(-5, 5) for y in x]


        GA = GenAlg(fitcalc1, x, xbounds, popMult=100, bitsPerGene=9, mutation=(1./9.), crossover=0.65, crossN=2, direction='min', maxGens=60, hammingDist=False)
        results = GA.run()
        print "*** DONE ***"
        #print results
        plt.ioff()
        #generate pareto frontier numerically
        x1_ = np.arange(-5., 0., 0.05)
        x2_ = np.arange(-5., 0., 0.05)
        x3_ = np.arange(-5., 0., 0.05)

        pfn = []
        for x1 in x1_:
            for x2 in x2_:
                for x3 in x3_:
                    pfn.append(fitcalc1([x1,x2,x3]))

        pfn.sort(key=lambda x:x[0])
        
        plt.figure()
        i = 0
        for x in results:
            plt.scatter(x[1][0], x[1][1], 20, c='r')

        plt.scatter([x[0] for x in pfn], [x[1] for x in pfn], 1.0, c='b', alpha=0.1)
        plt.xlim([-20,-1])
        plt.ylim([-12, 2])
        plt.draw()

def fit_plot_unlabeled_data(unlabeled_data_x, labeled_data_x, labeled_data_y, fit_order, data_type, other_data_list, other_data_name):
    output = open('predictions.csv','wb')
    coeffs = np.polyfit(labeled_data_x, labeled_data_y, fit_order) #does poly git to nth deg on labeled data
    fit_eq = np.poly1d(coeffs) #Eqn from fit
    predicted_y = fit_eq(unlabeled_data_x)
    i = 0
    writer = csv.writer(output,delimiter=',')
    header = [str(data_type),str(other_data_name),'Predicted_Num_Inc']
    writer.writerow(header)
    while i < len(predicted_y):
        output_data = [unlabeled_data_x[i],other_data_list[i],predicted_y[i]]
        writer.writerow(output_data)
        print 'For '+str(data_type)+' of: '+str(unlabeled_data_x[i])+', Predicted Number of Incidents is: '+str(predicted_y[i])
        i = i + 1
    plt.scatter(unlabeled_data_x, predicted_y, color='blue', label='Predicted Number of Incidents')
    fit_line_x = np.arange(min(unlabeled_data_x), max(unlabeled_data_x), 1)
    plt.plot(fit_line_x, fit_eq(fit_line_x), color='red',linestyle='dashed',label=' Order '+str(fit_order)+' Polynomial Fit')
#____Use below line to plot actual data also!! 
    #plt.scatter(labeled_data_x, labeled_data_y, color='green', label='Actual Incident Report Data')
    plt.title('Predicted Number of 311 Incidents by '+str(data_type))
    plt.xlabel(str(data_type))
    plt.ylabel('Number of 311 Incidents')
    plt.grid()
    plt.xlim([min(unlabeled_data_x)-1500, max(unlabeled_data_x)+1500])
    plt.legend(loc='upper left')
    plt.show()

def plot_q(model='cem', r_min=0.0, r_max=6371.0, dr=1.0):
    """
    Plot a radiallysymmetric Q model.

    plot_q(model='cem', r_min=0.0, r_max=6371.0, dr=1.0):

    r_min=minimum radius [km], r_max=maximum radius [km], dr=radius
    increment [km]

    Currently available models (model): cem, prem, ql6
    """
    import matplotlib.pylab as plt

    r = np.arange(r_min, r_max + dr, dr)
    q = np.zeros(len(r))

    for k in range(len(r)):

        if model == 'cem':
            q[k] = q_cem(r[k])
        elif model == 'ql6':
            q[k] = q_ql6(r[k])
        elif model == 'prem':
            q[k] = q_prem(r[k])

    plt.plot(r, q, 'k')
    plt.xlim((0.0, r_max))
    plt.xlabel('radius [km]')
    plt.ylabel('Q')
    plt.show()

def plot_average(filenames, save_plot=True, show_plot=False, dpi=100):

    ''' Plot Signal average from a list of averaged files. '''

    fname = get_files_from_list(filenames)

    # plot averages
    pl.ioff()  # switch off (interactive) plot visualisation
    factor = 1e15
    for fnavg in fname:
        name = fnavg[0:len(fnavg) - 4]
        basename = os.path.splitext(os.path.basename(name))[0]
        print fnavg
        # mne.read_evokeds provides a list or a single evoked based on condition.
        # here we assume only one evoked is returned (requires further handling)
        avg = mne.read_evokeds(fnavg)[0]
        ymin, ymax = avg.data.min(), avg.data.max()
        ymin *= factor * 1.1
        ymax *= factor * 1.1
        fig = pl.figure(basename, figsize=(10, 8), dpi=100)
        pl.clf()
        pl.ylim([ymin, ymax])
        pl.xlim([avg.times.min(), avg.times.max()])
        pl.plot(avg.times, avg.data.T * factor, color='black')
        pl.title(basename)

        # save figure
        fnfig = os.path.splitext(fnavg)[0] + '.png'
        pl.savefig(fnfig, dpi=dpi)

    pl.ion()  # switch on (interactive) plot visualisation

def draw_lineplot(x, y, title="title", xlab="x", ylab="y", odir="", xlim=None, ylim=None, outfmt='eps'):

  if len(x) == 0 or len(y) == 0:
    return;
  #fi

  plt.cla();
  plt.plot(x, y, marker='x');
  plt.xlabel(xlab);
  plt.ylabel(ylab);
  plt.title(title);

  if xlim == None:
    xmin = min(x);
    xmax = max(x);
    xlim = [xmin, xmax];
  #fi

  if ylim == None:
    ymin = min(y);
    ymax = max(y);
    ylim = [ymin, ymax];
  #fi

  plt.xlim(xlim);
  plt.ylim(ylim);

  plt.savefig('%s%s.%s' % (odir + ('/' if odir else ""), '_'.join(title.split(None)), outfmt), format=outfmt);

  return '%s.%s' % ('_'.join(title.split(None)), outfmt), title;

def plot_prediction_accuracy(x, y):
    plt.scatter(x, y, c='g', alpha=0.5)
    plt.title('Logistic Regression')
    plt.xlabel('r')
    plt.ylabel('Prediction Accuracy')
    plt.xlim(0,200)
    plt.show()

def plot_eye(Nodes,axes = None):
    """
    
    Create a movie of eye growth. To be used with EyeGrowthFunc

    :param Nodes: structure containing nodes
    :type Nodes: struct
    :param INTSTEP: time step used for integration
    :type INTSTEP: int

    :returns: plot handle for Node plot.  Used to update during for loop.

    .. note::
       Called in EyeGrowthFunc
    """
    
    
    #set plotting parameters:
    if axes == None:
        fig = plt.figure(figsize=(10, 8))
        axes = fig.add_subplot(111,aspect='equal')
        plt.xlim([-13, 13])
        plt.ylim([-13, 13])
        
    axes.plot(np.r_[ Nodes['x'][0],Nodes['x'][0,0] ] * Nodes['radius'], 
            np.r_[ Nodes['y'][0], Nodes['y'][0,0] ] * Nodes['radius'], 
            '-ok', markerfacecolor = 'k',linewidth = 4, markersize = 10)
                
    axes = pf.TufteAxis(axes,['left','bottom'])
    #axes.set_axis_bgcolor('w')


    return axes

def plot_tuning_curves(direction_rates, title):
    """
    This function takes the x-values and the y-values  in units of spikes/s 
    (found in the two columns of direction_rates) and plots a histogram and 
    polar representation of the tuning curve. It adds the given title.
    """
    x = direction_rates[:,0]
    y = direction_rates[:,1]
    plt.figure()
    plt.subplot(2,2,1)
    plt.bar(x,y,width=45,align='center')
    plt.xlim(-22.5,337.5)
    plt.xticks(x)
    plt.xlabel('Direction of Motion (degrees)')
    plt.ylabel('Firing Rate (spikes/s)')
    plt.title(title)   
        
        
    
    plt.subplot(2,2,2,polar=True)
    r = np.append(y,y[0])
    theta = np.deg2rad(np.append(x, x[0]))
    plt.polar(theta,r,label='Firing Rate (spikes/s)')
    plt.legend(loc=8)
    plt.title(title)

def plot_lift_data(lift_data, with_ellipses=True):
    np.random.seed(42113)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    alpha = [l['fit']['alpha'] for l in lift_data.values()]
    alpha_error = [l['fit']['alpha_error'] for l in lift_data.values()]
    beta = [l['fit']['beta'] for l in lift_data.values()]
    beta_error = [l['fit']['beta_error'] for l in lift_data.values()]
    message_class = lift_data.keys()

    num = len(beta)
    beta_jitter = np.random.randn(num)
    np.random.seed(None)
    beta = np.array(beta) + beta_jitter*0.0

    ax.plot(beta, alpha, color='red', linestyle='', marker='o', markersize=10)
    if not with_ellipses:
        ax.errorbar(beta, alpha, xerr=beta_error, yerr=alpha_error, linestyle='')
    else:
        for x, y, xerr, yerr, in zip(beta, alpha, beta_error, alpha_error):
            width = 2*xerr
            height = 2*yerr
            ellipse = patches.Ellipse((x, y), width, height,
                                      angle=0.0, linewidth=2,
                                      fill=True, alpha=0.15, color='gray')
            ax.add_patch(ellipse)

    for a, b, c in zip(alpha, beta, message_class):
        ax.annotate(c, xy=(b, a), xytext=(b+2, a+.01), fontsize=17)
    plt.xlim(0, max(beta)+30)
    plt.ylim(0, 0.9)
    plt.xlabel('Duration (days)')
    plt.ylabel('Initial Lift')
    plt.show()

def unteraufgabe_g():
    # Sampling punkte
    x = np.linspace(0.0,1.0,1000)
    N = np.arange(2,16)

    LU = np.ones_like(N,dtype=np.floating)
    LT = np.ones_like(N,dtype=np.floating)

    # Approximiere Lebesgue-Konstante
    for i,n in enumerate(N):
        ################################################################
        #
        # xU = np.linspace(0.0,1.0,n)
        #
        # LU[i] = ...
        #
        # j = np.arange(n+1)
        # xT = 0.5*(np.cos((2.0*j+1.0)/(2.0*(n+1.0))*np.pi) + 1.0)
        #
        # LT[i] = ...
        #
        ################################################################
        continue

    # Plot
    plt.figure()
    plt.semilogy(N,LU,"-ob",label=r"Aequidistante Punkte")
    plt.semilogy(N,LT,"-og",label=r"Chebyshev Punkte")
    plt.grid(True)
    plt.xlim(N.min(),N.max())
    plt.xlabel(r"$n$")
    plt.ylabel(r"$\Lambda^{(n)}$")
    plt.legend(loc="upper left")
    plt.savefig("lebesgue.eps")

def plot(all_models):

    import matplotlib.pylab as plt
    import numpy.random
    plt.close("all")
    plt.figure()
    plt.subplot(211)
    alt = np.arange(0., 500., 2.)
    sza = 0.

    for m in all_models:
        d = m(alt, sza)
        plt.plot(ne_to_fp(d)/1E6, alt,lw=2)
        # plt.plot(m(alt, sza),alt,lw=2)

    plt.ylim(0., 400.)
    plt.ylabel('Altitude / km')
    # plt.xlabel(r'$n_e / cm^{-3}$')
    plt.xlabel(r'$f / MHz$')
    plt.subplot(212)
    for m in all_models:
        delay, freq = m.ais_response()
        plt.plot(freq/1E6, delay*1E3, lw=2.)

    plt.hlines(-2 * np.amax(alt) / speed_of_light_kms * 1E3, *plt.xlim(), linestyle='dashed')
    # plt.vlines(ne_to_fp(1E5)/1E6, *plt.ylim())
    # plt.hlines(  -(500-150) * 2 / speed_of_light_kms * 1E3, *plt.xlim())
    plt.ylim(-10,0)
    plt.ylabel('Delay / ms')
    plt.xlim(0, 7)
    plt.xlabel('f / MHz')
    plt.show()

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 plot(self,file):
        cds = CaseDataset(file, 'bson')
        data = cds.data.driver('driver').by_variable().fetch()
        cds2 = CaseDataset('../output/therm_mc_20141110173851.bson', 'bson')
        data2 = cds2.data.driver('driver').by_variable().fetch()
        
        #temp
        temp_boundary_k = data['hyperloop.temp_boundary']
        temp_boundary_k.extend(data2['hyperloop.temp_boundary'])
        temp_boundary = [((x-273.15)*1.8 + 32) for x in temp_boundary_k]
        #histogram
        n, bins, patches = plt.hist(temp_boundary, 100, normed=1, histtype='stepfilled')
        plt.setp(patches, 'facecolor', 'b', 'alpha', 0.75)

        #stats
        mean = np.average(temp_boundary)
        std = np.std(temp_boundary)
        percentile = np.percentile(temp_boundary,99.5)
        print "mean: ", mean, " std: ", std, " 99.5percentile: ", percentile
        x = np.linspace(50,170,150)
        plt.plot(x,mlab.normpdf(x,mean,std), color='black', lw=2)
        plt.xlim([60,160])
        plt.ylabel('Probability', fontsize=18)
        plt.xlabel(u'Equilibrium Temperature, \N{DEGREE SIGN}F', fontsize=18)
        #plt.show()
        plt.tight_layout()
        plt.savefig('../output/histo.pdf', dpi=300)

def plot_behavior_count(agent_id, behavior_count):
    """Create a plot for the behavior count.

    Args:
        agent_id: The identifier of the agent.
        behavior_count: The count of the agent behaviors.
    """
    fig = plt.figure()
    ax = fig.add_subplot(111)
    plt.xlabel(u'Número de jogos')
    plt.ylabel(u'Probabilidades de selecionar comportamento')
    plt.title(u'Probabilidades do agente %d selecionar comportamento'
              % agent_id)

    plt.xlim([0, 115])
    plt.ylim([-0.1, 1.1])

    data = np.array([b for b in behavior_count.values()])
    prob = data/np.sum(data, axis=0)

    for i, behavior in enumerate(behavior_count):
        coeff = calculate_regression_coefficients(prob[i], degree=4)
        regression = [calculate_regression_y(x, coeff)
                      for x in range(len(prob[i]))]
        ax.plot(regression, label=T[behavior],
                c=COLOR_TABLE[COLOR_LIST[i]], linewidth=2.0)
        ax.scatter(range(len(prob[i])), prob[i], c=COLOR_TABLE[COLOR_LIST[i]],
                   alpha=1)

    ax.legend()

def plot_fullstack( binning = np.linspace(0,10,1), myquery='', plotvar = default_plot_variable, \
                    scalefactor = 1., user_ylim = None):

    fig = plt.figure(figsize=(10,6))
    plt.grid(True)
    lasthist = 0
    myhistos = gen_histos(binning=binning,myquery=myquery,plotvar=plotvar,scalefactor=scalefactor)
    for key, (hist, bins) in myhistos.iteritems():

      plt.bar(bins[:-1],hist,
              width=bins[1]-bins[0],
              color=colors[key],
              bottom = lasthist,
              edgecolor = 'k',
              label='%s: %d Events'%(labels[key],sum(hist)))
      lasthist += hist
     

    plt.title('CCSingleE Stacked Backgrounds',fontsize=25)
    plt.ylabel('Events',fontsize=20)
    if plotvar == '_e_nuReco' or plotvar == '_e_nuReco_better':
        xstring = 'Reconstructed Neutrino Energy [GeV]' 
    elif plotvar == '_e_CCQE':
        xstring = 'CCQE Energy [GeV]'
    else:
        xstring = plotvar
    plt.xlabel(xstring,fontsize=20)
    plt.legend()
    plt.xticks(list(plt.xticks()[0]) + [binning[0]])
    plt.xlim([binning[0],binning[-1]])

def plot_prob_effector(sens, fpr, xmax=1, baserate=0.1):
    """Plots a line graph of P(effector|positive test) against
    the baserate of effectors in the input set to the classifier.
        
    The baserate argument draws an annotation arrow
    indicating P(pos|+ve) at that baserate
    """
    assert 0.1 <= xmax <= 1, "Max x axis value must be in range [0,1]"
    assert 0.01 <= baserate <= 1, "Baserate annotation must be in range [0,1]"
    baserates = pylab.arange(0, 1.05, xmax * 0.005)  
    probs = [p_correct_given_pos(sens, fpr, b) for b in baserates]
    pylab.plot(baserates, probs, 'r')
    pylab.title("P(eff|pos) vs baserate; sens: %.2f, fpr: %.2f" % (sens, fpr))
    pylab.ylabel("P(effector|positive)")
    pylab.xlabel("effector baserate")
    pylab.xlim(0, xmax)
    pylab.ylim(0, 1)
    # Add annotation arrow
    xpos, ypos = (baserate, p_correct_given_pos(sens, fpr, baserate))
    if baserate < xmax:
        if xpos > 0.7 * xmax:
            xtextpos = 0.05 * xmax
        else:
            xtextpos = xpos + (xmax-xpos)/5.
        if ypos > 0.5:
            ytextpos = ypos - 0.05
        else:
            ytextpos = ypos + 0.05
        pylab.annotate('baserate: %.2f, P(pos|+ve): %.3f' % (xpos, ypos), 
                       xy=(xpos, ypos), 
                       xytext=(xtextpos, ytextpos),
                       arrowprops=dict(facecolor='black', shrink=0.05))
    else:
        pylab.text(0.05 * xmax, 0.95, 'baserate: %.2f, P(pos|+ve): %.3f' %
                   (xpos, ypos))

 def plot(x,y,field,filename,c=200):
     plt.figure()
     # define grid.
     xi = np.linspace(min(x),max(x),100)
     yi = np.linspace(min(y),max(y),100)
     # grid the data.
     si_lin = griddata((x, y), field, (xi[None,:], yi[:,None]), method='linear')
     si_cub = griddata((x, y), field, (xi[None,:], yi[:,None]), method='linear')
     print np.min(field)
     print np.max(field)
     plt.subplot(211)
     # contour the gridded data, plotting dots at the randomly spaced data points.
     CS = plt.contour(xi,yi,si_lin,c,linewidths=0.5,colors='k')
     CS = plt.contourf(xi,yi,si_lin,c,cmap=plt.cm.jet)
     plt.colorbar() # draw colorbar
     # plot data points.
     #    plt.scatter(x,y,marker='o',c='b',s=5)
     plt.xlim(min(x),max(x))
     plt.ylim(min(y),max(y))
     plt.title('Lineaarinen interpolointi')
     #plt.tight_layout()
     plt.subplot(212)
     # contour the gridded data, plotting dots at the randomly spaced data points.
     CS = plt.contour(xi,yi,si_cub,c,linewidths=0.5,colors='k')
     CS = plt.contourf(xi,yi,si_cub,c,cmap=plt.cm.jet)
     plt.colorbar() # draw colorbar
     # plot data points.
     #    plt.scatter(x,y,marker='o',c='b',s=5)
     plt.xlim(min(x),max(x))
     plt.ylim(min(y),max(y))
     plt.title('Kuubinen interpolointi')
     plt.savefig(filename)

def plotFirstTacROC(dataset):
    import matplotlib.pylab as plt
    from os.path import join
    from src.utils import PROJECT_DIR

    plt.figure(figsize=(6, 6))
    time_sampler = TimeSerieSampler(n_time_points=12)
    evaluator = Evaluator()
    time_series_idx = 0
    methods = {
        "cross_correlation": "Cross corr.   ",
        "kendall": "Kendall        ",
        "symbol_mutual": "Symbol MI    ",
        "symbol_similarity": "Symbol sim.",
    }
    for method in methods:
        print method
        predictor = SingleSeriesPredictor(good_methods[method], time_sampler)
        prediction = predictor.predictAllInstancesCombined(dataset, time_series_idx)
        roc_auc, fpr, tpr = evaluator.evaluate(prediction)
        plt.plot(fpr, tpr, label=methods[method] + " (auc = %0.3f)" % roc_auc)
    plt.legend(loc="lower right")
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.05])
    plt.xlabel("False Positive Rate")
    plt.ylabel("True Positive Rate")
    plt.grid()
    plt.savefig(join(PROJECT_DIR, "output", "firstTACROC.pdf"))