def transmission(path,g,src,dest):
 global disp_count
 global bar_colors
 global edge_colors
 global node_colors
 k=0
 j=0
 list_of_edges = g.edges() 
 for node in path :
  k=path.index(node)
  disp_count = disp_count + 1
  if k != (len(path)-1):
   k=path[k+1]
   j=list_of_edges.index((node,k))
   initialize_edge_colors(j)
   #ec[disp_count].remove(-3000) 
   pylab.subplot(121) 
   nx.draw_networkx(g,pos = nx.circular_layout(g),node_color= node_colors,edge_color = edge_colors)
   pylab.annotate("Source",node_positions[src])
   pylab.annotate("Destination",node_positions[dest])
   pylab.title("Transmission")

   he=initializeEnergies(disp_count)
   print he
   pylab.subplot(122)
   pylab.bar(left=[1,2,3,4,5],height=[300,300,300,300,300],width=0.5,color = ['w','w','w','w','w'],linewidth=0) 
   pylab.bar(left=[1,2,3,4,5],height=initializeEnergies(disp_count),width=0.5,color = 'b') 
   pylab.title("Node energies")
   #pylab.legend(["already passed throgh","passing" , "yet to pass"])
   pylab.xlabel('Node number')
   pylab.ylabel('Energy') 
   pylab.suptitle('Leach Protocol', fontsize=12)
   pylab.pause(2)
  else :
     return

def group_plots(ylist, ncols, x = None,
		titles = None,
		suptitle = None,
		ylabels = None,
		figsize = None,
		sameyscale = True,
                order='C',
		imkw={}):
    import pylab as pl
    nrows = np.ceil(len(ylist)/float(ncols))
    figsize = ifnot(figsize, (2*ncols,2*nrows))
    fh, axs = pl.subplots(int(nrows), int(ncols),
                          sharex=True,
                          sharey=bool(sameyscale),
                          figsize=figsize)
    ymin,ymax = data_range(ylist)
    axlist = axs.ravel(order=order)
    for i,f in enumerate(ylist):
	x1 = ifnot(x, range(len(f)))
        _im = axlist[i].plot(x1,f,**imkw)
	if titles is not None:
            pl.setp(axlist[i], title = titles[i])
	if ylabels is not None:
            pl.setp(axlist[i], ylabel=ylabels[i])
    if suptitle:
        pl.suptitle(suptitle)
    return

    def log_posterior(self,theta):

        model_g1 , model_g2, limit_mask , _ , _ = self.draw_model(theta)

        likelihood = self.log_likelihood(model_g1,model_g2,limit_mask)
        prior = self.log_prior(theta)
        if not np.isfinite(prior):
            posterior = -np.inf
        else:
            # use no info from prior for now
            posterior = likelihood 

        if logger.level == logging.DEBUG:
            n_progress = 10
        elif logger.level == logging.INFO:
            n_progress = 1000
        if self.n_model_evals % n_progress == 0:

            logger.info('%7d post=% 2.8e like=% 2.8e prior=% 2.4e M200=% 6.3e ' % (self.n_model_evals,posterior,likelihood,prior,theta[0]))

        if np.isnan(posterior):
            import pdb; pdb.set_trace()

        if self.save_all_models:

            self.plot_residual_g1g2(model_g1,model_g2,limit_mask)

            pl.suptitle('model post=% 10.8e M200=%5.2e' % (posterior,theta[0]) )
            filename_fig = 'models/res2.%04d.png' % self.n_model_evals
            pl.savefig(filename_fig)
            logger.debug('saved %s' % filename_fig)
            pl.close()


        return posterior

def plot_weights(init_som, final_som, title=['SOM init', 'SOM final'], dim_lab=None):
	'''
	Function to plot neural weights before and after the training for each dimension.
	'''	
	
	assert init_som.shape == final_som.shape
	
	n, d = init_som.shape
	width = np.int(np.sqrt(n))
	
	if dim_lab is None:
		dim_lab = ['w' + str(i) for i in xrange(d)]

	fig = plt.figure()
	
	for lab, i in zip(dim_lab, xrange(d)):
			
		# plot weights before training
		plt.suptitle(title[0], fontsize = 14)
		ax = fig.add_subplot(2, d, i+1)
		img = init_som[:, i].reshape(width, width)
		ax.imshow(img, interpolation='nearest')
		plt.title(lab)

		# same weights after training

		ax = fig.add_subplot(2, d, (i+1) + d)
		if i==int(d/2.): plt.title(title[1])
		img_f = final_som[:, i].reshape(width, width)
		ax.imshow(img_f, interpolation='nearest')

    def plot_bases(self, autoscale=True, stampsize=None):
        import pylab as plt

        N = len(self.psfbases)
        cols = int(np.ceil(np.sqrt(N)))
        rows = int(np.ceil(N / float(cols)))
        plt.clf()
        plt.subplots_adjust(hspace=0, wspace=0)

        cut = 0
        if stampsize is not None:
            H, W = self.shape
            assert H == W
            cut = max(0, (H - stampsize) / 2)

        ima = dict(interpolation="nearest", origin="lower")
        if autoscale:
            mx = self.psfbases.max()
            ima.update(vmin=-mx, vmax=mx)
        nil, xpows, ypows = self.polynomials(0.0, 0.0, powers=True)
        for i, (xp, yp, b) in enumerate(zip(xpows, ypows, self.psfbases)):
            plt.subplot(rows, cols, i + 1)

            if cut > 0:
                b = b[cut:-cut, cut:-cut]
            if autoscale:
                plt.imshow(b, **ima)
            else:
                mx = np.abs(b).max()
                plt.imshow(b, vmin=-mx, vmax=mx, **ima)
            plt.xticks([])
            plt.yticks([])
            plt.title("x^%i y^%i" % (xp, yp))
        plt.suptitle("PsfEx eigen-bases")

def get_cav(c_mat,nchan,scaling=False):
        #Compute Cav by taking diags, averaging and reforming the matrix
        diags =[c_mat.diagonal(count) for count in xrange(nchan-1, -nchan,-1)]
        cav=n.zeros_like(c_mat)

        for count,count_chan in enumerate(range(nchan-1,-nchan,-1)):
                cav += n.diagflat( n.mean(diags[count]).repeat(len(diags[count])), count_chan)

        if scaling:
            temp=cav.copy()
            for count in xrange(nchan):
                cav[count,:] *= n.sqrt(c_mat[count,count]/temp[count,count])
                cav[:,count] *= n.sqrt(c_mat[count,count]/temp[count,count])

        if not n.allclose(cav.T.conj(),cav):
            print 'Cav is not Hermitian'
            print 'bl'

        if PLOT and False:
            p.subplot(131); capo.arp.waterfall(c_mat,mode='real',mx=3,drng=6);
            p.title('C')
            p.subplot(132); capo.arp.waterfall(cav,mode='real',mx=3,drng=6); p.colorbar()
            p.title('Cav')
            p.subplot(133); capo.arp.waterfall(diff,mode='real',mx=500,drng=500); p.colorbar()
            p.title('Abs Difference')
            p.suptitle('%d_%d'%a.miriad.bl2ij(bl))
            p.show()

        return cav

def _show_plots(target, fitted, wt, wo, corrected):
    tau_NP, tau_P, attenuator, rate = target
    tau_NP_f, tau_P_f, attenuation_f, rate_f = fitted

    # Plot the results
    sim_pars = (r'Sim $\tau_{NP}=%g\,{\rm %s}$,  $\tau_{P}=%g\,{\rm %s}$,  ${\rm attenuator}=%g$'
               )%(tau_NP, DEADTIME_UNITS, tau_P, DEADTIME_UNITS, attenuator)
    fit_pars = (r'Fit $\tau_{NP}=%s$,  $\tau_P=%s$,  ${\rm attenuator}=%.2f$'
               )%(
                   ("%.2f"%tau_NP_f[0] if np.inf > tau_NP_f[1] > 0 else "-"),
                   ("%.2f"%tau_P_f[0] if np.inf > tau_P_f[1] > 0 else "-"),
                   1./attenuation_f[0],
               )
    title = '\n'.join((sim_pars, fit_pars))
    import pylab
    pylab.subplot(211)
    #pylab.errorbar(rate, rate_f[0], yerr=rate_f[1], fmt='c.', label='fitted rate')
    #mincident = np.linspace(rate[0], rate[-1], 400)
    #munattenuated = expected_rate(mincident, tau_NP_f[0], tau_P_f[0])
    #mattenuated = expected_rate(mincident/attenuator, tau_NP_f[0], tau_P_f[0])
    #minc = np.hstack((mincident, 0., mincident))
    #mobs = np.hstack((munattenuated, np.NaN, mattenuated))
    #pylab.plot(minc, mobs, 'c-', label='expected rate')
    pylab.errorbar(rate, uval(corrected), yerr=udev(corrected), fmt='r.', label='corrected rate')
    _show_rates(rate, wo, wt, attenuator, tau_NP_f[0], tau_P_f[0])
    pylab.subplot(212)
    _show_droop(rate, wo, wt, attenuator)
    pylab.suptitle(title)

    #pylab.figure(); _show_inversion(wo, tau_P_f, tau_NP_f)
    pylab.show()

def compare_expvaluecollections(coll1, coll2, show=True, **kwargs):
    """
    Plot all subsystems of two ExpectationValueCollections.

    *Arguments*
        * *coll1*
            First :class:`pycppqed.expvalues.ExpectationValue.Collection`.

        * *coll2*
            Second :class:`pycppqed.expvalues.ExpectationValue.Collection`.

        * *show* (optional):
            If True pylab.show() is called finally. This means a plotting
            window will pop up automatically. (Default is True)

        * Any other arguments that the pylab plotting command can use.

    This function allows a fast comparison between two sets of expectation
    values that were obtained by different calculations.
    """
    import pylab
    s1 = coll1.subsystems
    s2 = coll2.subsystems
    assert len(s1) == len(s2)
    for i in range(len(s1)):
        pylab.figure()
        _compare_expvaluesubsystems(s1.values()[i], s2.values()[i],
                                        show=False, **kwargs)
        title = "%s vs. %s" % (s1.keys()[i], s2.keys()[i])
        if hasattr(pylab, "suptitle"): # For old versions not available.
            pylab.suptitle(title)
            pylab.gcf().canvas.set_window_title(title)
    if show:
        pylab.show()

 def plot_prototypes(self):
         """ plots each of the components as a prototype (sum of fitted b-splines) and returns a dictionary of figures """
         figs = {}
         for h in np.unique(self.design.header):
                 fig = pl.figure()
                 splines = self.design.get(h)
                 if 'rate' in h:
                         pl.plot(np.sum(self.design.get(h)[:self.design.trial_length] * self.beta[self.design.getIndex(h)],1))
                         pl.title(h)
                 elif len(splines.shape) == 1 or (splines.shape[0] == 1):
                         pl.plot(np.sum(self.design.get(h) * self.beta[self.design.getIndex(h)],1),'o')
                         pl.title(h)
                 elif len(splines.shape) == 2:
                         pl.plot(np.sum(self.design.get(h) * self.beta[self.design.getIndex(h)],1))
                         pl.title(h)
                 elif len(splines.shape) == 3:
                         slices = np.zeros(splines.shape)
                         for (i, ind) in zip(range(splines.shape[0]),self.design.getIndex(h)):
                                 slices[i,:,:] = splines[i,:,:] * self.beta[ind]
                         pl.imshow(slices.sum(axis=0),cmap='jet')
                         figs[h + '_sum'] = fig
                         fig = pl.figure()
                         for i in range(len(slices)):
                                 pl.subplot(np.ceil(np.sqrt(slices.shape[0])),np.ceil(np.sqrt(slices.shape[0])),i+1)
                                 pl.imshow(slices[i],vmin=np.percentile(slices,1),vmax=np.percentile(slices,99),cmap='jet')
                         pl.suptitle(h)
                         figs[h] = fig
                 else:
                         pl.plot(np.sum(self.design.get(h) * self.beta[self.design.getIndex(h)],1))
                         pl.title(h)
                 figs[h] = fig
         return figs

def ks_gof(mcmc, samples, format="png"):
    """Runs ks_2samp test and plots Observed/Expected vs. Simulated/Expected"""
    size = len(samples)
    #Check for bedrock data
    for sample in samples:
        if isinstance(sample, BedrockSample):
            size = len(samples)-1
    #Set size and create grid
    if size == 1:
        fig = plt.figure(figsize=(10, 10))
    else:
        fig = plt.figure(figsize=(6, 10))
    grid = ag.axes_grid.Grid(fig, 111, nrows_ncols = (size, 1), axes_pad = 1, share_x=False, share_y=False, label_mode = "all" )
    j = 0 
    for sample in samples:
        if isinstance(sample, DetritalSample):
            #obs = mcmc.get_node("ObsAge_" + sample.name)
            sim = mcmc.get_node("SimAge_" + sample.name)
            exp = mcmc.get_node("ExpAge_" + sample.name)
            d_simExp=[]
            d_obsExp=[]
            #import ipdb; ipdb.set_trace()
            for i in range(len(exp.trace()[:,-1])):
                D, P = sp.stats.ks_2samp(mcmc.trace(exp)[i,-1], mcmc.trace(sim)[i,-1])
                d_simExp.append(D)
                D, P = sp.stats.ks_2samp(mcmc.trace(exp)[i,-1], sample.ages)
                d_obsExp.append(D)
                
            #The test statistics generated from ks_2samp plot as a grid. The following adds
            #random uniform noise for a better visual representation on the plot.
            noise=(0.5/len(sample.ages))
            for i in range(len(d_simExp)):
                d_simExp[i] = d_simExp[i] + np.random.uniform(-noise, noise, 1)
                d_obsExp[i] = d_obsExp[i] + np.random.uniform(-noise, noise, 1)

            #Calculate p-value (The proportion of test statistics above the y=x line)
            count=0
            for i in range(len(d_simExp)):
                if (d_simExp[i]>d_obsExp[i]):
                    count=count+1
            count=float(count)
            p_value=(count/len(d_simExp))
    
            #Plot
            grid[j].scatter(d_obsExp, d_simExp, color='gray', edgecolors='black')
            grid[j].set_xlabel('Observed/Expected', fontsize='small')
            grid[j].set_ylabel('Simulated/Expected', fontsize='small')
            label = sample.name + ", " + "p-value="+"%f"%p_value
            grid[j].set_title(label, fontsize='medium')

            #Add a y=x line to plot
            if (max(d_simExp)>=max(d_obsExp)):
                grid[j].plot([0,max(d_simExp)],[0,max(d_simExp)], color='black')
            else:
                grid[j].plot([0,max(d_obsExp)],[0,max(d_obsExp)], color='black')
            j+=1
            
            
    plt.suptitle('Kolmogorov-Smirnov Statistics', fontsize='large')
    fig.savefig("KS_test."+format)

def expvaluecollection(evc, show=True, **kwargs):
    """
    Visualize a :class:`pycppqed.expvalues.ExpectationValueCollection`.

    *Usage*
        >>> import numpy as np
        >>> T = np.linspace(0,10)
        >>> d = (np.sin(T), np.cos(T))
        >>> titles = ("<x>", "<y>")
        >>> evc = pycppqed.expvalues.ExpectationValueCollection(d, T, titles)
        >>> evc.plot()

    *Arguments*
        * *evc*
            A :class:`pycppqed.expvalues.ExpectationValueCollection`.

        * *show* (optional):
            If True pylab.show() is called finally. This means a plotting
            window will pop up automatically. (Default is True)

        * Any other arguments that the pylab plotting command can use.
    """
    if evc.subsystems:
        import pylab
        for sysname, data in evc.subsystems.iteritems():
            pylab.figure()
            _expvalues(data.evtrajectories, show=False, **kwargs)
            if hasattr(pylab, "suptitle"): # For old versions not available.
                pylab.suptitle(sysname)
                pylab.gcf().canvas.set_window_title(sysname)
        if show:
            pylab.show()
    else:
        titles = ["(%s) %s" % (i, title) for i, title in enumerate(evc.titles)]
        _expvalues(evc.evtrajectories, titles, show=show, **kwargs)

	def NeutralLinguagram(self, M, savename, start=1):
		fakeX = []
		for i in range(len(M)):
			xs = []
			for j in range(1,33):
				xs.append(j)
			fakeX.append(xs)
		
		x1 = array(fakeX)
		y1 = array(M)
		Z = []
		for i in range(start, (len(M)+start)):
			zs = []
			for j in range(32):
				zs.append(i)
			Z.append(zs)
		z1 = array(Z)
		
		fig = p.figure()
		ax = Axes3D(fig)
		ax.plot_surface(z1, -x1, y1, rstride=1, cstride=1, cmap=cm.jet)
		ax.view_init(90,-90)
		p.suptitle(savename[:-4])
		#p.show()
		
		p.savefig(savename, format = 'png')

    def createCoefGraph(data, nFig, lim, ymin):
        plt.figure(nFig)
        plt.suptitle('Coef')
        nBase = data.shape[0]
        nSubCols = nBase / 10
        if nSubCols > 0:
            nSubRows = nBase / nSubCols
        else:
            nSubRows = nBase 
            nSubCols = 	1
        # print data.shape

        # サンプリング周波数とシフト幅によって式を変える必要あり
        timeLine = [i * 1024 / 8000.0 for i in range(data.shape[1])]
        # print len(timeLine)
        for i in range(nBase):
            plt.subplot(nSubRows, nSubCols, i + 1)
            plt.tick_params(labelleft='off', labelbottom='off')
            # FIXME: Arguments of X
            # plt.plot(timeLine, data[i,:])
            if lim:
                plt.ylim(ymin=ymin)
            plt.plot(timeLine, data[i,:])
        # Beacuse I want to add lable in bottom, xlabel is declaration after loop.
        plt.tick_params(labelleft='off', labelbottom="on")
        plt.xlabel('time [ms]')

def plot_multiplot_histogram(data):
    xcol=4
    ycol=6
    for index,attrib in enumerate(attribList):
        byAttrib=get_byAttrib(data,attrib)
        P.suptitle('Attribute Histograms')
        ax = P.subplot(xcol,ycol,index+1)
        plot_histogram(byAttrib,attrib=attribDict[attrib],bool_labels=False)
        for item in ax.get_xticklabels():
            item.set_fontsize(0)
        for item in ax.get_yticklabels():
            item.set_fontsize(0)
        if index % ycol == 0:
            P.ylabel('Probability')
            for item in ax.get_yticklabels():
                item.set_fontsize(8)
        if index > (xcol-1)*ycol-1:
            P.xlabel('Rank')
            for item in ax.get_xticklabels():
                item.set_fontsize(8)
        P.xlim([1,24])
        P.ylim([0,0.25])
        ax.yaxis.label.set_size(10)
        ax.xaxis.label.set_size(10)
        if np.sum(byAttrib[0:7,1])>2*np.sum(byAttrib[16:23,1]):
            P.text(20,0.20,'+')
        elif np.sum(byAttrib[16:23,1])>2*np.sum(byAttrib[0:7,1]):
            P.text(20,0.20,'-')
    P.subplots_adjust(hspace=.50)
    P.subplots_adjust(wspace=.50)
    P.subplots_adjust(bottom=0.1)

def check_negative_offsets(test_arcs):
  d = 0.4
  # Offset inward then outward would normally erode sharp features, but we can use a positive offset to generate a shape with no sharp features
  outer = offset_arcs(test_arcs, d*1.5) # Ensure features big enough to not disappear if we inset/offset again
  inset = offset_arcs(outer, -d)
  reset = offset_arcs(inset, d)
  outer_area = circle_arc_area(outer)
  inset_area = circle_arc_area(inset)
  reset_area = circle_arc_area(reset)
  assert inset_area < outer_area # Offset by negative amount should reduce area
  if not fuzzy_arcs_equal(outer, reset, 2):
    if 0:
      import pylab
      area_error = abs(outer_area - reset_area)
      pylab.suptitle("Offset: %s, Area error:%s" % (sub_d, area_error))
      print "outer_area:",outer_area
      print "inset_area:",inset_area
      print "reset_area:",reset_area
      subplot_arcs(outer, 131, "outer", full=False)
      subplot_arcs(inset, 132, "inset", full=False)
      subplot_arcs(reset, 133, "reset", full=False)
      pylab.show()
    assert False
  # Check that a large negative offset leaves nothing
  empty_arcs = offset_arcs(test_arcs, -100.)
  assert len(empty_arcs) == 0

 def _barplot_engine(self, suptitle, numIter=None, saveas=None, show=True):
     '''
     code to support a barchart visualization of the results, with one
     bar chart per iteration.
     '''
     res = self.annotation_dat if numIter is None else self.annotation_dat[0:numIter]
     
     pl.figure(num=1)
         
     #AS is complete Annotation Set, max_count is the maximum number
     # of genes that appears in any single annotation entry
     (AS, max_count) = self._common_Y()
     
     for (i, dat, genelist) in res:
         if len(dat) < 1: continue
         pl.subplot(2, math.ceil( len(res)/2.0), i+1)
         title = "Iteration %d (%d Genes)"%((i+1), len(genelist))
         flag1=False if i in [0,5] else True
         
         self._plot_GO_bar_chart(i, title, annotation_set=AS, xlim=[0,max_count],
                               suppress_ylab=flag1, suppress_desc=False, desc_filter=True,
                               color="yellow", show=False)  
     pl.subplots_adjust(left=0.1, right=0.95, top=0.90, bottom=0.05, wspace=0.1, hspace=0.2)
     if not suptitle is None:
         pl.suptitle("Ontology Annotations per Iteration: %s"%suptitle)
     if not saveas is None:
         pl.savefig(saveas)   
     if show: pl.show()

def simulationDelayedTreatment():
    """
    Runs simulations and make histograms for problem 5.
    Runs multiple simulations to show the relationship between delayed treatment
    and patient outcome.
    Histograms of final total virus populations are displayed for delays of 300,
    150, 75, 0 timesteps (followed by an additional 150 timesteps of
    simulation).    
    """
    histBins = [i*50 for i in range(11)]
    delays = [0, 75, 150, 300]
    subplot = 1
    for d in delays:
        results = []
        for t in xrange(NUM_TRIALS):
            results.append(runSimulation(d))
        pylab.subplot(2, 2, subplot)
        subplot += 1
        pylab.hist(results, bins=histBins, label='delayed ' + str(d))
        pylab.xlim(0, 500)
        pylab.ylim(0, NUM_TRIALS)
        pylab.ylabel('number of patients')
        pylab.xlabel('total virus population')
        pylab.title(str(d) + ' time step delay')
        popStd = numpy.std(results)
        popMean = numpy.mean(results)
        ## print str(d)+' step delay standard deviation: '+str(popStd)
        ## print str(d)+' step mean: '+str(popMean)
        ## print str(d)+' step CV: '+str(popStd / popMean)
        ## print str(d) + ' step delay: ' + str(results)
    pylab.suptitle('Patient virus populations after 150 time steps when ' +\
                   'prescription\n' +\
                   'is applied after delays of 0, 75, 150, 300 time steps')
    pylab.show()    

def run_demo(with_plots=True):
    """
    Demonstrates how to use PyFMI for simulation of 
    Co-Simulation FMUs (version 1.0).
    """
    fmu_name = O.path.join(path_to_fmus,'bouncingBall.fmu')
    model = load_fmu(fmu_name)
    
    res = model.simulate(final_time=2.)
    
    # Retrieve the result for the variables
    h_res = res['h']
    v_res = res['v']
    t     = res['time']

    assert N.abs(res.final('h') - (0.0424044)) < 1e-2
    
    # Plot the solution
    if with_plots:
        # Plot the height
        fig = P.figure()
        P.clf()
        P.subplot(2,1,1)
        P.plot(t, h_res)
        P.ylabel('Height (m)')
        P.xlabel('Time (s)')
        # Plot the velocity
        P.subplot(2,1,2)
        P.plot(t, v_res)
        P.ylabel('Velocity (m/s)')
        P.xlabel('Time (s)')
        P.suptitle('FMI Bouncing Ball')
        P.show()

    def plotGeometry(self,SNP,PCAD,title):
        fig=plt.figure(figsize=self.figsize, dpi=self.dpi);plt.ioff()
        SNP.loc['Reference0']=np.zeros(SNP.shape[1],dtype=int)
        plt.subplot(1,2,1)
        D=pd.DataFrame(None,index=SNP.index,columns=SNP.index)
        for i in SNP.index:
            for j in SNP.index:
                D.loc[i,j]= sum(np.logical_xor(SNP.loc[i],SNP.loc[j]))
        D=D.astype(float)
        im=plt.imshow(D,interpolation='nearest',cmap='Reds')
        plt.gca().xaxis.tick_top()
        x=np.arange(D.index.shape[0])
        plt.yticks(x,map(lambda x: x.replace('mdio','') ,D.index.values))
        plt.xticks(x,map(lambda x: x.replace('mdio','') ,D.columns))
        plt.colorbar(im)
        plt.gca().tick_params(axis='both', which='major', labelsize=10)
        plt.title('Pairwise Hamming Distance',y=1.03)

        
        plt.subplot(1,2,0)
        D=PCAD.astype(float)
        im=plt.imshow(D,interpolation='nearest',cmap='Reds')
        plt.gca().xaxis.tick_top()
        x=np.arange(D.index.shape[0])
        plt.yticks(x,map(lambda x: x.replace('mdio','') ,D.index.values))
        plt.xticks(x,map(lambda x: x.replace('mdio','') ,D.columns))
        plt.colorbar(im)
        plt.gca().tick_params(axis='both', which='major', labelsize=10)
        plt.title('Euclidean Distance in PC3',y=1.03)
        plt.suptitle('Figure {}. {}'.format(self.fignumber, title),fontsize=self.titleSizeSup); self.pdf.savefig(fig);self.fignumber+=1

    def createBasisGraph(self, data):
        plt.figure(self.nFig)
        plt.suptitle('Basis')
        nBase = data.shape[1]
        # The cols number of subplot is 
        nSubCols = nBase / 10
        if nSubCols > 0:
            if nBase % 2 == 0:
                nSubRows = nBase / nSubCols
            else:
                nSubRows = nBase / nSubCols + 1
        else:
            nSubRows = nBase 
            nSubCols = 1
        # freqList = np.fft.fftfreq(513, d = 1.0 / 44100)

        for i in range(nBase):
            nowFig = self.nFig + (i / nSubRows) + 1
            # Because Index of graph is start by 1, The Graph index start from i + 1.
            plt.subplot(nSubRows, nSubCols, i + 1)
            plt.tick_params(labelleft='off', labelbottom='off')

            # FIXME
            #plt.ylabel(self.st5[i%12] + str(i/12  + 1))
            plt.ylabel(str(i))
            plt.plot(data[:,i])
        # Beacuse I want to add lable in bottom, xlabel is declaration after loop.
        plt.tick_params(labelleft='off', labelbottom='on')
        plt.xlabel('frequency [Hz]')

        #self.nFig += nowFig
        self.nFig += 1