def plot_ex(ion, line):
    ''' Plot the collisions cross-section with electrons
    '''
    global cst, iplot
    iplot += 1
    ofname = format(iplot, '0>4')+'_ex.pdf'
    T = 5000.
    x0  = cst / line.lbd * Cst.Q / Cst.K / T
    x_vec_log = pl.linspace(-3, 4, 100)
    x_vec = 10**x_vec_log
    q0 = [sigma_weak_seaton(i, ion, line, T) for i in x_vec]
    q1 = [sigma_strong_seaton(i, ion, line, T) for i in x_vec]
    q = [sigma_seaton(i, ion, line, T) for i in x_vec]
    ll = line.lower
    ul = line.upper
    #title='Na I: '+str(ll.cfg)+' '+str(ll.term)+' [g='+str(ll.g)+'] <=> '+str(ul.cfg)+' '+str(ul.term)+' [g='+str(ul.g)+'], Ro = '+format(orbital_radius(ion, line.lower), '4.2f')+' a$_0$, f = '+str(line.f)
    title = SPE+' '+DEG+': '+str(ll.cfg)+' '+str(ll.term)+' [g='+str(ll.g)+'] '+str(ul.cfg)+' '+str(ul.term)+' [g='+str(ul.g)+'] Ro = '+str(orbital_radius(ion, line.lower))+' f = '+str(line.f)
    pl.figure(figsize=(12, 6), dpi=100, facecolor='w', edgecolor='k')
    pl.plot(x_vec*Cst.K*T/Cst.Q, q0, 'k', lw=1, label='Weak coupling')
    pl.plot(x_vec*Cst.K*T/Cst.Q, q1, 'k', lw=2, label='Strong coupling')
    pl.plot(x_vec*Cst.K*T/Cst.Q, q, 'r+', label='IPM (Seaton 1962)')
    pl.semilogx()
    #pl.semilogy()
    pl.xlim([1e-2, 1e4])
    pl.xlabel('E [eV]')
    pl.ylabel('Cross-section [$\pi a_0^2$]')
    pl.legend(frameon=False, loc=0)
    pl.title(title)
    #bbox_inches='tight'
    pl.savefig(ofname, dpi=100, format='pdf', orientation='landscape', papertype='a4')
    pl.close()

def plot():
    for time_filename, size_filename, f_label, f_style in lines_to_plot:
        time_file = open('parser_results/' + time_filename)
        size_file = open('parser_results/' + size_filename)
        times = {}
        for mark, _, _ in swarmsize_lines:
            times[mark] = []
        for line in time_file:
            time = float(line.split()[0])
            size = int(size_file.next().split()[0])
            for mark, _, _ in swarmsize_lines:
                if size <= mark:
                    times[mark].append(time)
                    break
        for mark, m_label, m_style in reversed(swarmsize_lines):
            cum_values = cdf.cdf(times[mark])            
            x = [cum_value[1] for cum_value in cum_values]
            y = [cum_value[0] for cum_value in cum_values]
            pylab.semilogx(x, y, f_style+m_style,
                           label=f_label+' '+m_label,
                           markevery=markevery)
    pylab.legend(loc='lower right')
    pylab.savefig(output_filename)
#    pylab.close()
    

    print 'Output saved to:', output_filename

def sums2nloop(model_input, images_input, iterations):
    """Run through a loop of calculating the summed S/N for a set of data.
    
    This is useful if you'd like to see how well we'd do as we keep adding
    more and more data together in order to attempt to extract a signal.
    
    """
    if images_input == "simdata":
        header = model_input.simdata()
    else:
        header = getheaderinfo(images_input)

    # iterations = header['numexp']

    snvec = numpy.zeros((iterations, 2))

    for i in range(0, iterations):
        snvec[i] = models2n(i + 1, model_input, images_input)

    print snvec

    # GRAB EXPOSURE TIMES
    timearr = header["tstart"][0:iterations]  # should be error bar from tstart to tstop

    # Plot the results
    pylab.semilogx(timearr, snvec[:, 1], linestyle="", marker="o")
    pylab.semilogx(timearr, snvec[:, 0], linestyle="", marker="x")
    pylab.xlim(50, 200000)
    pylab.xlabel("Time")
    pylab.ylabel("Summed S/N (not a lightcurve!)")
    plottitle = "Model: " + model_input.name + ", Data: " + images_input
    pylab.title(plottitle)

def filterresponse(b,a,fs=44100,scale='log',**kwargs):
    w, h = freqz(b,a)
    pl.subplot(2,1,1)
    pl.title('Digital filter frequency response')
    pl.plot(w/max(w)*fs/2, 20 * np.log10(abs(h)),**kwargs)
    pl.xscale(scale)
#    if scale=='log':
#        pl.semilogx(w/max(w)*fs/2, 20*np.log10(np.abs(h)), 'k')
#    else:
#        pl.plot(w/max(w)*fs/2, 20*np.log10(np.abs(h)), 'k')
        
    pl.ylabel('Gain (dB)')
    pl.xlabel('Frequency (rad/sample)')
    pl.axis('tight')    
    pl.grid()    
    

    pl.subplot(2,1,2)
    angles = np.unwrap(np.angle(h))
    if scale=='log':
        pl.semilogx(w/max(w)*fs/2, angles, **kwargs)
    else:
        pl.plot(w/max(w)*fs/2, angles, **kwargs)
        

    pl.ylabel('Angle (radians)')
    pl.grid()
    pl.axis('tight')
    pl.xlabel('Frequency (rad/sample)')

 def noise_data(self, chan=0, m=1, plot=False):
     NFFT = self.NFFT
     N = self.N
     lendata = self.lendata
     Fs = self.Fs
     pst = np.zeros(NFFT)
     cnt = 0
     while cnt < m:
         data, addr = self.r.get_data_udp(N*lendata, demod=True)
         if self.use_r2 and ((not np.all(np.diff(addr)==2**21/N)) or data.shape[0]!=256*lendata):
             print "bad"
         else:
             ps, f = mlab.psd(data[:,chan], NFFT=NFFT, Fs=Fs/self.r.nfft)
             pst += ps
             cnt += 1
     pst /= cnt
     pst = 10.*np.log10(pst)
     if plot:
         ind = f > 0
         pl.semilogx(f[ind], pst[ind])
         pl.xlabel('Hz')
         pl.ylabel('dB/Hz')
         pl.title('chan data psd')
         pl.grid()
         pl.show()
     return f, pst

def test_cross_validation():
    np.random.seed(1)
    n_features = 3
    func = get_test_model(n_features=n_features)
    dataset_train, labels_train = get_random_dataset(func, 
        n_datapoints=5e3, n_features=n_features, noise_level=1.)
    alphas = 10**np.linspace(-3, 2, 10)

    # Fit scikit lasso
    lasso_scikit = lasso.SKLasso(alpha=0.001, max_iter=1000)
    lasso_scikit.fit_CV(dataset_train, labels_train[:,0], alphas=alphas,
                        n_folds=5)

    # Fit tf lasso
    gen_lasso = gl.GeneralizedLasso(alpha=0.001, max_iter=1000,
                                    link_function=None)
    gen_lasso.fit_CV(dataset_train, labels_train[:,0], alphas=alphas,
                     n_folds=5)

    pl.figure()
    pl.title('CV test. TF will have higher errors, since it is not exact')
    pl.semilogx(alphas, gen_lasso.alpha_mse, 'o-', label='tf')
    pl.semilogx(alphas, lasso_scikit.alpha_mse, '*-', label='scikit')
    pl.legend(loc='best')
    pl.xlabel('alpha')
    pl.ylabel('cost')
    pl.show()

def histo(denlist):

    for d in denlist:
        h,bin_edges = N.histogram(N.log10(d.flatten()),normed=True,range=(-2,2),bins=50)
        M.semilogx(10.**(bin_edges),h)

    M.show()

def plotResistance():
    
    dic = makeListResistanceSims(3)
    
    pl.close("all")
    
    colors = []
    
    for k,key in enumerate(np.sort(dic.keys())):
        pl.figure(k,(11,7))
        pl.title("starting property violation: %s"%key)
        
        
        low = 0
        high = 0
        
        if rootDir == '/Users/maithoma/work/compute/pgames_resistance/':
            range = np.arange(key,0.11,0.005)
            iter = 2000
             
        elif rootDir == '/Users/maithoma/work/compute/pgames_resistance2/':
            range = np.arange(0.043,0.062,0.002)
            iter = 4000
                
        color = np.linspace(0.7,0.3,len(range))
    
        for i,s in enumerate(range):
            for j,jx in enumerate(dic[key]):
                
                if j==0:
                    label = str(s)
                elif j>0:
                    label = '_nolegend_'
                    
                
                filename = 'iter_%s_l_100_h_100_d_0.500_cl_0.500_ns_4_il_1.000_q_0.000_M_5_m_0.000_s_%.4f_%s.csv'%(iter,s,jx)
                print filename     
                #if s < 0.04 or s > 0.065:
                #    continue
                
                try:
                    x,y = coopLevel(filename)
                    #print key,i,jx,s,y[-1],filename
                    if float(y[-1]) < 0.7:
                        pl.semilogx(x,y,alpha=1,lw= 1.,label=label,color=[color[i],color[i],1])
                        #pl.semilogx(x,y,alpha=1,lw= 1.,label=label,color='c')
                        low +=1
                    elif float(y[-1]) > 0.7:
                        pl.semilogx(x,y,alpha=1,lw= 1.,label=label,color=[color[i],1,color[i]])
                        #pl.semilogx(x,y,alpha=1,lw= 1.,label=label,color='m')
                        high +=1
                except:
                    continue

        #pl.text(100000,0.85,"high:%s \n low:%s"%(high,low))       
        pl.legend(loc=0)
        pl.xlabel("Iterations (log10)")
        pl.ylabel("Cooperation Level")
        pl.xlim(xmax=5*10**8)
        pl.ylim(0.,1)

def clust(Nclust=25,loadfile='/media/HD1_/Documents/AG3C/Results/formated_data.p',writefile='/media/HD1_/Documents/AG3C/Results/cluster_kmeans.p',plots=True):

    
    #load the formated data
    (avout,stdout,av_reference)=pickle.load(open(loadfile,'rb'))
    
    TFav=[] 
    TFref=[] 
    TFstd=[] 
    
    #do k-means clustering
    (o,n)=kmeans2(np.array(avout),Nclust)
    
    #save the output
    pickle.dump((o,n),open(writefile,'wb'))

    if plots==True:
        #make the plots
        t=[3.,4.,5.,6.,8.,24.,48,7*24.,14*24.]
        plt.figure()
        for i in range(Nclust):
            
            idx=np.where(n==i)[0]
           
            plt.subplot(np.ceil(np.sqrt(Nclust)),np.ceil(np.sqrt(Nclust)),i+1)
            plt.hold(True)
            for j in range(len(idx)):
                plt.semilogx(t,avout[idx[j]])
                plt.ylim([0,1])
                
        plt.show()

def chunked_timing(X, Y, axis=1, metric="euclidean", **kwargs):
    sizes = [20, 50, 100,
             200, 500, 1000,
             2000, 5000, 10000,
             20000, 50000, 100000,
             200000]

    t0 = time.time()
    original(X, Y, axis=axis, metric=metric, **kwargs)
    t1 = time.time()
    original_timing = t1 - t0

    chunked_timings = []

    for batch_size in sizes:
        print("batch_size: %d" % batch_size)
        t0 = time.time()
        chunked(X, Y, axis=axis, metric=metric, batch_size=batch_size,
                **kwargs)
        t1 = time.time()
        chunked_timings.append(t1 - t0)

    import pylab as pl
    pl.semilogx(sizes, chunked_timings, '-+', label="chunked")
    pl.hlines(original_timing, sizes[0], sizes[-1],
              color='k', label="original")
    pl.grid()
    pl.xlabel("batch size")
    pl.ylabel("execution time (wall clock)")
    pl.title("%s %d / %d (axis %d)" % (
        str(metric), X.shape[0], Y.shape[0], axis))
    pl.legend()
    pl.savefig("%s_%d_%d_%d" % (str(metric), X.shape[0], Y.shape[0], axis))
    pl.show()

def extrapolate(n, y, tolerance=1e-15, plot=False, call_show=True):
    "Extrapolate functional value Y from sequence of values (n, y)."

    # Make sure we have NumPy arrays
    n = array(n)
    y = array(y)

    # Create initial "bound"
    Y0 = 0.99*y[-1]
    Y1 = 1.01*y[-1]

    # Compute initial interior points
    phi = (sqrt(5.0) + 1.0) / 2.0
    Y2 = Y1 - (Y1 - Y0) / phi
    Y3 = Y0 + (Y1 - Y0) / phi

    # Compute initial values
    F0, e, nn, ee, yy = _eval(n, y, Y0)
    F1, e, nn, ee, yy = _eval(n, y, Y1)
    F2, e, nn, ee, yy = _eval(n, y, Y2)
    F3, e, nn, ee, yy = _eval(n, y, Y3)

    # Solve using direct search (golden ratio fraction)
    while Y1 - Y0 > tolerance:

        if F2 < F3:
            Y1, F1 = Y3, F3
            Y3, F3 = Y2, F2
            Y2 = Y1 - (Y1 - Y0) / phi
            F2, e, nn, ee, yy = _eval(n, y, Y2)
        else:
            Y0, F0 = Y2, F2
            Y2, F2 = Y3, F3
            Y3 = Y0 + (Y1 - Y0) / phi
            F3, e, nn, ee, yy = _eval(n, y, Y3)

        print Y0, Y1

    # Compute reference value
    Y = 0.5*(Y0 + Y1)

    # Print results
    print
    print "Reference value:", Y

    # Plot result
    if plot:
        pylab.figure()
        pylab.subplot(2, 1, 1)
        pylab.title("Reference value: %g" % Y)
        pylab.semilogx(n, y, 'b-o')
        pylab.semilogx(nn, yy, 'g--')
        pylab.subplot(2, 1, 2)
        pylab.loglog(n, e, 'b-o')
        pylab.loglog(nn, ee, 'g--')
        pylab.grid(True)
        if call_show:
            pylab.show()

    return Y

def flipPlot(minExp, maxExp):
    """Assumes minExp and maxExp positive integers; minExp < maxExp
    Plots results of 2**minExp to 2**maxExp coin flips"""
    ratios = []
    diffs = []
    xAxis = []
    for exp in range(minExp, maxExp + 1):
        xAxis.append(2 ** exp)
    for numFlips in xAxis:
        numHeads = 0
        for n in range(numFlips):
            if random.random() < 0.5:
                numHeads += 1
        numTails = numFlips - numHeads
        ratios.append(numHeads / float(numTails))
        diffs.append(abs(numHeads - numTails))
            
    pylab.title('Difference Between Heads and Tails')
    pylab.xlabel('Number of Flips')
    pylab.ylabel('Abs(#Heads - #Tails)')
    pylab.rcParams['lines.markersize'] = 10
    pylab.semilogx()
    pylab.semilogy()
    pylab.plot(xAxis, diffs, 'bo')
    pylab.figure()
    pylab.title('Heads/Tails Ratios')
    pylab.xlabel('Number of Flips')
    pylab.ylabel('Heads/Tails')
    pylab.plot(xAxis, ratios)

def flipPlot(minExp,maxExp):
    '''假定minExp和maxExp是正整数，并且minExp＜maxExp，
       绘制出2**minExp到2**maxExp次抛硬币的结果'''
    ratios=[]
    diffs=[]
    xAxis=[]
    for exp in range(minExp,maxExp+1):
        xAxis.append(2**exp)
    for numFlips in xAxis:
        numHeads=0
        for n in range(numFlips):
            if random.random()<0.5:
                numHeads+=1
        numTails=numFlips-numHeads
        ratios.append(numHeads/float(numTails))
        diffs.append(abs(numHeads-numTails))
    pylab.title('Difference Between Heads and Tails')
    pylab.xlabel('Number of Flips')
    pylab.semilogx()
    pylab.semilogy()
    pylab.ylabel('Abs(#Heads-#Tails)')
    pylab.plot(xAxis,diffs,'bo')
    pylab.figure()
    pylab.title('Heads/Tails Ratios')
    pylab.xlabel('Number of Flips')
    pylab.semilogx()
    pylab.ylabel('#Heads/#Tails')
    pylab.plot(xAxis,ratios,'bo')

    def plot(self, marker='o', color='red', m=0, M=6000):
        """Plots the position / height of the peaks in the well

        .. plot::
            :include-source:

            from fragment_analyser import Line, fa_data
            l = Line(fa_data("alternate/peaktable.csv"))
            well = l.wells[0]
            well.plot()

        """
        import pylab
        if len(self.df) == 0:
            print("Nothing to plot (no peaks)")
            return
        x = self.df['Size (bp)'].astype(float).values
        y = self.df['RFU'].astype(float).values

        pylab.stem(x, y, marker=marker, color=color)
        pylab.semilogx()
        pylab.xlim([1, M])
        pylab.ylim([0, max(y)*1.2])
        pylab.grid(True)
        pylab.xlabel("size (bp)")
        pylab.ylabel("RFU")
        return x, y

def flipPlot(minExp,maxExp):
	ratios = []
	diffs = []
	xAxis = []

	for exp in range(minExp,maxExp+1):
		xAxis.append(2 ** exp)
	print "xAxis: ", xAxis
			
	for numFlips in xAxis:
		numHeads = 0
		for n in range(numFlips):
			if random.random() < 0.5:
				numHeads += 1
		numTails = numFlips - numHeads
		ratios.append(numHeads/float(numTails))
		diffs.append(abs(numHeads - numTails))

	pylab.figure()
	pylab.title('Difference Between Heads and Tails')
	pylab.xlabel('Number of Flips')
	pylab.ylabel('Abs(#Heads - #Tails')
	pylab.plot(xAxis, diffs, 'bo') #do not connect, show dot
	pylab.semilogx()
	pylab.semilogy()
	pylab.figure()
	pylab.plot(xAxis, ratios, 'bo') #do not connect, show dot
	pylab.title('Heads/Tails Ratios')
	pylab.xlabel('Number of Flips')
	pylab.ylabel('Heads/Tails')
	pylab.semilogx()

def plotBode(fAg, lowerFreq, higherFreq = None):
    """'  plot Bode diagram using matplotlib
        Default frequency width is 3 decades
    '"""
    import pylab as py
    #import pdb; pdb.set_trace()
    if ( higherFreq == None):
        rangeAt = 1000.0  # 3 decade
    else:
        assert higherFreq > lowerFreq
        rangeAt = float(higherFreq)/lowerFreq

    N = 128
    lstScannAngFreqAt = [2*sc.pi*1j*lowerFreq*sc.exp(
                                                sc.log(rangeAt)*float(k)/N)
                                for k in range(N)]

    t              = [        lowerFreq*sc.exp(sc.log(rangeAt)*float(k)/N)
                                for k in range(N)]

    py.subplot(211)
    py.loglog( t, [(abs(fAg(x))) for x in lstScannAngFreqAt]  )
    py.ylabel('gain')
    py.grid(True)

    py.subplot(212)
    py.semilogx( t, [sc.arctan2(fAg(zz).imag, fAg(zz).real)
                        for zz in lstScannAngFreqAt])
    py.ylabel('phase')
    py.grid(True)
    py.gca().xaxis.grid(True, which='minor')  # minor grid on too
    py.show()

 def plotResults(self, titlestr="", ylimits=[0.5,1.05], plotfunc = pl.semilogx, ylimitsB=[0,101],
                  legend_loc=3, show=True ):
     pl.figure(num=None, figsize=(15,5))
     xvals = range(1, (1+len(self.removed)) )
     #Two subplots. One the left is the test accuracy vs. iteration
     pl.subplot(1,2,1)
     plotfunc(xvals, self.test_acc_list, "b", label="Test Accuracy")    
     pl.hold(True)
     plotfunc(xvals, self.getRollingAvgTestAcc(window_size=10), "r", label="Test Acc (rolling avg)")
     plotfunc(xvals, self.getRollingAvgTrainAcc(window_size=10), "g--", label="Train Acc (rolling avg)")
     pl.ylim(ylimits)
     if titlestr == "":
         pl.title("Iterative Feature Removal")
     else:
         pl.title(titlestr)
     pl.ylabel("Test Accuracy")
     pl.xlabel("Iteration")
     pl.legend(loc=legend_loc) #3=lower left
     pl.hold(False)
     
     #second subplot. On the right is the number of features removed per iteration
     pl.subplot(1,2,2)
     Ns = [ len(lst) for lst in self.removed ]
     pl.semilogx(xvals, Ns, "bo", label="#Features per Iteration")
     pl.xlabel("Iteration")
     pl.ylabel("Number of Features Selected")
     pl.title("Number of Features Removed per Iteration")
     pl.ylim(ylimitsB)
     
     pl.subplots_adjust(left=0.05, bottom=0.15, right=0.95, top=0.90, wspace=0.20, hspace=0.20)
     if show: pl.show()

def flipPlot(minExp, maxExp, numTrials):
    meanRatios = []
    meanDiffs = []
    ratiosSDs =  []
    diffsSDs =  []
    xAxis = []
    for exp in range(minExp, maxExp + 1):
        xAxis.append(2**exp)
    for numFlips in xAxis:
        ratios = []
        diffs = []
        for t in range(numTrials):
            numHeads, numTails = runTrial(numFlips)
            ratios.append(numHeads/float(numTails))
            diffs.append(abs(numHeads - numTails))
        meanRatios.append(sum(ratios)/numTrials)
        meanDiffs.append(sum(diffs)/numTrials)
        ratiosSDs.append(stdDev(ratios))
        diffsSDs.append(stdDev(diffs))
    pylab.plot(xAxis, meanRatios, 'bo')
    pylab.title('Mean Heads/Tails Ratios ('
                + str(numTrials) + ' Trials)')
    pylab.xlabel('Number of Flips')
    pylab.ylabel('Mean Heads/Tails')
    pylab.semilogx()
    pylab.figure()
    pylab.plot(xAxis, ratiosSDs, 'bo')
    pylab.title('SD Heads/Tails Ratios ('
                + str(numTrials) + ' Trials)')
    pylab.xlabel('Number of Flips')
    pylab.ylabel('Standard Deviation')
    pylab.semilogx()
    pylab.semilogy()

def print_ROC(multi_result, n_imgs, save_filename = None, show_curve = True,
              xmin = None, ymin = None, xmax = None, ymax = None,
              grid_major = False, grid_minor = False):
    points = []
    n_imps = float(n_imgs)
    for result in multi_result:
        tp = float(result.n_true_positives())
        fp = float(result.n_false_positives())
        fn = float(result.n_false_negatives())
        #n_imgs = float(len(result.images))
        points.append((fp / n_imgs, tp / (tp + fn)))
    points.sort()
    if save_filename != None:
        f = open(save_filename, "w")
        cPickle.dump(points, f)
        f.close()
    if show_curve:
        pylab.semilogx([a[0] for a in points], [a[1] for a in points])
        if xmin != None:
            pylab.axis(xmin = xmin)
        if xmax != None:
            pylab.axis(xmax = xmax)
        if ymin != None:
            pylab.axis(ymin = ymin)
        if ymax != None:
            pylab.axis(ymax = ymax)
        pylab.xlabel("False positives per image (FPPI)")
        pylab.ylabel("Detection rate")
        pylab.show()

 def plot_mass_magnitude(self, mag, savefile=None):
     """
     Make a standard mass-luminosity relation plot for this isochrone.
     
     Parameters
     ----------
     mag : string
         The name of the magnitude column to be plotted.
     savefile : string (default None)
         If a savefile is specified, then the plot will be saved to that file. 
     """
     plt.clf()
     plt.semilogx(self.points['mass'], self.points[mag], 'k.')
     plt.gca().invert_yaxis()
     plt.xlabel(r'Mass (M$_\odot$)')
     plt.ylabel(mag + ' (mag)')
     
     fmt_title = 'logAge={0:.2f}, d={1:.2f} kpc, AKs={2:.2f}'
     plt.title(fmt_title.format(self.points.meta['LOGAGE'],
                               self.points.meta['DISTANCE']/1e3,
                               self.points.meta['AKS']))
     
     if savefile != None:
         plt.savefig(savefile)
     
     return