def plotAgainstGFP(self, extradataA = [], extradataG = [], intensity = [], seq = []):
        fig1 = pylab.figure(figsize = (25, 10))
        print len(self.GFP)
        for i in xrange(min(len(data.cat), 3)):
            print len(self.GFP[self.categories == i])
            vect = []
            pylab.subplot(1,3,i+1)
            #pylab.hist(self.GFP[self.categories == i], bins = 20, color = data.colors[i])
            pop = self.GFP[self.categories == i]
            pylab.plot(self.GFP[self.categories == i], self.angles[self.categories == i], data.colors[i]+'o', markersize = 8)#, label = data.cat[i])
            print "cat", i, "n pop", len(self.GFP[(self.categories == i) & (self.GFP > -np.log(12.5))])
            x = np.linspace(np.min(self.GFP[self.categories == i]), np.percentile(self.GFP[self.categories == i], 80),40)
            #fig1.canvas.mpl_connect('pick_event', onpick)
            for j in x:
                vect.append(np.median(self.angles[(self.GFP > j) & (self.categories == i)]))

            pylab.plot([-4.5, -0.5], [vect[0], vect[0]], data.colors[i], label = "mediane de la population entiere", linewidth = 5)
            print vect[0], vect[np.argmax(x > -np.log(12.5))]
            pylab.plot([-np.log(12.5), -0.5], [vect[np.argmax(x > -np.log(12.5))] for k in  [0,1]], data.colors[i], label = "mediane de la population de droite", linewidth = 5, ls = '--')
            pylab.axvline(x = -np.log(12.5), color = 'm', ls = '--', linewidth = 3)
            pylab.xlim([-4.5, -0.5])
            pylab.legend(loc = 2, prop = {'size':17})

            pylab.title(data.cat[i].split(',')[0], fontsize = 24)
            pylab.xlabel('score GFP', fontsize = 20)
            pylab.ylabel('Angle (degre)', fontsize = 20)
            pylab.tick_params(axis='both', which='major', labelsize=20)
            pylab.ylim([-5, 105])
            ##pylab.xscale('log')
        pylab.show()

def plotDirections(aabb=(),mask=0,bins=20,numHist=True,noShow=False,sphSph=False):
	"""Plot 3 histograms for distribution of interaction directions, in yz,xz and xy planes and
	(optional but default) histogram of number of interactions per body. If sphSph only sphere-sphere interactions are considered for the 3 directions histograms.

	:returns: If *noShow* is ``False``, displays the figure and returns nothing. If *noShow*, the figure object is returned without being displayed (works the same way as :yref:`yade.plot.plot`).
	"""
	import pylab,math
	from yade import utils
	for axis in [0,1,2]:
		d=utils.interactionAnglesHistogram(axis,mask=mask,bins=bins,aabb=aabb,sphSph=sphSph)
		fc=[0,0,0]; fc[axis]=1.
		subp=pylab.subplot(220+axis+1,polar=True);
		# 1.1 makes small gaps between values (but the column is a bit decentered)
		pylab.bar(d[0],d[1],width=math.pi/(1.1*bins),fc=fc,alpha=.7,label=['yz','xz','xy'][axis])
		#pylab.title(['yz','xz','xy'][axis]+' plane')
		pylab.text(.5,.25,['yz','xz','xy'][axis],horizontalalignment='center',verticalalignment='center',transform=subp.transAxes,fontsize='xx-large')
	if numHist:
		pylab.subplot(224,polar=False)
		nums,counts=utils.bodyNumInteractionsHistogram(aabb if len(aabb)>0 else utils.aabbExtrema())
		avg=sum([nums[i]*counts[i] for i in range(len(nums))])/(1.*sum(counts))
		pylab.bar(nums,counts,fc=[1,1,0],alpha=.7,align='center')
		pylab.xlabel('Interactions per body (avg. %g)'%avg)
		pylab.axvline(x=avg,linewidth=3,color='r')
		pylab.ylabel('Body count')
	if noShow: return pylab.gcf()
	else:
		pylab.ion()
		pylab.show()

 def imshow_box(f,im, x,y,s):
     '''imshow_box(f,im, x,y,s)
     f: figure
     im: image
     x: center coordinate for box
     y: center coord
     s: box shape, (width, height)
     '''
     global coord
     P.figure(f.number)
     P.clf();
     P.imshow(im);
     P.axhline(y-s[1]/2.)
     P.axhline(y+s[1]/2.)
     P.axvline(x-s[0]/2.)
     P.axvline(x+s[0]/2.)
     xy=crop(m,s,y,x)
     coord=(0.5*(xy[2]+xy[3]), 0.5*(xy[0]+xy[1]))
     P.title(str('x: %d y: %d' % (x,y)));        
     P.figure(999);
     P.imshow(master[xy[0]:xy[1],xy[2]:xy[3]])
     P.title('Master');
     P.figure(998);
     df=(master[xy[0]:xy[1],xy[2]:xy[3]]-slave)
     P.imshow(np.abs(df))
     P.title(str('RMS: %0.6f' % np.sqrt((df**2.).mean()) ));        

def _show_rates(rate, wo, wt, attenuator, tau_NP, tau_P):
    import pylab

    #pylab.figure()
    pylab.errorbar(rate, wt[0], yerr=wt[1], fmt='g.', label='attenuated')
    pylab.errorbar(rate, wo[0], yerr=wo[1], fmt='b.', label='unattenuated')

    pylab.xscale('log')
    pylab.yscale('log')
    pylab.xlabel('incident rate (counts/second)')
    pylab.ylabel('observed rate (counts/second)')
    pylab.legend(loc='best')
    pylab.grid(True)
    pylab.plot(rate, rate/attenuator, 'g-', label='target')
    pylab.plot(rate, rate, 'b-', label='target')

    Ipeak, Rpeak = peak_rate(tau_NP=tau_NP, tau_P=tau_P)
    if rate[0] <= Ipeak <= rate[-1]:
        pylab.axvline(x=Ipeak, ls='--', c='b')
        pylab.text(x=Ipeak, y=0.05, s=' %g'%Ipeak,
                   ha='left', va='bottom',
                   transform=pylab.gca().get_xaxis_transform())
    if False:
        pylab.axhline(y=Rpeak, ls='--', c='b')
        pylab.text(y=Rpeak, x=0.05, s=' %g\n'%Rpeak,
                   ha='left', va='bottom',
                   transform=pylab.gca().get_yaxis_transform())

def plotLDDecaySpaceTime2d(ld0=1):
    T=np.arange(0,1500+1,500)
    L=1e6+1
    pos=500000
    r=2*1e-8
    s=0.01;x0=0.005
    # s=0.05;x0=1e-3

    positions=np.arange(0,L,1000)
    dist=abs(positions - pos)
    def getColor(n):
        color=['k','red','blue','g','m','c','coral']
        return color[:n]
    plt.figure(figsize=(25,12))
    for t,color in zip(T,getColor(len(T))):
        LD(t,ld0,s,x0,r,dist,positions).plot(ax=plt.gca(),linewidth=2, color=color,label='t={}'.format(t))
    for t,color in zip(T,getColor(len(T))):
        if not t :continue
        pd.Series(ld0*np.exp(-r*t*(dist)),index=positions).plot(ax=plt.gca(),style='--',color=color,linewidth=2)
    plt.legend(map(lambda x: 't={}'.format(x),T),loc='best');

    plt.axvline(500000,color='k',linewidth=2)
    plt.gca().axvspan(475000, 525000, alpha=0.25, color='black')
    plt.grid();plt.ylim([-0.05,1.1]);plt.xlabel('Position');plt.ylabel('LD to Position 500K');plt.title('Space-Time Decay of LD under Neutral Evolution');
    # plt.savefig(Simulation.paperFiguresPath+'LDDecay2d')
    plt.show()

def simFlips(numFlips, numTrials):      # performs and displays the simulation result
    diffs = []                          # diffs to know if there was a fair Trial. It has the absolute differences of heads and tails in each trial
    for i in xrange(0, numTrials):           
        heads, tails = flipTrial(numFlips)
        diffs.append(abs(heads - tails))
                
    diffs = pylab.array(diffs)          # create an array of diffs
    diffMean = sum(diffs)/len(diffs)    # average of absolute differences of heads and tails from each trial
    diffPercent = (diffs/float(numFlips)) * 100     # create an array of percentage of each diffs from its no. of flips.
    percentMean = sum(diffPercent)/len(diffPercent)     # create a percent mean of all diffPercents in the array
    
    pylab.hist(diffs)                   # displays the distribution of elements in diffs array
    pylab.axvline(diffMean, color = 'r', label = 'Mean')
    pylab.legend()
    titleString = str(numFlips) + ' Flips, ' + str(numTrials) + ' Trials'
    pylab.title(titleString)
    pylab.xlabel('Difference between heads and tails')
    pylab.ylabel('Number of Trials')
    
    pylab.figure()
    pylab.plot(diffPercent)
    pylab.axhline(percentMean, color = 'r', label = 'Mean')
    pylab.legend()
    pylab.title(titleString)
    pylab.xlabel('Trial Number')
    pylab.ylabel('Percent Difference between heads and tails')

def esat_comparison_plot(t=_np.linspace(173.15, 373.15, 20), std="Hyland_Wexler", percent=True, log=False):
    import pylab
    import brewer2mpl

    methods = list(esat(0, method="return"))
    methods.remove(std)
    print len(methods)
    pylab.rcParams.update({"axes.color_cycle": brewer2mpl.get_map("Paired", "qualitative", 12).mpl_colors})
    y = esat(t, method=std, info=False)
    i = 0
    style = "-"
    for im in methods:
        if i > 11:
            style = "--"
        print im
        if percent:
            pylab.plot(t, 100.0 * (esat(t, method=im, info=False) - y) / y, lw=2, ls=style)
        else:
            pylab.plot(t, esat(t, method=im, info=False) - y, lw=2, ls=style)
        i += 1
    pylab.legend(methods, loc="upper right", fontsize=8)
    pylab.xlabel("Temperature [K]")
    if percent:
        # pylab.semilogy()
        pylab.ylabel("Water Vapor Pressure Difference [%]")
    else:
        pylab.ylabel("Water Vapor Pressure [Pa]")
    pylab.title("Comparison of Water Vapor Calculations Ref:" + std)
    pylab.xlim(_np.round(t[0]), _np.round(t[-1]))
    pylab.grid()
    pylab.axvline(x=273.15, color="k")
    if log:
        pylab.yscale("log")

def Cross(x0=0.0, y0=0.0, clr='black', ls='dashed', lw=1, zorder=0):
    """
    Draw cross through zero
    =======================
    """
    axvline(x0, color=clr, linestyle=ls, linewidth=lw, zorder=zorder)
    axhline(y0, color=clr, linestyle=ls, linewidth=lw, zorder=zorder)

def test():
    if 0:
        from pandas import DataFrame
        X = np.linspace(0.01, 1.0, 10)
        Y = np.log(X)
        Y -= Y.min()
        Y /= Y.max()
        Y *= 0.95

        df = DataFrame({'X': X, 'Y': Y})
        P = Pareto(df, 'X', 'Y')

        data = []
        for val in np.linspace(0,1,15):
            data.append(dict(val=val, x=P.lookup_x(val), y=P.lookup_y(val)))
            pl.axvline(val, alpha=.5)
            pl.axhline(val, alpha=.5)
        dd = DataFrame(data)
        pl.scatter(dd.y, dd.val, lw=0, c='r')
        pl.scatter(dd.val, dd.x, lw=0, c='g')
        print dd

        #P.scatter(c='r', lw=0)
        P.show_frontier(c='r', lw=4)
        pl.show()

    X,Y = np.random.normal(0,1,size=(2, 30))

    for maxX in [0,1]:
        for maxY in [0,1]:
            pl.figure()
            pl.title('max x: %s, max y: %s' % (maxX, maxY))
            pl.scatter(X,Y,lw=0)
            show_frontier(X, Y, maxX=maxX, maxY=maxY)
            pl.show()

    def dofit(xs, slambdas, alpha, sinbeta, gamma, delta, band, y):
        xs = np.array(xs)
        ok = np.isfinite(slambdas)
        lsf = Fit.do_fit_wavelengths(xs[ok], lines[ok], alpha, sinbeta, 
                gamma, delta, band, y, error=slambdas[ok])
        (order, pixel_y, alpha, sinbeta, gamma, delta) = lsf.params

        print "order alpha    sinbeta   gamma  delta"
        print "%1.0i %5.7f %3.5f %3.2e %5.2f" % (order, alpha, sinbeta, 
                gamma, delta)

        ll = Fit.wavelength_model(lsf.params, pix)
        if DRAW:
            pl.figure(3)
            pl.clf()
            pl.plot(ll, spec)
            pl.title("Pixel %4.4f" % pos)
            for lam in lines:
                pl.axvline(lam, color='r')

            pl.draw()
    
        return [np.abs((
            Fit.wavelength_model(lsf.params, xs[ok]) - lines[ok]))*1e4, 
                lsf.params]

    def test_band_unpolarized(self):
        """polarizedのPROCAR用"""
        # import seaborn
        path = os.path.join(self.path, 'unpolarized', 'C_P63mmmc')
        dos = collect_vasp.Doscar(os.path.join(path, 'DOSCAR_polarized'))
        dos.get_data()
        e_fermi = dos.fermi_energy  # DOSCARからのEf
        band = collect_vasp.Procar(os.path.join(path, 'PROCAR_band'),
                                   is_polarized=False)
        band['energy'] = band['energy'] - e_fermi
        band.output_keys = ['kpoint_id', 'energy', "spin"]
        kp_label = band.get_turning_kpoints_pmg(os.path.join(path, 'KPOINTS_band'))

        # plot
        plt = pylab.figure(figsize=(7, 7/1.618))
        ax1 = plt.add_subplot(111)
        ax1.set_ylabel("Energy from $E_F$(meV)")
        ax1.set_xlabel("BZ direction")

        # 枠
        for i in kp_label[0]:
            pylab.axvline(x=i, ls=':', color='gray')
        pylab.axhline(y=0, ls=':', color='gray')

        ax1.scatter(band['kpoint_id'], band['energy'], s=3, c='blue',
                    linewidths=0)
        ax1.set_xlim(0, band['kpoint_id'][-1])

        pylab.xticks(kp_label[0], kp_label[1])
        pylab.show()

def simpleExample():
	# Create a signal
	N = 256
	signal = numpy.zeros(N,numpy.complex128)
	# Make our signal 1 in the middle, zero elsewhere
	signal[N/2] = 1.0
	# plot our signal
	pylab.figure()
	pylab.plot(abs(signal))
	# Do the GFT on the signal
	SIGNAL = gft.gft1d(signal,'gaussian')
	# plot the magnitude of the signal
	pylab.figure()
	pylab.plot(abs(SIGNAL),'b')
	
	# get the partitions
	partitions = gft.partitions(N)
	# for each partition, draw a line on the graph.  Since the partitions only indicate the
	# positive frequencies, we need to draw partitions on both sides of the DC
	for i in partitions[0:len(partitions)/2]:
		pylab.axvline((N/2+i),color='r',alpha=0.2)
		pylab.axvline((N/2-i),color='r',alpha=0.2)
		
	# finally, interpolate the GFT spectrum and plot a spectrogram
	pylab.figure()
	pylab.imshow(abs(gft.gft1dInterpolateNN(SIGNAL)))

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 measure_tae():
   print "Measuring initial perception of all orientations..."
   before=test_all_orientations(0.0,0.0)
   pylab.figure(figsize=(5,5))
   vectorplot(degrees(before.keys()),  degrees(before.keys()),style="--") # add a dashed reference line
   vectorplot(degrees(before.values()),degrees(before.keys()),\
             title="Initial perceived values for each orientation")

   print "Adapting to pi/2 gaussian at the center of retina for 90 iterations..."
   for p in ["LateralExcitatory","LateralInhibitory","LGNOnAfferent","LGNOffAfferent"]:
      # Value is just an approximate match to bednar:nc00; not calculated directly
      topo.sim["V1"].projections(p).learning_rate = 0.005

   inputs = [pattern.Gaussian(x = 0.0, y = 0.0, orientation = pi/2.0,
                     size=0.088388, aspect_ratio=4.66667, scale=1.0)]
   topo.sim['Retina'].input_generator.generators = inputs
   topo.sim.run(90)


   print "Measuring adapted perception of all orientations..."
   after=test_all_orientations(0.0,0.0)
   before_vals = array(before.values())
   after_vals  = array(after.values())
   diff_vals   = before_vals-after_vals # Sign flipped to match conventions

   pylab.figure(figsize=(5,5))
   pylab.axvline(90.0)
   pylab.axhline(0.0)
   vectorplot(wrap(-90.0,90.0,degrees(diff_vals)),degrees(before.keys()),\
             title="Difference from initial perceived value for each orientation")

def plot(kmv):
    py.scatter([d / float(2 ** 32 - 1)
                for d in kmv.data[:-1]], [0] * (len(kmv.data) - 1), alpha=0.25)
    py.axvline(x=(kmv.data[-2] / float(2 ** 32 - 1)), c='r')
    py.gca().get_yaxis().set_visible(False)
    py.gca().get_xaxis().set_ticklabels([])
    py.gca().get_xaxis().set_ticks([x / 10. for x in xrange(11)])

    def draw_cell_size_distribution_diagram(self, diagram_file_name, cell_list):

        cell_distribution_hashmap = {}
        nr_of_points_list = []
        for cell in cell_list:
            value = cell_distribution_hashmap.get(cell.get_nr_of_points(), 0)
            cell_distribution_hashmap[cell.get_nr_of_points()] = value + 1
            nr_of_points_list.append(cell.get_nr_of_points())

        index_max = cell_distribution_hashmap.keys()[-1] + 3
        index_min = 0

        value_list = []
        for index in range(index_min, index_max):
            value_list.append(cell_distribution_hashmap.get(index, 0))


        median_value = numpy.median(nr_of_points_list)
        print '[i] calculated median value: %f' % median_value

        mean_value =  float(sum(value_list)) / len(value_list)
        print '[i] calculated mean value: %f' % mean_value

        pylab.plot(value_list, 'r',  label='cell size distribution')
        pylab.axvline(x=median_value, label='median')

        pylab.xlabel('cell size [pixel]')
        pylab.ylabel('cell count [cells]')
        pylab.title('cell size distribution')
        pylab.grid(True)
        pylab.legend(loc= 'best')

        pylab.savefig(self.log_path + '/' + diagram_file_name)

	def plot_statistic_vs_burnin(self, param, points=20, mean=True, error=False, showpoints=True):
		samples = self.samples[parameters[param]]
		burn = np.linspace(0, len(samples), points)[:-1].astype(int)
		means = []
		stds=[]
		for b in burn:
			if mean:
				means.append( np.mean(samples[b:]) )

			if error:
				stds.append(np.std(samples[b:]))


		if mean:
			plt.plot(burn, means, "-", color="purple", lw=2.0, label="mean")
			if showpoints:
				plt.plot(burn, means, "o", color="purple")
			plt.ylabel("mean %s"%labels[param])

		if error:
			plt.plot(burn, stds, "-", color="midnightblue", lw=2.0, label="error")
			if showpoints:
				plt.plot(burn, stds, "o", color="midnightblue")
			plt.ylabel("error %s"%labels[param])

		if error and mean:
			plt.ylabel("error or mean %s"%labels[param])
			plt.legend(loc="upper left")


		plt.xlabel("points discarded")
		

		plt.axvline(len(samples), color="k", lw=2)

 def plot(self, gs):
     unit_len = self.show_len * 1. / 5.
     if self.s.now - self.show_len < 0:
         return 
         
     price = self.price[0][self.s.now - self.show_len : self.s.now]
     profile_range = [price.min(), price.max() + 1]
     floor, ceil = profile_range[0] - 1, profile_range[1] + 1
         
     d = self.output(3, profile_range)
     
     ax = plt.subplot(gs)
     plt.plot(price)
     day_begin = np.where(self.s.history['time_in_ticks'][self.s.now - self.show_len : self.s.now] == 0)[0]
     for x in day_begin:
         plt.axvline(x, color='r', linestyle=':')
     y = self.smoothed_pivot_profile[floor : ceil]
     plt.barh(np.arange(floor, ceil) - 0.5, y * unit_len, 1.0, label=self.name,
              alpha=0.2, color='r', edgecolor='none')
     
     last_price = int(get(self.price))
     support = last_price + int(round((d['S_offset']) * self.volatility))
     resistance = last_price + int(round((d['R_offset']) * self.volatility))
     highlighted = [support, resistance]
     plt.barh(np.array(highlighted) - 0.5, self.smoothed_pivot_profile[highlighted] * unit_len, 1.0,
              alpha=1.0, color='r', edgecolor='none')
     ax.set_xticks(np.arange(0, self.show_len * 1.22, unit_len))
     ax.xaxis.grid(b=True, linestyle='--')
     ax.yaxis.grid(b=False)
     plt.legend(loc='upper right')
     return ax

    def plot(self):
        f = pylab.figure(figsize=(8,4))
        co = [] #colors container
        for zScore, r in itertools.izip(self.zScores, self.log2Ratio):
            if zScore < self.pCut:
                if r > 0:
                    co.append(Colors().greenColor)
                elif r < 0:
                    co.append(Colors().redColor)
                else:
                    raise Exception
            else:
                co.append(Colors().blueColor)

        #print "Probability this is from a normal distribution: %.3e" %stats.normaltest(self.log2Ratio)[1]
        ax = f.add_subplot(121)
        pylab.axvline(self.meanLog2Ratio, color=Colors().redColor)
        pylab.axvspan(self.meanLog2Ratio-(2*self.stdLog2Ratio), 
                      self.meanLog2Ratio+(2*self.stdLog2Ratio), color=Colors().blueColor, alpha=0.2)
        his = pylab.hist(self.log2Ratio, bins=50, color=Colors().blueColor)
        pylab.xlabel("log2 Ratio %s/%s" %(self.sampleNames[1], self.sampleNames[0]))
        pylab.ylabel("Frequency")
        
        ax = f.add_subplot(122, aspect='equal')
        pylab.scatter(self.genes1, self.genes2, c=co, alpha=0.5)        
        pylab.ylabel("%s RPKM" %self.sampleNames[1])
        pylab.xlabel("%s RPKM" %self.sampleNames[0])
        pylab.yscale('log')
        pylab.xscale('log')
        pylab.tight_layout()

    def main(self):
        global weights, densities, weighted_densities
        plt.figure()

        cluster = clusters.SingleStation()
        self.station = cluster.stations[0]

        R = np.linspace(0, 100, 100)
        densities = []
        weights = []
        for E in np.linspace(1e13, 1e17, 10000):
            relative_flux = E ** -2.7
            Ne = 10 ** (np.log10(E) - 15 + 4.8)
            self.ldf = KascadeLdf(Ne)
            min_dens = self.calculate_minimum_density_for_station_at_R(R)

            weights.append(relative_flux)
            densities.append(min_dens)
        weights = np.array(weights)
        densities = np.array(densities).T

        weighted_densities = (np.sum(weights * densities, axis=1) /
                              np.sum(weights))
        plt.plot(R, weighted_densities)
        plt.yscale('log')
        plt.ylabel("Min. density [m^{-2}]")
        plt.xlabel("Core distance [m]")
        plt.axvline(5.77)
        plt.show()