def plot(coefs_files, digit=0, mixture=0, axis=0):
    import numpy as np
    import matplotlib.pylab as plt
    import amitgroup as ag

    for coefs_file in coefs_files:
        coefs_data = np.load(coefs_file)
        var = coefs_data['prior_var']
        samples = coefs_data['samples']
        llh_var = coefs_data['llh_var']

        var_flat = ag.util.wavelet.smart_flatten(var[digit,mixture,axis])
        last_i = len(var_flat)-1
        plt.xlim((0, last_i))

        #imdef = ag.util.DisplacementFieldWavelet((32, 32), 'db4', penalty=100

        if len(coefs_files) == 1:
            add = ""
        else:   
            add = " ({0})".format(coefs_file.name)

        plt.subplot(121)
        plt.semilogy(1/var_flat, label="ML"+add)
        plt.legend(loc=0)
        plt.xlabel('Coefficient')
        plt.ylabel('Precision $\lambda$')
        plt.xlim((0, 63))

        plt.subplot(122)
        plt.imshow(1/llh_var[digit,mixture], interpolation='nearest')
        plt.xlabel("Likelihood precision $\lambda'$")
        plt.colorbar()
    plt.show()

def plot_grid_experiment_results(grid_results, params, metrics):
    global plt
    params = sorted(params)
    grid_params = grid_results.grid_params
    plt.figure(figsize=(8, 6))
    for metric in metrics:
        grid_params_shape = [len(grid_params[k]) for k in sorted(grid_params.keys())]
        params_max_out = [(1 if k in params else 0) for k in sorted(grid_params.keys())]
        results = np.array([e.results.get(metric, 0) for e in grid_results.experiments])
        results = results.reshape(*grid_params_shape)
        for axis, included_in_params in enumerate(params_max_out):
            if not included_in_params:
                results = np.apply_along_axis(np.max, axis, results)

        print results
        params_shape = [len(grid_params[k]) for k in sorted(params)]
        results = results.reshape(*params_shape)

        if len(results.shape) == 1:
            results = results.reshape(-1,1)
        import matplotlib.pylab as plt

        #f.subplots_adjust(left=.2, right=0.95, bottom=0.15, top=0.95)
        plt.imshow(results, interpolation='nearest', cmap=plt.cm.hot)
        plt.title(str(grid_results.name) + " " + metric)

        if len(params) == 2:
            plt.xticks(np.arange(len(grid_params[params[1]])), grid_params[params[1]], rotation=45)
        plt.yticks(np.arange(len(grid_params[params[0]])), grid_params[params[0]])
        plt.colorbar()
        plt.show()

 def plot_bernoulli_matrix(self, show_npfs=False):
   """
   Plot the heatmap of the Bernoulli matrix 
   @self
   @show_npfs - Highlight NPFS detections [Boolean] 
   """
   matrix = self.Bernoulli_matrix
   if show_npfs == False:
     plot = plt.imshow(matrix)
     plot.set_cmap('hot')
     plt.colorbar()
     plt.xlabel("Bootstraps")
     plt.ylabel("Feature")
     plt.show()
   else:
     for i in self.selected_features:
       for k in range(len(matrix[i])):
         matrix[i,k] = .5
     plot = plt.imshow(matrix)
     plot.set_cmap('hot')
     plt.xlabel("Bootstraps")
     plt.ylabel("Feature")
     plt.colorbar()
     plt.show()
   return None

    def plotscalefactors(self):
        '''
        Plots the distribution of scalefactors between consecutive images.
        This is to check manually that the hdr image assembly is done correctly.
        '''
        import matplotlib.pylab as plt

        print('Scalefactors for HDR-assembling are', self.scalefactors)
        print('Standard deviations for Scalefactors are', self.scalefactorsstd)

        # Plots the scalefactordistribution and scalefactors for each pixelvalue
        self.scaleforpix = range(len(self.raw))
        for n in range(len(self.raw) - 1):
            A = self.getimgquotient(n)
            B = self._getrealimg(n + 1)
            fig = plt.figure()
            plt.xlabel('x [pixel]')
            plt.ylabel('y [pixel]')
            fig.set_size_inches(10, 10)
            plt.imshow(A)
            plt.clim([0.95 * A[np.isfinite(A)].min(), 1.05 * A[np.isfinite(A)].max()])
            plt.colorbar()
            fig = plt.figure()
            linplotdata = np.array([B[np.isfinite(B)].flatten(), A[np.isfinite(B)].flatten()])
            plt.plot(linplotdata[0, :], linplotdata[1, :], 'ro')

    def XXtest5_regrid(self):
        srcF = cdms2.open(sys.prefix + \
                              '/sample_data/so_Omon_ACCESS1-0_historical_r1i1p1_185001-185412_2timesteps.nc')
        so = srcF('so')[0, 0, ...]
        clt = cdms2.open(sys.prefix + '/sample_data/clt.nc')('clt')
        dstData = so.regrid(clt.getGrid(), 
                            regridTool = 'esmf', 
                            regridMethod='conserve')

        if self.pe == 0:
            dstDataMask = (dstData == so.missing_value)
            dstDataFltd = dstData * (1 - dstDataMask)
            zeroValCnt = (dstData == 0).sum()
            if so.missing_value > 0:
                dstDataMin = dstData.min()
                dstDataMax = dstDataFltd.max()
            else:
                dstDataMin = dstDataFltd.min()
                dstDataMax = dstData.max()
                zeroValCnt = (dstData == 0).sum()
            print 'Number of zero valued cells', zeroValCnt
            print 'min/max value of dstData: %f %f' % (dstDataMin, dstDataMax)                   
            self.assertLess(dstDataMax, so.max())
            if False:
                pylab.figure(1)
                pylab.pcolor(so, vmin=20, vmax=40)
                pylab.colorbar()
                pylab.title('so')
                pylab.figure(2)
                pylab.pcolor(dstData, vmin=20, vmax=40)
                pylab.colorbar()
                pylab.title('dstData')

def reconstructContactMap(map, datavec):
    """ Plots a given vector as a contact map

    Parameters
    ----------
    map : np.ndarray 2D
        The map from a MetricData object
    datavec : np.ndarray
        The data we want to plot in a 2D map
    """
    map = np.array(map, dtype=int)
    atomidx = np.unique(map.flatten()).astype(int)
    mask = np.zeros(max(atomidx)+1, dtype=int)
    mask[atomidx] = range(len(atomidx))

    # Create a new map which maps from vector indexes to matrix indexes
    newmap = np.zeros(np.shape(map), dtype=int)
    newmap[:, 0] = mask[map[:, 0]]
    newmap[:, 1] = mask[map[:, 1]]

    contactmap = np.zeros((len(atomidx), len(atomidx)))
    for i in range(len(datavec)):
        contactmap[newmap[i, 0], newmap[i, 1]] = datavec[i]
        contactmap[newmap[i, 1], newmap[i, 0]] = datavec[i]

    from matplotlib import pylab as plt
    plt.imshow(contactmap, interpolation='nearest', aspect='equal')
    plt.colorbar()
    #plt.axis('off')
    #plt.tick_params(axis='x', which='both', bottom='off', top='off', labelbottom='off')
    #plt.tick_params(axis='y', which='both', left='off', right='off', labelleft='off')
    plt.show()

 def test2_3x4_to_5x7_cart(self):
     # Test non-periodic grid returning double grid resolution
     roESMP = CdmsRegrid(self.fromGrid3x4, self.toGrid5x7,
                         dtype = self.data3x4.dtype,
                         srcGridMask = self.data3x4.mask,
                         regridTool = 'ESMP',
                         periodicity = 0,
                         regridMethod = 'Conserve', 
                         coordSys = 'cart')
     diag = {'srcAreas':0, 'dstAreas':0, 
             'srcAreaFractions':0, 'dstAreaFractions':0}
     ESMP5x7 = roESMP(self.data3x4, diag = diag)
     dstMass = (ESMP5x7 * diag['dstAreas']).sum()
     srcMass = (self.data3x4 * diag['srcAreas'] \
                             * diag['srcAreaFractions']).sum()
     if False:
         pylab.figure(1)
         pylab.pcolor(self.data3x4)
         pylab.colorbar()
         pylab.title('original self.data3x4')
         pylab.figure(2)
         pylab.pcolor(ESMP5x7)
         pylab.colorbar()
         pylab.title('interpolated ESMP5x7')
     self.assertLess(abs(srcMass - dstMass), self.eps)
     self.assertEqual(self.data3x4[0,0], ESMP5x7[0,0])
     self.assertEqual(1.0, ESMP5x7[0,0])
     self.assertEqual(0.25, ESMP5x7[1,1])
     self.assertEqual(0.0, ESMP5x7[2,2])

def plot_block( files, m, n ):
    """
        - files is directory for cell
        - m starting index, n ending index
        Use pylab.imshow() to plot a frame (npy array). 
        Note: The boundary should be removed, thus there will not be a halo
        """
    #Get all frames
    if os.path.isdir (files):
        fdir = files + '/'
        dlist = os.listdir(fdir)
        frames = []
        for f in dlist:
            if f.endswith('npy') and not os.path.isdir(fdir+f):
                frames.append(f)
    else:
        print 'Error - input is not a directory'
        return
    #Sort frames
    frames . sort(key=natural_key)
    stack = []
    #load npy arrays
    for ind in xrange(len(frames)):
        if ind >= m and ind <= n:
            stack.append(numpy.load(files + frames[ind]))
    #Create 3d array
    d = numpy.dstack(stack)
    d = numpy.rollaxis(d,-1)
    fig = plt.figure()
    plt.title("RBC Stack")
    plt.imshow(d[1])
    plt.colorbar()
    plt.show()

def plot_prof(obj):
  '''Function that plots a vertical profile of scattering from
  TOC to BOC.
  Input: rt_layer instance object.
  Output: plot of scattering.
  '''
  I = obj.Inodes[:,1,:] \
      / -obj.mu_s
  y = np.linspace(obj.Lc, 0., obj.K+1)
  x = obj.views*180./np.pi
  xm = np.array([])
  for i in np.arange(0,len(x)-1):
    nx = x[i] + (x[i+1]-x[i])/2.
    xm = np.append(xm,nx)
  xm = np.insert(xm,0,0.)
  xm = np.append(xm,180.)
  xx, yy = np.meshgrid(xm, y)
  plt.pcolormesh(xx,yy,I)
  plt.colorbar()
  plt.title('Canopy Fluxes')
  plt.xlabel('Exit Zenith Angle')
  plt.ylabel('Cumulative LAI (0=soil)')
  plt.arrow(135.,3.5,0.,-3.,head_width=5.,head_length=.2,\
      fc='k',ec='k')
  plt.text(140.,2.5,'Downwelling Flux',rotation=90)
  plt.arrow(45.,.5,0.,3.,head_width=5.,head_length=.2,\
      fc='k',ec='k')
  plt.text(35.,2.5,'Upwelling Flux',rotation=270)
  plt.show()

def plotWeightChanges():
    if f.usestdp:
        # create plot
        figh = figure(figsize=(1.2*8,1.2*6))
        figh.subplots_adjust(left=0.02) # Less space on left
        figh.subplots_adjust(right=0.98) # Less space on right
        figh.subplots_adjust(top=0.96) # Less space on bottom
        figh.subplots_adjust(bottom=0.02) # Less space on bottom
        figh.subplots_adjust(wspace=0) # More space between
        figh.subplots_adjust(hspace=0) # More space between
        h = axes()

        # create data matrix
        wcs = [x[-1][-1] for x in f.allweightchanges] # absolute final weight
        wcs = [x[-1][-1]-x[0][-1] for x in f.allweightchanges] # absolute weight change
        pre,post,recep = zip(*[(x[0],x[1],x[2]) for x in f.allstdpconndata])
        ncells = int(max(max(pre),max(post))+1)
        wcmat = zeros([ncells, ncells])

        for iwc,ipre,ipost,irecep in zip(wcs,pre,post,recep):
            wcmat[int(ipre),int(ipost)] = iwc *(-1 if irecep>=2 else 1)

        # plot
        imshow(wcmat,interpolation='nearest',cmap=bicolormap(gap=0,mingreen=0.2,redbluemix=0.1,epsilon=0.01))
        xlabel('post-synaptic cell id')
        ylabel('pre-synaptic cell id')
        h.set_xticks(f.popGidStart)
        h.set_yticks(f.popGidStart)
        h.set_xticklabels(f.popnames)
        h.set_yticklabels(f.popnames)
        h.xaxif.set_ticks_position('top')
        xlim(-0.5,ncells-0.5)
        ylim(ncells-0.5,-0.5)
        clim(-abs(wcmat).max(),abs(wcmat).max())
        colorbar()

    def plot(self, normalized=False, backend=None, ax=None, **kwargs):
        """
        Plots confusion matrix
        """
        df = self.to_dataframe(normalized)

        try:
            cmap = kwargs['cmap']
        except:
            cmap = plt.cm.gray_r

        title = self.title

        if normalized:
            title += " (normalized)"

        if backend is None:
            backend = self.backend

        if backend == Backend.Matplotlib:
            #if ax is None:
            fig, ax = plt.subplots(figsize=(9, 8))
            plt.imshow(df, cmap=cmap, interpolation='nearest') # imshow / matshow
            ax.set_title(title)

            tick_marks_col = np.arange(len(df.columns))
            tick_marks_idx = tick_marks_col.copy()

            ax.set_yticks(tick_marks_idx)
            ax.set_xticks(tick_marks_col)
            ax.set_xticklabels(df.columns, rotation=45, ha='right')
            ax.set_yticklabels(df.index)
            
            ax.set_ylabel(df.index.name)
            ax.set_xlabel(df.columns.name)

            (N_min, N_max) = (0, self.max())
            if N_max > COLORBAR_TRIG:
                plt.colorbar() # Continuous colorbar
            else:
                # Discrete colorbar
                ax2 = fig.add_axes([0.93, 0.1, 0.03, 0.8])
                bounds = np.arange(N_min, N_max + 2, 1)
                norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
                cb = mpl.colorbar.ColorbarBase(ax2, cmap=cmap, norm=norm, spacing='proportional', ticks=bounds, boundaries=bounds, format='%1i')

            return(ax)


        elif backend == Backend.Seaborn:
            import seaborn as sns
            ax = sns.heatmap(df, **kwargs)
            return(ax)
            # You should test this yourself
            # because I'm facing an issue with Seaborn under Mac OS X (2015-04-26)
            # RuntimeError: Cannot get window extent w/o renderer
            #sns.plt.show()

        else:
            raise(NotImplementedError("backend=%r not allowed" % backend))

def plotConn():
    # Create plot
    figh = figure(figsize=(8,6))
    figh.subplots_adjust(left=0.02) # Less space on left
    figh.subplots_adjust(right=0.98) # Less space on right
    figh.subplots_adjust(top=0.96) # Less space on bottom
    figh.subplots_adjust(bottom=0.02) # Less space on bottom
    figh.subplots_adjust(wspace=0) # More space between
    figh.subplots_adjust(hspace=0) # More space between
    h = axes()
    totalconns = zeros(shape(f.connprobs))
    for c1 in range(size(f.connprobs,0)):
        for c2 in range(size(f.connprobs,1)):
            for w in range(f.nreceptors):
                totalconns[c1,c2] += f.connprobs[c1,c2]*f.connweights[c1,c2,w]*(-1 if w>=2 else 1)
    imshow(totalconns,interpolation='nearest',cmap=bicolormap(gap=0))

    # Plot grid lines
    hold(True)
    for pop in range(f.npops):
        plot(array([0,f.npops])-0.5,array([pop,pop])-0.5,'-',c=(0.7,0.7,0.7))
        plot(array([pop,pop])-0.5,array([0,f.npops])-0.5,'-',c=(0.7,0.7,0.7))

    # Make pretty
    h.set_xticks(range(f.npops))
    h.set_yticks(range(f.npops))
    h.set_xticklabels(f.popnames)
    h.set_yticklabels(f.popnames)
    h.xaxis.set_ticks_position('top')
    xlim(-0.5,f.npops-0.5)
    ylim(f.npops-0.5,-0.5)
    clim(-abs(totalconns).max(),abs(totalconns).max())
    colorbar()

def plot(rho, u, uLB, tau, rho_history, zdjecia, image, nx, maxIter ):
#    plt.figure(figsize=(15,15))
#    plt.subplot(4, 1, 1)
#    plt.imshow(u[1,:,0:50],vmin=-uLB*.15, vmax=uLB*.15, interpolation='none')#,cmap=cm.seismic
#    plt.colorbar()
    plt.rcParams["figure.figsize"] = (15,15)
    plt.subplot(5, 1, 1)
    plt.imshow(sqrt(u[0]**2+u[1]**2),vmin=0, vmax=uLB*1.6)#,cmap=cm.seismic
    plt.colorbar()
    plt.title('tau = {:f}'.format(tau))      
    
    plt.subplot(5, 1, 2)
    plt.imshow(u[0,:,:30],  interpolation='none')#,cmap=cm.seismicvmax=uLB*1.6,
    plt.colorbar()
    plt.title('tau = {:f}'.format(tau))  
    
    plt.subplot(5, 1, 3)
    plt.imshow(rho, interpolation='none' )#,cmap=cm.seismic
    plt.title('rho')   
    
    plt.subplot(5, 1,4)
    plt.title(' history rho')
    plt.plot(linspace(0,len(rho_history),len(rho_history)),rho_history)
    plt.xlim([0,maxIter])   
    
    plt.subplot(5, 1,5)
    plt.title(' u0 middle develop')
    plt.plot(linspace(0,nx,len(u[0,20,:])), u[1,20,:])
    plt.tight_layout()                  
    
    plt.savefig(path.join(zdjecia,'f{0:06d}.png'.format(image)))
    plt.clf();
        

def show_slice(data, slice_no=0, clim=None):
    """ Visualize the reconstructed slice.

    Parameters
    -----------
    data : ndarray
        3-D matrix of stacked reconstructed slices.

    slice_no : scalar, optional
        The index of the slice to be imaged.
    """
    plt.figure(figsize=(7, 7))
    if len(data.shape) is 2:
        plt.imshow(data,
                   interpolation='none',
                   cmap='gray')
        plt.colorbar()
    elif len(data.shape) is 3:
        plt.imshow(data[slice_no, :, :],
                   interpolation='none',
                   cmap='gray')
        plt.colorbar()
    if clim is not None:
        plt.clim(clim)
    plt.show()

	def drown_field(self,data,mask,drown):
		''' drown_field is a wrapper around the fortran code fill_msg_grid.
		depending on the output geometry, applies land extrapolation on 1 or N levels'''
		if self.geometry == 'surface':
			for kz in _np.arange(self.nz):
				tmpin = data[kz,:,:].transpose()
				if self.debug and kz == 0:
					tmpin_plt = _np.ma.masked_values(tmpin,self.xmsg)
					_plt.figure() ; _plt.contourf(tmpin_plt.transpose(),40) ; _plt.colorbar() ; 
					_plt.title('normalized before drown')
				if drown == 'ncl':
					tmpout = _fill.mod_poisson.poisxy1(tmpin,self.xmsg, self.guess, self.gtype, \
					self.nscan, self.epsx, self.relc)
				elif drown == 'sosie':
					tmpout = _mod_drown_sosie.mod_drown.drown(self.kew,tmpin,mask[kz,:,:].T,\
					nb_inc=200,nb_smooth=40)
				data[kz,:,:] = tmpout.transpose()
				if self.debug and kz == 0:
					_plt.figure() ; _plt.contourf(tmpout.transpose(),40) ; _plt.colorbar() ; 
					_plt.title('normalized after drown') 
					_plt.show()
		elif self.geometry == 'line':
			tmpin = data[:,:].transpose()
			if drown == 'ncl':
				tmpout = _fill.mod_poisson.poisxy1(tmpin,self.xmsg, self.guess, self.gtype, \
				self.nscan, self.epsx, self.relc)
			elif drown == 'sosie':
				tmpout = _mod_drown_sosie.mod_drown.drown(self.kew,tmpin,mask[:,:].T,\
				nb_inc=200,nb_smooth=40)
			data[:,:] = tmpout.transpose()
		return data

def plot_image(*positional_parameters,title="TITLE",xtitle=r"X",ytitle=r"Y",cmap=None,aspect=None,show=1):

    n_arguments = len(positional_parameters)
    if n_arguments == 1:
        z = positional_parameters[0]
        x = np.arange(0,z.shape[0])
        y = np.arange(0,z.shape[0])
    elif n_arguments == 2:
        z = positional_parameters[0]
        x = positional_parameters[1]
        y = positional_parameters[1]
    elif n_arguments == 3:
        z = positional_parameters[0]
        x = positional_parameters[1]
        y = positional_parameters[2]
    else:
        raise Exception("Bad number of inputs")


    fig = plt.figure()

    # cmap = plt.cm.Greys
    plt.imshow(z.T,origin='lower',extent=[x[0],x[-1],y[0],y[-1]],cmap=cmap,aspect=aspect)
    plt.colorbar()
    ax = fig.gca()
    ax.set_xlabel(xtitle)
    ax.set_ylabel(ytitle)

    plt.title(title)

    if show:
        plt.show()

    return fig

	def perform_interpolation(self,dataextrap,regridme,field_src,field_target,use_locstream):
		data = self.allocate()
		if self.geometry == 'surface':
			for kz in _np.arange(self.nz):
				field_src.data[:] = dataextrap[kz,:,:].transpose()
				field_target = regridme(field_src, field_target)
				if use_locstream:
					if self.nx == 1:
						data[kz,:,0] = field_target.data.copy()
					elif self.ny == 1:
						data[kz,0,:] = field_target.data.copy()
				else:
					data[kz,:,:] = field_target.data.transpose()[self.jmin:self.jmax+1, \
					                                             self.imin:self.imax+1]
					if self.debug and kz == 0:
						data_target_plt = _np.ma.masked_values(data[kz,:,:],self.xmsg)
						#data_target_plt = _np.ma.masked_values(field_target.data,self.xmsg)
						_plt.figure() ; _plt.contourf(data_target_plt[:,:],40) ; _plt.colorbar() ; 
						_plt.title('regridded') ; _plt.show()
		elif self.geometry == 'line':
			field_src.data[:] = dataextrap[:,:].transpose()
			field_target = regridme(field_src, field_target)
			if use_locstream:
				data[:,:] = _np.reshape(field_target.data.transpose(),(self.ny,self.nx))
			else:
				data[:,:] = field_target.data.transpose()[self.jmin:self.jmax+1,self.imin:self.imax+1]
		return data

def generate_matrix_visualization(known_words, bookmark_words):

    m = []
    for i in range(0,100):
        m.append([])
        for j in range (0, 100):
            if (i*100+j) in bookmark_words:
                m[i].append(0.65)
            elif (i*100+j) in known_words:
                m[i].append(0.4)
            else:
                m[i].append(0.2)
    m.reverse()
	
    # we need this next line to scale the color scheme
    m[0][0]=0
    matrix = numpy.matrix(m)

	fig = plt.figure()
	ax = fig.add_subplot(1,1,1)
	ax.set_aspect('equal')
	plt.imshow(matrix, interpolation='none')
	plt.set_cmap('hot')
	plt.colorbar()
	plt.show()

 def plot(x,y,field,filename,c=200):
     plt.figure()
     # define grid.
     xi = np.linspace(min(x),max(x),100)
     yi = np.linspace(min(y),max(y),100)
     # grid the data.
     si_lin = griddata((x, y), field, (xi[None,:], yi[:,None]), method='linear')
     si_cub = griddata((x, y), field, (xi[None,:], yi[:,None]), method='linear')
     print np.min(field)
     print np.max(field)
     plt.subplot(211)
     # contour the gridded data, plotting dots at the randomly spaced data points.
     CS = plt.contour(xi,yi,si_lin,c,linewidths=0.5,colors='k')
     CS = plt.contourf(xi,yi,si_lin,c,cmap=plt.cm.jet)
     plt.colorbar() # draw colorbar
     # plot data points.
     #    plt.scatter(x,y,marker='o',c='b',s=5)
     plt.xlim(min(x),max(x))
     plt.ylim(min(y),max(y))
     plt.title('Lineaarinen interpolointi')
     #plt.tight_layout()
     plt.subplot(212)
     # contour the gridded data, plotting dots at the randomly spaced data points.
     CS = plt.contour(xi,yi,si_cub,c,linewidths=0.5,colors='k')
     CS = plt.contourf(xi,yi,si_cub,c,cmap=plt.cm.jet)
     plt.colorbar() # draw colorbar
     # plot data points.
     #    plt.scatter(x,y,marker='o',c='b',s=5)
     plt.xlim(min(x),max(x))
     plt.ylim(min(y),max(y))
     plt.title('Kuubinen interpolointi')
     plt.savefig(filename)

    def showladderpca(self):

        import mdp
        from matplotlib import pylab as plt
        import pprint
        import math

        cerr('calculating PCA & plotting')
        peak_sizes = sorted(list([ x.rtime for x in self.alleles ]))
        #peak_sizes = sorted( peak_sizes )[:-5]
        #pprint.pprint(peak_sizes)
        #comps = algo.simple_pca( peak_sizes )
        #algo.plot_pca(comps, peak_sizes)

        from fatools.lib import const
        std_sizes = const.ladders['LIZ600']['sizes']

        x = std_sizes
        y = [ x * 0.1 for x in peak_sizes ]

        D = np.zeros( (len(y), len(x)) )
        for i in range(len(y)):
            for j in range(len(x)):
                D[i,j] = math.exp( ((x[j] - y[i]) * 0.001) ** 2 )

        pprint.pprint(D)
        im = plt.imshow(D, interpolation='nearest', cmap='Reds')
        plt.gca().invert_yaxis()
        plt.xlabel("STD")
        plt.ylabel("PEAK")
        plt.grid()
        plt.colorbar()
        plt.show()