def plotHousing(impression):
    f=open('midWestHousingPrices.txt','r')
    labels,prices=[],[]
    for line in f:
        year,quarter,price=line.split(' ')
        label=year[2:4]+'\n'+quarter[:]
        labels.append(label)
        prices.append(float(price)/1000)
    quarters=pylab.arange(len(labels))
    width=0.8
    if impression=='flat':
        pylab.semilogy()
    pylab.bar(quarters,prices,width)
    pylab.xticks(quarters+width/2.0,labels)
    pylab.title('Housing Prices in U.S. Midwest')
    pylab.xlabel('Quarter')
    pylab.ylabel('Average Price ($1000\'s)')
    if impression=='flat':
        pylab.ylim(10,10**3)
    elif impression=='volatile':
        pylab.ylim(180,220)
    elif impression=='fair':
        pylab.ylim(150,250)
    else:
        raise ValueError

    def plot(self, outputDirectory):
        """
        Plot both the raw kinetics data and the Arrhenius fit versus 
        temperature. The plot is saved to the file ``kinetics.pdf`` in the
        output directory. The plot is not generated if ``matplotlib`` is not
        installed.
        """
        # Skip this step if matplotlib is not installed
        try:
            import pylab
        except ImportError:
            return

        Tlist = 1000.0/numpy.arange(0.4, 3.35, 0.05)
        klist = numpy.zeros_like(Tlist)
        klist2 = numpy.zeros_like(Tlist)
        for i in range(Tlist.shape[0]):
            klist[i] = self.reaction.calculateTSTRateCoefficient(Tlist[i])
            klist2[i] = self.reaction.kinetics.getRateCoefficient(Tlist[i])

        order = len(self.reaction.reactants)
        klist *= 1e6 ** (order-1)
        klist2 *= 1e6 ** (order-1)

        pylab.semilogy(1000.0 / Tlist, klist, 'ok')
        pylab.semilogy(1000.0 / Tlist, klist2, '-k')
        pylab.xlabel('1000 / Temperature (1000/K)')
        pylab.ylabel('Rate coefficient ({0})'.format(self.kunits))
        pylab.savefig(os.path.join(outputDirectory, 'kinetics.pdf'))
        pylab.close()

def plotEventFlop(library, num, eventNames, sizes, times, events, filename = None):
  from pylab import legend, plot, savefig, semilogy, show, title, xlabel, ylabel
  import numpy as np

  arches = sizes.keys()
  bs     = events[arches[0]].keys()[0]
  data   = []
  names  = []
  for event, color in zip(eventNames, ['b', 'g', 'r', 'y']):
    for arch, style in zip(arches, ['-', ':']):
      if event in events[arch][bs]:
        names.append(arch+'-'+str(bs)+' '+event)
        data.append(sizes[arch][bs])
        data.append(1e-3*np.array(events[arch][bs][event])[:,1])
        data.append(color+style)
      else:
        print 'Could not find %s in %s-%d events' % (event, arch, bs)
  semilogy(*data)
  title('Performance on '+library+' Example '+str(num))
  xlabel('Number of Dof')
  ylabel('Computation Rate (GF/s)')
  legend(names, 'upper left', shadow = True)
  if filename is None:
    show()
  else:
    savefig(filename)
  return

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 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 qdisk_plot(root):

	# some labels
	ylabels = ["Heating", r"$N_{\mathrm{hit}}$", r"$N_{\mathrm{hit}}/N_{\mathrm{tot}}$", 
	r"$T_{\mathrm{heat}}$", r"$T_{\mathrm{irrad}}$", r"$W_{\mathrm{irrad}}$"]

	log_lin = [1,0,0,1,1,1]


	p.figure(figsize=(9,10))

	disk_diag = "diag_%s/%s.disk.diag" % (root, root)

	# read the disk_diag file
	a = ascii.read(disk_diag)

	# cyce through the physical quantities and plot for each annulus
	for j, name in enumerate(a.colnames[3:]):
		
		p.subplot(3,2,j+1)
		p.plot(a[name], ls="steps", c="k", linewidth=2)
		p.ylabel(ylabels[j])
		p.xlabel("Annulus")

		if log_lin[j]:
			p.semilogy()

	p.savefig("qdisk_%s.png" % root, dpi=300)

def make_plots():
    work_dir = "/u/cmutnik/work/upperSco_copy/finished/"

    # Read in data
    files = glob.glob(work_dir + "*.fits")

    specs = []

    for ff in range(len(files)):
        spec = fits.getdata(files[ff])

        if ff == 0:
            tot0 = spec[1].sum()

        spec[1] *= tot0 / spec[1].sum()

        specs.append(spec)

    # Plot
    plt.clf()
    for ff in range(len(files)):
        legend = files[ff].split("/")[-1]
        plt.semilogy(specs[ff][0], specs[ff][1], label=legend)

    plt.legend(loc="lower left")
    plt.xlim(0.7, 2.55)

    return

def test_calcpow():
    N1 = 128
    N2 = 128
    t1 = numpy.arange(N1)
    t2 = numpy.arange(N2)
    y1 = numpy.sin(t1*16.*numpy.pi/N1) + numpy.cos(t1*64.*numpy.pi/N1)
    y2 = numpy.sin(t2*16.*numpy.pi/N2) + numpy.sin(t2*32.*numpy.pi/N1)
    
    x = y1[:,None]*y2[None,:]
    x += 0.1*numpy.random.normal(size=(N1,N2))
    
    dt = 2.0

    ell,Pl = calcpow(x,dt,Nl=100)

    pylab.figure()
    pylab.imshow(x)
    pylab.colorbar()

    pylab.figure()
    pylab.semilogy(ell,Pl)
    
    i = numpy.argmax(Pl)
    print "scale of Pmax: %.3g arcmin" % (180.*60./ell[i])
    
    pylab.show()

def plot_track_props(tracks, nx, ny, len_cutoff=20):
    pl.ioff()
    wdist = wraparound_dist(nx, ny)
    val_fig = pl.figure()
    area_fig = pl.figure()
    psn_fig = pl.figure()
    delta_vals = []
    delta_dists = []
    for tr in tracks:
        if len(tr) < len_cutoff:
            continue
        idxs, regs = zip(*tr)
        delta_vals.extend([abs(regs[idx].val - regs[idx + 1].val) for idx in range(len(regs) - 1)])
        dists = [wdist(regs[i].loc, regs[i + 1].loc) for i in range(len(regs) - 1)]
        delta_dists.extend([abs(dists[idx] - dists[idx + 1]) for idx in range(len(dists) - 1)])
        pl.figure(val_fig.number)
        pl.plot(idxs, [reg.val for reg in regs], "s-", hold=True)
        pl.figure(area_fig.number)
        pl.semilogy(idxs, [reg.area for reg in regs], "s-", hold=True)
        pl.figure(psn_fig.number)
        pl.plot(idxs[:-1], dists, "s-", hold=True)
    pl.figure(val_fig.number)
    pl.savefig("val_v_time.pdf")
    pl.figure(area_fig.number)
    pl.savefig("area_v_time.pdf")
    pl.figure(psn_fig.number)
    pl.savefig("psn_v_time.pdf")
    pl.figure()
    pl.hist(delta_vals, bins=pl.sqrt(len(delta_vals)))
    pl.savefig("delta_vals.pdf")
    pl.figure()
    pl.hist(delta_dists, bins=pl.sqrt(len(delta_dists)))
    pl.savefig("delta_dists.pdf")
    pl.close("all")

def demo():
    import pylab

    # The module normalize is not part of the osrefl code base.
    from reflectometry.reduction import normalize

    from .examples import ng7 as dataset
    spec = dataset.spec()[0]
    water = WaterIntensity(D2O=20,probe=spec.probe)
    spec.apply(normalize())
    theory = water.model(spec.Qz,spec.detector.wavelength)

    pylab.subplot(211)
    pylab.title('Data normalized to water scattering (%g%% D2O)'%water.D2O)
    pylab.xlabel('Qz (inv Ang)')
    pylab.ylabel('Reflectivity')
    pylab.semilogy(spec.Qz,theory,'-',label='expected')
    scale = theory[0]/spec.R[0]
    pylab.errorbar(spec.Qz,scale*spec.R,scale*spec.dR,fmt='.',label='measured')

    spec.apply(water)
    pylab.subplot(212)
    #pylab.title('Intensity correction factor')
    pylab.xlabel('Slit 1 opening (mm)')
    pylab.ylabel('Incident intensity')
    pylab.yscale('log')
    pylab.errorbar(spec.slit1.x,spec.R,spec.dR,fmt='.',label='correction')

    pylab.show()

def new_draw_parcel_trace(Tb, PLCL, Press):

	# Convert Pressures to log scale
	Pfact = np.multiply(skewness,np.log10(np.divide(1000., Press)))

	parcelT = []
	flag = 1

	for p in range(len(Press)):
		if Press[p] >= PLCL:
			newTB = ((Tb + 273.) * (Press[p]/Press[0]) ** (287.04/1004.)) - 273.
			parcelT.append(newTB)
		else:
			if flag:
				if p == 0:
					moists = draw_moist_adiabats(0, 1, Tb, 0)
				else:
					moists = draw_moist_adiabats(0,1,parcelT[p-1], (p - 1 + len(press_levels) - len(Press)))
				for m in moists:
					parcelT.append(m)
				flag = 0


	minlen = min(len(parcelT), len(Pfact))
	
	dry_parcel_trace = np.add(parcelT[:minlen], Pfact[:minlen])



	pylab.semilogy(dry_parcel_trace,Press[:minlen],\
		basey=10, color = 'brown', linestyle = 'dotted',\
		linewidth = 1.5)

def showGrowth(lower, upper):
    log = []
    linear = []
    quadratic = []
    logLinear = []
    exponential = []
    for n in range(lower, upper+1):
        log.append(math.log(n, 2))
        linear.append(n)
        logLinear.append(n*math.log(n, 2))
        quadratic.append(n**2)
        exponential.append(2**n)
    pylab.plot(log, label = 'log')
    pylab.plot(linear, label = 'linear')
    pylab.legend(loc = 'upper left')
    pylab.figure()
    pylab.plot(linear, label = 'linear')
    pylab.plot(logLinear, label = 'log linear')
    pylab.legend(loc = 'upper left')
    pylab.figure()
    pylab.plot(logLinear, label = 'log linear')
    pylab.plot(quadratic, label = 'quadratic')
    pylab.legend(loc = 'upper left')
    pylab.figure()
    pylab.plot(quadratic, label = 'quadratic')
    pylab.plot(exponential, label = 'exponential')
    pylab.legend(loc = 'upper left')
    pylab.figure()
    pylab.plot(quadratic, label = 'quadratic')
    pylab.plot(exponential, label = 'exponential')
    pylab.semilogy()
    pylab.legend(loc = 'upper left')
    return

def demo_perfidious(n):
    plt.figure()

    r = (np.arange(n)+1)/float(n+1)

    bases = [(PowerBasis(), "Power"),
             (ChebyshevBasis(interval=(1./(n+1),n/float(n+1))),"Chebyshev"), 
             (LagrangeBasis(interval=(1./(n+1),n/float(n+1))),"Lagrange"), 
             (LagrangeBasis(r),"Specialized Lagrange")]

    xs = np.linspace(0,1,50*n)
    
    for (i,(b,l)) in enumerate(bases):
        p = b.from_roots(r)
        plt.subplot(len(bases),1,i+1)
        plt.semilogy(xs,np.abs(p(xs)),label=l)
        plt.xlim(0,1)
        plt.ylim(min=1)
        
        for j in range(n):
            plt.axvline((j+1)/float(n+1),linestyle=":",color="black")
        plt.legend(loc="best")
    print b.points
    print p.coefficients
    plt.subplot(len(bases),1,1)
    plt.title('The "perfidious polynomial" for n=%d' % n)

def hanning_standard_plot(data, rate):
    sample_length = len(data)
    k = arange(sample_length)
    period = sample_length / rate
    freqs = (k / period)[range(sample_length / 2)] #right-side frequency range
    Y = (fft(data * np.hanning(sample_length)) / sample_length)[range(sample_length / 2)]
    semilogy(freqs, abs(Y))

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.ylabel('Abs(#Heads - #Tails)')
    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.semilogx()
    pylab.semilogy()
    pylab.plot(xAxis, ratios, 'bo')

def standard_plot(data, rate):
    sample_length = len(data)
    k = arange(sample_length)
    period = sample_length / rate
    freqs = (k / period)[range(sample_length / 2)] #right-side frequency range
    Y = (fft(data) / sample_length)[range(sample_length / 2)]
    semilogy(freqs, abs(Y)) # plotting the spectrum

def plotDependencyEPS():

    """Plot thoretical dependency between n_components and eps"""
    
    # range of admissible distortions
    eps_range = np.linspace(0.01, 0.99, 100)

    # range of number of samples to embed
    n_samples_range = np.logspace(2, 6, 5)
    colors = pl.cm.Blues(np.linspace(0.3, 1.0, len(n_samples_range)))

    pl.figure()

    for n_samples, color in zip(n_samples_range, colors):
        min_n_components = johnson_lindenstrauss_min_dim(n_samples, \
                                                         eps=eps_range)
        pl.semilogy(eps_range, min_n_components, color=color)

    pl.legend(["n_samples = %d" % n for n in n_samples_range], \
              loc="upper right")

    pl.xlabel("Distortion eps")
    pl.ylabel("Minimum number of dimensions")
    pl.title("Johnson-Lindenstrauss bounds:\nn_components vs eps")
    pl.show()

	def plotHist(self):
		p.figure()
		p.semilogy(self.counts.keys(), self.counts.values(), '.')
		p.xlabel('Log-return')
		p.ylabel('Count')
		p.title(self.symbol)
		p.show()

def testPlot2(trials=51, maxsteps=5000):
    f = FunctionWrapper(trials, OptimumJumper(StochQuad(noiseLevel=10, curvature=1), jumptime=1000, jumpdist_std=1))
    for aclass, aparams in [#(SGD, {'learning_rate':0.1}),
                            #(SGD, {'learning_rate':0.01}),
                            #(AveragingSGD, {'learning_rate':0.01}),
                            #(AveragingSGD, {'learning_rate':0.01, 'fixedDecay':0.1}),
                            #(AveragingSGD, {'learning_rate':0.01, 'fixedDecay':0.1}),
                            #(AveragingSGD, {'learning_rate':0.1}),
                            #(AveragingSGD, {'learning_rate':1.0}),
                            (AveragingOracle, {}),
                            (AveragingOracle, {"fixedDecay":0.1}),
                            #(AveragingOracle, {"fixedDecay":0.01}),
                            (AdaptivelyAveragingOracle, {}),
                            #(AdaGrad, {'init_lr':0.3}),
                            #(Amari, {'init_lr':0.1, 'time_const':100}),
                            #(RMSProp, {'init_lr':0.1}),
                            (OracleSGD, {}),
                            #(vSGD, {'verbose':False}),
                            #(vSGDfd, {}),
                            ]:
        ls = lossTraces(fwrap=f, aclass=aclass, dim=trials,
                        maxsteps=maxsteps, algoparams=aparams)
        plotWithPercentiles(ls, algo_colors[aclass], aclass.__name__)
    pylab.semilogy()
    pylab.xlim(0, maxsteps)
    pylab.legend()
    pylab.show()