def timeseries(bs, dt=1):
    if not hasattr(bs, '__iter__'):
        bs = [bs,]
    nrow = len(bs)
    for (i, b) in enumerate(bs):
        p = params(b)
        plt.subplot(nrow, 1, 1+i)
        v = pop212(b).sum(axis=1) # stack over channels
        nt = len(v)
        dt = min(dt, nt)
        nt = nt - np.fmod(nt, dt) # fit time segments after decimation
        v = v[:nt].reshape((nt/dt, -1)).sum(axis=1) # clip to multiple of dt and stack
        t = p.dtvec[:nt].reshape((-1, dt)).mean(axis=1) + p.T/2.
        amp = np.abs(v)
        phase = np.angle(v)
        plt.plot(t, amp, 'b.-')
        plt.ylim(0, plt.ylim()[1])
        plt.gca().set_yticklabels([])
        plt.twinx()
        plt.plot(t, phase, 'r.-')
        plt.ylim(-np.pi, np.pi)
        plt.gca().set_yticklabels([])
        putil.rmgaps(1e6, 2.0)
        plt.xlim(0, p.T)
        plt.gca().add_artist(AnchoredText(p.baseline, loc=1, frameon=False, borderpad=0))
    plt.setp(plt.gcf(), figwidth=8, figheight=2+nrow)
    plt.tight_layout()
    plt.subplots_adjust(hspace=0)

def log_plot():
    start = int(np.floor(np.log(min(median)) / np.log(10))) + 3
    end = int(np.ceil(np.log(max(median)) / np.log(10))) + 3

    xs = []
    ticks = []
    for i in range(start, end + 1):
        xs.append(10 ** (i - 3))

        if i % 3 == 0:
            ticks.append('{}s'.format(prefix[i / 3]))
        else:
            ticks.append(str(10 ** (i % 3)))

    plt.barh(pos, median, align='center', height=0.25, left=1e-3,
             color=bar_color, lw=0)
    plt.errorbar(median, pos, ecolor=error_bar_color, fmt=None, xerr=err)

    plt.grid(True)
    plt.xlabel('Time')
    plt.xlim(min(xs), max(xs))
    plt.xscale('log')
    plt.xticks(xs, ticks)
    plt.ylim(ymax=size)
    plt.yticks(pos, language)

    plt.twinx()
    plt.ylim(ymax=size)
    plt.yticks(pos, relative)

    plt.savefig('plots/{}.png'.format(pid), bbox_inches='tight')
    plt.clf()

def main():

    # switch the commented lines here to alternate between CV testing and making kaggle submission
    x_train, x_test, y_train, y_test = load_data_cv('data/train_32.npy')
    #x_train, y_train, x_test = load_data_test('data/train_32.npy', 'data/test_32.npy')

    model = build_model()

    print("Starting training")
    # batch iterator with 300 epochs
    train_loss = []
    valid_loss = []
    valid_acc = []
    for i in range(300):
        loss = batch_iterator(x_train, y_train, 128, model)
        train_loss.append(loss)
        valid_avg = model.evaluate(x_test, y_test, show_accuracy = True, verbose = 0)
        valid_loss.append(valid_avg[0])
        valid_acc.append(valid_avg[1])
        print 'epoch:', i, 'train loss:', np.round(loss, decimals = 4), 'valid loss:', np.round(valid_avg[0], decimals = 4), 'valid acc:', np.round(valid_avg[1], decimals = 4)

    train_loss = np.array(train_loss)
    valid_loss = np.array(valid_loss)
    valid_acc = np.array(valid_acc)
    sns.set_style("whitegrid")
    pyplot.plot(train_loss, linewidth = 3, label = 'train loss')
    pyplot.plot(valid_loss, linewidth = 3, label = 'valid loss')
    pyplot.legend(loc = 2)
    pyplot.ylim([0,4.5])
    pyplot.twinx()
    pyplot.plot(valid_acc, linewidth = 3, label = 'valid accuracy', color = 'r')
    pyplot.grid()
    pyplot.ylim([0,1])
    pyplot.legend(loc = 1)
    pyplot.show()

def fidelity_vs_power(folder='', x_axis_par='par_ro_Ex_power'):
    
    if folder == '':
        folder = os.getcwd()

    allfiles = os.listdir(folder)
    fidfiles = [f for f in allfiles if ('totalfid' in f and '.dat' in f)]

    x_axis = [] #in this case powers
    maxfid = []
    maxfid_t = []

    for f in fidfiles:
        fn, ext = os.path.splitext(f)
        idx = int(fn[fn.find('+_')+2:])
        basepath = os.path.join(folder, PREFIX+'-'+str(idx))
        parfile = basepath+'_'+PARAMS_SUFFIX+'.dat'  
        x_axis.append(loadtxt(parfile)[get_param_column(parfile,x_axis_par)]*1e9) #pwr in nW
        
        fiddat = loadtxt(f)
        maxidx = argmax(fiddat[1,:])
        maxfid.append(fiddat[1,maxidx])
        maxfid_t.append(fiddat[0,maxidx])
    
    
    fig = plt.figure()
    plt.plot(x_axis, maxfid, 'ro', label='max F')
    plt.xlabel('P [nW]')
    plt.ylabel('max. F')
    plt.legend()
    
    plt.twinx()
    plt.plot(x_axis, maxfid_t, 'bo')
    plt.ylabel('best ro-time')
    plt.savefig('fidelity_vs_power.png')

def plt_twin( axis, tick0=None, tick=None ) : 
	'''
	Add x-top or y-right axis

	axis:
		['x' | 'y']

	tick0:
		Must between [0, 1]
	'''
	if (str(axis).lower() not in ['x', 'y']) : Raise(Warning, "axis='"+str(axis)+"' not in ['x', 'y']. Do nothing !")
	axis = str(axis).lower()
	#--------------------------------------------------
	if (tick0 is not None and tick is not None) : 
		tick0 = npfmt(tick0, float)
		if (tick0.min()<0 or tick0.max()>1) : 
			Raise(Warning, 'tick0.(min,max)=(%.1f, %.1f) out of [0, 1]. Do nothing !' % (tick0.min(), tick0.max()))
		else : 
			if   (axis == 'x') : 
				plt.twiny()
				plt.xticks(tick0, tick)
			elif (axis == 'y') : 
				plt.twinx()
				plt.yticks(tick0, tick)
	#--------------------------------------------------
	elif (tick0 is None and tick is None) : 
		if   (axis == 'x') : plt.tick_params(axis='x', which='both', labeltop='on', labelbottom='on')
		elif (axis == 'y') : plt.tick_params(axis='y', which='both', labelleft='on', labelright='on')
	#--------------------------------------------------
	else : Raise(Warning, 'tick0, tick must both ==None or !=None, now one is None but the other is not. Do nothing !')

    def plot_filter_characteristics(self):
        w, h = freqz(self.freq_filter.num, self.freq_filter.denom)
        plt.figure(1)
        plt.subplot(2,1,1)
        plt.hold(True)
        powa = plt.plot((self.filter_parameters.sample_rate*0.5/pi)*w, abs(h),'b-', label = 'Char. amplitudowa')
        plt.title('Charakterystyki filtru')
        plt.xlabel('Czestotliwosc [Hz]')
        plt.ylabel('Amplituda')


        plt.twinx(ax=None)
        angles = unwrap(angle(h))
        plt.znie = plot((self.filter_parameters.sample_rate*0.5/pi)*w,angles, 'g-', label = 'Char. fazowa')
        plt.ylabel('Faza')

        plt.grid()
        tekst = powa + znie
        wybierz = [l.get_label() for l in tekst]

        plt.legend(tekst, wybierz, loc='best')
    ########################################################################################################################
        plt.subplot(2,1,2)

        w2, gd = group_delay((num, denom))

        plt.plot((sample_rate*0.5/pi)*w2, gd)
        plt.grid()
        plt.xlabel('Czestotliwosc [Hz]')
        plt.ylabel('Opoznienie grupowe [probki]')
        plt.title('Opoznienie grupowe filtru')

        plt.show()

 def pcs(self):
     self.pesobj = PES.PES(self.X,self.Y,self.S,self.D,self.lb.flatten(),self.ub.flatten(),self.para['kindex'],self.para['mprior'],self.para['sprior'],DH_SAMPLES=self.para['DH_SAMPLES'], DM_SAMPLES=self.para['DM_SAMPLES'], DM_SUPPORT=self.para['DM_SUPPORT'],DM_SLICELCBPARA=self.para['DM_SLICELCBPARA'],mode=self.para['SUPPORT_MODE'])
     xmin = self.reccomend()
     plt.figure(1)
     plt.plot(xmin[0],xmin[1],'r.')
     print xmin
     plt.figure(2)
     plt.subplot(4,1,1)
     ns = 6000
     sup = sp.linspace(-1,1,ns)
     for i in xrange(2):
         X = sp.vstack([xmin for k in xrange(ns)])
         print X.shape
         for j in xrange(ns):
             X[j,i] = sup[j]
         [m,v] = self.pesobj.G.infer_diag_post(X,[[sp.NaN]]*ns)
         s = sp.sqrt(v)
         plt.subplot(4,1,2*i+1)
         plt.fill_between(sup,(m-2*s).flatten(),(m+2*s).flatten(), facecolor='lightblue',edgecolor='lightblue')
         plt.plot(sup,m.flatten())
         [m,v] = self.pesobj.G.infer_diag_post(X,[[i]]*ns)
         s = sp.sqrt(v)
         plt.subplot(4,1,2*i+2)
         plt.fill_between(sup,(m-2*s).flatten(),(m+2*s).flatten(), facecolor='lightblue',edgecolor='lightblue')
         plt.plot(sup,m.flatten(),'r')
         p = sp.exp(-0.5*(m**2)/v)
         
         
         
         plt.twinx().plot(sup,p.flatten(),'g')
     return

def laserPlot(filenameLst):
	'''
	Plot the temporal and spectral profile of a
	mode-locked fiber laser simulation
	'''

	nbrPlots = len(filenameLst)
	for i in arange(len(filenameLst)):
		results = load(filenameLst[i])
		t = results['t']
		nt = results['nt']
		T = results['T']
		archivePass = results['archivePass']
		nu_inst_out3 = results['nu_inst_out3']
		nu_inst_out4 = results['nu_inst_out4']
		spectre_out = results['spectre_out']
		wavelength = results['wavelength']


		# Graph
		plt.figure(figsize=(12,9))

		ax3 = plt.subplot(221)
		plt.plot(t, pow(abs(archivePass[0]),2), color="black")
		plt.ylabel("$|u(z,T)|^2$ [W]")
		plt.xlabel("$T/T_0$")
		plt.xlim([-T/2,T/2])
		plt.grid(True)
		ax4 = plt.twinx()
		plt.plot(t[0:nt-1], nu_inst_out3)
		plt.ylabel("Chirp")
		ax4.yaxis.tick_right()
		plt.ylim([-1.5,1.5])

		ax5 = plt.subplot(223)
		plt.semilogy(t, pow(abs(archivePass[0]),2), color="black")
		plt.ylabel("$|u(z,T)|^2$ [dBm]")
		plt.xlabel("$T/T_0$")
		plt.xlim([-T/2,T/2])
		plt.grid(True)
		ax4 = plt.twinx()
		plt.plot(t[0:nt-1], nu_inst_out3)
		plt.ylabel("Chirp")
		ax4.yaxis.tick_right()
		plt.ylim([-1.5,1.5])

		ax7 = plt.subplot(222)
	 	plt.plot(wavelength, spectre_out, color="black")
		plt.xlabel("Wavelength [nm]")
		plt.grid(True)

		ax8 = plt.subplot(224)
	 	plt.semilogy(wavelength, spectre_out, color="black")
		plt.xlabel("$T/T_0$")
		plt.xlabel("Wavelength [nm]")
		plt.grid(True)

	plt.show()

def makeFig():
    plt.title('recycle')
    plt.plot(xArray, '-ro', label='paper')
    plt.legend(loc='upper left')
    plt2=plt.twinx()
    plt2.plot(yArray, 'b^-', label='plastic')
    plt2.legend(loc='upper center')
    plt3=plt.twinx()
    plt3.plot(zArray, 'gD-', label='can')
    plt3.legend(loc='lower left')

def trialPlot(dm, soa, _show=show, err=True, minSmp=200, suffix='', padding=0,
	diff=True):

	"""
	A pupil-trace plot for the full trial epoch.

	Arguments:
	dm				--	A DataMatrix.
	soa				--	The SOA to select.

	Keyword arguments:
	_show			--	Indicates whether the plot should be shown.
						(default=True)
	err				--	Indicates whether error bars should be drawn.
						(default=True)
	suffix			--	A suffix to identify the trace. (default='')
	padding			--	A padding time to be added to the traceLen. (default=0)
	diff			--	Indicates whether the difference trace should be plotted
						as well. (default=True)
	"""

	assert(soa in dm.unique('soa'))
	if _show:
		Plot.new(size=Plot.ws)
		plt.title('SOA: %d ms' % (soa+55))
	plt.axhline(1, linestyle='--', color='black')
	dm = dm.select('soa == %d' % soa)
	# Determine the trace length and create the trace plot
	traceLen = soa + 105 + padding
	traceParams = trialParams.copy()
	traceParams['traceLen'] = traceLen
	tracePlot(dm, traceParams=traceParams, err=err, suffix='.%d%s' % (soa, \
		suffix), minSmp=minSmp)
	# Cue
	plt.axvspan(0, cueDur, color=blue[1], alpha=.2)
	# Target. Take into account to cue duration in determining the target onset.
	targetOnset = soa+55
	plt.axvspan(targetOnset, targetOnset+targetDur, color=blue[1], alpha=.2)
	plt.xlim(0, 2550)
	plt.legend(frameon=False)
	plt.xlabel('Time since cue onset (ms)')
	plt.ylabel('Pupil size (norm.)')
	plt.yticks([1,1.025, 1.05])
	plt.xticks(range(0, 2501, 500))
	if diff:
		plt.ylim(.92, 1.07)
		plt.axhline(diffY, linestyle='--', color='black')
		plt.twinx()
		plt.tick_params(axis="y")
		plt.ylim(.92, 1.07)
		plt.yticks([.925,.95, .975], [-.025, 0, .025])
	else:
		plt.ylim(.98, 1.07)
	if _show:
		Plot.save('trialPlot.%d' % soa, 'trialPlot', show=show)

def plot_label_by_success(setup_name, school, context_name=None, term_type=None, legend=True, linetype='-', show_data_size=True, set_order=None):
    data = load_ratings_with_contexts()
    if set_order is not None:
        data = data[data['practice_set_order'] == set_order]
    data = data[(data['experiment_setup_name'] == setup_name)]
    if context_name is not None:
        data = data[data['context_name'] == context_name]
    if term_type is not None:
        data = data[data['term_type'] == term_type]
    if school is not None:
        school_usage = load_school_usage().reset_index().rename(columns={'ip_address': 'school', 'user_id': 'user'})
        data = pandas.merge(data, school_usage, on='user', how='inner')
        data = data[data['school'] == school]
    data = data[data['error_rate'].apply(lambda x: x % 10 == 0)]

    def _apply(group):
        result = []
        for label in group['label'].unique():
            mean = binomial_confidence_mean(group['label'] == label)
            result.append({
                'label': label,
                'learners': 100 * mean[0],
                'learners_min': 100 * mean[1][0],
                'learners_max': 100 * mean[1][1],
            })
        return pandas.DataFrame(result)
    to_plot = data.groupby(['experiment_setup_name', 'error_rate']).apply(_apply).reset_index().sort_values(by=['label', 'error_rate'])
    for i, (label, label_data) in enumerate(to_plot.groupby('label')):
        plt.plot(
            label_data['error_rate'],
            label_data['learners'],
            linetype,
            label=label.split('-')[-1],
            color=output.palette()[i],
            marker='.',
            markersize=20
        )
        plt.fill_between(
            label_data['error_rate'],
            label_data['learners_min'],
            label_data['learners_max'],
            color=output.palette()[i], alpha=0.35
        )
    if legend:
        plt.legend(ncol=3, loc='upper left', frameon=True)
    plt.ylabel('Label (%)')
    plt.xlabel('Real error rate')
    plt.gca().xaxis.grid(True)
    plt.gca().yaxis.grid(True)
    if show_data_size:
        plt.twinx()
        size = data.groupby('error_rate').apply(len).reset_index().rename(columns={0: 'size'})
        plt.plot(size['error_rate'], size['size'], '.-', color='gray')
        plt.ylabel('Data size')
    plt.ylim(0, 70)

	def sampleExponentials(self, n_samples = 1000):
		la = np.load("results/fitDistLog.npy")
		x = np.linspace(0, 1, 100)
		av = np.average(la[:, :2], axis=0)
		cov = np.cov(la[:, :2].T)
		s = np.random.multivariate_normal(av, cov, size=n_samples)
		vals = s[:, 0]*np.exp(np.outer(x, s[:, 1]))
		print vals.shape
		bins, _, _ = self.getBinnedData(2, 100)
		'''for i in xrange(100):
			if len(bins[i]) != 0:
				plt.hist(bins[i], bins=50, alpha=0.6, normed=True)
			plt.hist(vals[i], bins=50, alpha=0.6, color='r', normed=True)
			plt.savefig("results/comparefullhist{0:03d}.pdf".format(i))
			plt.clf()'''
		plt.subplot(231)
		hist, xedges, yedges = np.histogram2d(la[:, 0], la[:, 1], bins=20, normed=True)
		plt.imshow(hist.T, interpolation = 'nearest', cmap=cm.Blues, aspect = 'auto',
				extent=[np.min(xedges), np.max(xedges), np.min(yedges), np.max(yedges)],
				origin='lower')
		plt.subplot(232)
		samples = np.random.multivariate_normal(av, cov, size=10*n_samples)
		hist_sp, xedges_sp, yedges_sp = np.histogram2d(samples[ :, 0], 
									samples[ :, 1], bins=50, normed=True)
		plt.imshow(hist_sp.T, interpolation = 'bicubic', cmap=cm.Reds, aspect = 'auto',
				extent=[np.min(xedges_sp), np.max(xedges_sp), np.min(yedges_sp), np.max(yedges_sp)],
				origin='lower')
		plt.subplot(233)
		plt.pcolor(np.linspace(np.min(xedges), np.max(xedges), len(xedges)),
				np.linspace(np.min(yedges), np.max(yedges), len(yedges)), hist.T,
				alpha=0.8, cmap=cm.Blues)
		plt.pcolor(np.linspace(np.min(xedges_sp), np.max(xedges_sp), len(xedges_sp)),
				np.linspace(np.min(yedges_sp), np.max(yedges_sp), len(yedges_sp)), hist_sp.T,
				alpha=0.5, cmap=cm.Reds)
		plt.subplot(234)
		plt.hist2d(la[:, 0], la[:, 1], bins=20, cmap=cm.Blues)
		plt.subplot(235)
		_, bx, _ = plt.hist(la[:, 1], 20, alpha=0.6, normed=True)
		plt.hist(samples[:, 1], 20, color='r', alpha=0.6, normed=True)
		xp = np.linspace(np.min(bx), np.max(bx), 100)
		p = scipy.stats.norm.pdf(xp, loc=av[1], scale=np.sqrt(cov[1,1]))
		plt.plot(xp, p)
		plt.subplot(236)
		nx, bx, _ = plt.hist(la[:, 0], 20, alpha=0.6, normed=True)
		plt.hist(samples[:, 0], 20, color='r', alpha=0.6, normed=True)
		xp = np.linspace(np.min(bx), np.max(bx), 100)
		p = scipy.stats.norm.pdf(xp, loc=av[0], scale=np.sqrt(cov[0,0]))
		plt.plot(xp, p)
		plt.twinx()
		plt.plot(bx[1:]-np.diff(bx), np.cumsum(nx)*np.diff(bx), color='g')
		plt.savefig("results/comparehistdist.pdf")
		plt.clf()

def fidelity_vs_power(folder='',sweep_param='Ex_RO_amplitude'):
    
    if folder == '':
        folder = os.getcwd()

    allfiles = os.listdir(folder)
    fidfiles = [f for f in allfiles if ('totalfid' in f and '.dat' in f)]

    pow = []
    maxfid = []
    maxfid_t = []

    
    for f in fidfiles:
        fn, ext = os.path.splitext(f)
        idx = int(fn[fn.find('+_')+2:])-1
        print idx
        
        if idx < 10:
            SUFFIX = '-00'+str(idx)
        elif idx == 10:
            SUFFIX = '-0'+str(idx)
        elif idx < 100:
            SUFFIX = '-0'+str(idx)

        basepath = os.path.join(folder, PREFIX+SUFFIX)
        parfile = basepath+'_parameters_dict.npz'  
        param_dict=load(parfile)
        #pow.append(loadtxt(parfile)[get_param_column(parfile,'par_ro_Ex_power')]*1e9)
        #pow.append(int(idx/2)*1.)
        pow.append(param_dict[sweep_param])
        param_dict.close
        fiddat = loadtxt(f)
        maxidx = argmax(fiddat[1,:])
        maxfid.append(fiddat[1,maxidx])
        maxfid_t.append(fiddat[0,maxidx])
        #fiddat.close() 
    
    fig = plt.figure()
    plt.plot(pow, maxfid, 'ro', label='max F')
    plt.xlabel('P [nW]')
    plt.ylabel('max. F')
    plt.ylim(ymax=1.2*max(maxfid))
    plt.title('Maximum SSRO Fidelity and Optimal readout time vs' + sweep_param)
    plt.text(0.01*(max(pow)+min(pow)),1.15*max(maxfid),folder,fontsize='x-small')
    plt.legend()
    
    plt.twinx()
    plt.plot(pow, maxfid_t, 'bo')
    plt.ylabel('best ro-time')
    plt.ylim(ymax=1.2*max(maxfid_t))
    plt.savefig('fidelity_vs_power.png')

def makeFig(): #Create a function that makes the desired plot
    gs = gridspec.GridSpec(3, 3) #gridspec is created 3x3
    #First plot fig 1
    #Plot 1
    plt.subplot(gs[0, :])                                   #subplot position atributes
    plt.ylim([-30,50])                                      #Set y min and max values
    plt.title('Temperatura em Graus C')                     #Plot the title
    plt.grid(True)                                          #Turn on the grid
    plt.ylabel('Temp-1 c')                                  #Set ylabels
    plt.plot(tempF, 'ro-', label='temperatura em graus')    #plot the temperature
    plt.legend(loc='upper left')                            #plot the legend
    plt2=plt.twinx()                                        #Create a second y axis
    plt.ylim(-30,50)                                        #Set limits of second y axis- adjust to readings you are getting
    plt2.plot(tempF2, 'b^-', label='Temp-2 c')              #plot temperature array
    plt2.set_ylabel('Temp-2 c')                             #label second y axis
    plt2.ticklabel_format(useOffset=False)                  #Force matplotlib to NOT autoscale y axis
    plt2.legend(loc='upper right')                          #plot the legend
    #second plot same figure (same window)
    #Plot 2
    plt.subplot(gs[1, :])                                   #subplot position attributes
    plt.ylim([0,100])                                       #Set y min and max values
    plt.title('Humidade do Ar em Percentagem')              #Plot the title
    plt.grid(True)                                          #Turn on the grid
    plt.ylabel('Himidade-1 %')                              #Set ylabels
    plt.plot(humF1, 'ro-', label='Humidade')                #plot the temperature
    plt.legend(loc='upper left')                            #plot the legend
    plt2=plt.twinx()                                        #Create a second y axis
    plt.ylim(0,100)                                         #Set limits of second y axis- adjust to readings you are getting
    plt2.plot(humF2, 'b^-', label='Humidade-2 %')           #plot temperature array
    plt2.set_ylabel('Humidade-2 %')                         #label second y axis
    plt2.ticklabel_format(useOffset=False)                  #Force matplotlib to NOT autoscale y axis
    plt2.legend(loc='upper right')                          #plot the legend
    #Third plot same figure (same window)
    #Plot 3
    plt.subplot(gs[-1,0])                                   #subplot position atributes
    plt.ylim([0,100])                                       #Set y min and max values
    plt.title('Humidade do solo')                           #Plot the title
    plt.grid(True)                                          #Turn on the grid
    plt.ylabel('Himidade %')                                #Set ylabels
    plt.plot(moist, 'ro-', label='Humidade')                #plot the temperature
    plt.legend(loc='upper left')                            #plot the legend
    #Fourth plot same figure (same window)
    #Plot 4
    plt.subplot(gs[-1,-1])                                  #subplot position atributes
    plt.ylim([0,2000])                                      #Set y min and max values
    plt.title('Luminosidade')                               #Plot the title
    plt.grid(True)                                          #Turn on the grid
    plt.ylabel('Luminosidade (lux)')                        #Set ylabels
    plt.plot(lum, 'ro-', label='Luminosidade')              #plot the temperature
    plt.legend(loc='upper left')                            #plot the legend

 def plotsigma(self,emit=2.5e-6/7000*0.938,deltap=1.1e-4,**nargs):
   self.sigx =sqrt(self.betx*emit)*1000
   self.sigy =sqrt(self.bety*emit)*1000
   self.sigdx=self.dx*deltap*1000
   self.plot('sigx sigy sigdx',**nargs)
   ya,yb=pl.ylim()
   pl.twinx()
   bmax=max(self.betx.max(),self.bety.max())
   rng=range(0,int(_n.ceil(_n.log10(bmax)))+1)
   bval=_n.array([n*10**dd for dd in rng for n in [1,2,5] ])
   bval=bval[bval<bmax]
   pl.ylim(ya,yb)
   self._plot=_p.gcf()
   return t

    def payload_loss_figure(self):
        ps = []
        rs = []
        rs2 = []
        Omegas = [i * 0.001 for i in range(1, 201)]
        last_p = 0
        for Omega in Omegas:
            self.udp.Omega = Omega
            p, r = self.udp.get_payload_size()
            if p != last_p:
                print "Omega=%s, p=%s, r=%s" % (Omega, p, r)
                last_p = p
            ps.append(p)
            rs.append(r)

            r2 = self.udp.get_ratio_lower(p=4096)
            rs2.append(r2 / r)

        plt.clf()
        plt.cla()
        l, = plt.plot(Omegas, ps)
        plt.legend([l], ["Optimal\nPayload Size"], loc="upper left", prop={"size": 18})
        # print ps
        # plt.ylim(ymin=0.0, ymax=1)
        plt.ticklabel_format(style="sci", axis="x", scilimits=(3, 3))
        plt.ticklabel_format(style="sci", axis="y", scilimits=(-2, 0))

        plt.grid(True)
        plt.xlabel("Packet Loss Rate")
        plt.ylabel("Optimal Payload Size per Chunk")

        plt.twinx()
        #         l1, = plt.plot(Omegas, rs, "g")
        #         l2, = plt.plot(Omegas, rs2, "r")
        #         plt.ylabel("Max Goodput-to-Throughput Ratio")
        #         plt.legend([l1, l2], ["Max G2T Ratio", "Ratio of 4096-to-LB"], loc="upper right",  prop={'size':18})
        #
        l1, = plt.plot(Omegas, rs, "g")
        # l2, = plt.plot(Omegas, rs2, "r")
        plt.ylabel("Max Goodput-to-Throughput Ratio")
        plt.legend([l1], ["Max G2T Ratio"], loc="upper right", prop={"size": 18})

        plt.title("UDP: $\Delta$=%s, M=%s" % (self.udp.Delta, self.udp.M))

        name = "payload-loss-M%s-Delta%s" % (self.udp.M, self.udp.Delta)
        name.replace(".", "")
        plt.savefig(self.udp.out_dir + "/" + name + ".pdf")

        print "%s.pdf ends" % (name)

    def ex_twinplot(self, X, Y, **kwargs):
#        go_with_x()
#        ax = gca()

        # If type is tuple, use secondary verticle axis, otherwise plot multiple values on same scale
        if type(Y) == tuple:
            if len(Y)>2:
                raise Exception, "twinx plots only takes two arguments"
            else:
                for i in range(len(Y)):
                    self.ex_plot(X, Y[i], **kwargs)
                    twinx()
        else:
            for i in range(len(Y)):
                self.ex_plot(X, Y[i], **kwargs)

def agdd_plots(nplots, iplot, tbase, tmax, t_range, temp, agdd):
    """This function does the AGDD plots in a nplots vertical stack of plots.
    iplot is the current plot (from top to bottom, starting at 1), tbase and
    tmax are the values used for AGDD calculations. t_range is a temporal range
    (usually DoY) and temp and AGDD are extracted 2m temperature and AGDD"""
    plt.subplot(nplots, 1, iplot)
    # Put a grey area for the AGDD calculation bounds
    plt.axhspan(tbase, tmax, xmin=0, xmax=366, color="0.9")
    # Plot temperature
    plt.plot(t_range, temp, "-r", label="Tm")
    plt.ylabel("Mean Temp [degC]")
    plt.grid(True)
    plt.twinx()
    plt.plot(t_range, agdd, "-g", label="AGDD")
    plt.ylabel("AGDD [degC]")

def hist(nndist, **kwds):
    if 'output' in kwds:
        import matplotlib
        matplotlib.use('Agg')
    from matplotlib import pyplot as plt

    # handle formatting keywords
    linewidth = kwds.get('linewidth', 3)
    #plt.rc('text', usetex=True)
    plt.rc('font', family='serif')
    n, bins, patches = plt.hist(nndist, normed=True, bins=100)
    width = bins[1] - bins[0]
    ax1 = plt.gca()
    ax2 = plt.twinx()
    ax2.plot(bins[:-1], width*np.cumsum(n), 'r-', linewidth=linewidth)
    ax2.set_ylim(top=1.0)
    tics = ax1.get_yticks(); ax1.set_yticks(tics[1:])
    tics = ax2.get_yticks(); ax2.set_yticks(tics[1:])
    ax1.set_xlabel(r"d ($\mu$m)")
    ax1.set_ylabel(r"PDF")
    ax2.set_ylabel(r"CDF", color='r')
    for tl in ax2.get_yticklabels():
        tl.set_color('r')
    # check scalar descriptors
    median_dist = np.median(nndist)
    ax1.axvline(median_dist, color='gray', linewidth=linewidth)
    # handle keywords
    if 'title' in kwds:
        plt.title(kwds['title'])
    if 'output' in kwds:
        plt.draw()
        plt.savefig(kwds['output'])
    else:
        plt.show()

def hist_plot_two_y(m, n, p, mv, nv, pv, filename):
    x = np.arange(0,len(m))
    ax = plt.subplot(111)
    common_params = dict(bins=len(m), range=(0,len(m)), normed=False)

    patterns = ('//', '', '\\')
    _, _, patches = plt.hist((x,x,x), weights=(m,n,p), \
            label=["Autocorrelation", "HOS", "PCS-based HOS"], **common_params)
    for patch,pattern in zip(patches, patterns):
        [k.set_hatch(pattern) for k in patch]

    ax2 = plt.twinx()
    ax2.plot(map(lambda i:i+0.5, x), np.array(mv)/np.array(m), label="Autocorrelation", linestyle=':', color='b', marker='v', linewidth=1.5)
    ax2.plot(map(lambda i:i+0.5, x), np.array(nv)/np.array(n), label="HOS", linestyle='-.', color='g', marker='^', linewidth=1.5)
    ax2.plot(map(lambda i:i+0.5, x), np.array(pv)/np.array(p), label="PCS-based HOS", linestyle='-', color='r', marker='s', linewidth=1.5)
    formatter = FuncFormatter(to_percent)
    ax2.yaxis.set_major_formatter(formatter)
    ax2.yaxis.label.set_fontsize(15)
    ax2.set_ylabel("Ratio between variance and amplitude", labelpad=10)

    names = ["lag#0","lag#1","lag#2","lag#3","lag#4","lag#5","lag#6"]
    ax.xaxis.label.set_fontsize(15)
    ax.yaxis.label.set_fontsize(15)
    ax.set_xticks([k+0.5 for k in range(7)])
    #ax.set_xticklabels(names,rotation=30, rotation_mode="anchor", ha="right")
    ax.set_xticklabels(names)
    ax.set_ylabel("Amplitude of cumulants", labelpad=10)
    ax.set_xlabel("Lags of the model", labelpad=10)
    plt.tight_layout()
    plt.grid(axis='y')
    plt.legend(loc=9, fancybox=True, shadow=True)

    plt.savefig(filename, format='pdf')
    plt.show()