def plot_frontier(self,frontier_only=False,plot_samples=True) :
		""" Plot the frontier"""
		frontier = self.frontier
		frontier_energy = self.frontier_energy
		feat1,feat2 = self.feats	
	
		pl.figure()
		if not frontier_only :	
			ll_list1,ll_list2 = zip(*self.all_seq_energy)
			pl.plot(ll_list1,ll_list2,'b*')
		if plot_samples :
			ll_list1,ll_list2 = zip(*self.sample_seq_energy)
			pl.plot(ll_list1,ll_list2,'g*')
					
		pl.plot(*zip(*sorted(frontier_energy)),color='magenta',\
			marker='*',	linestyle='dashed')
		ctr = dict(zip(set(frontier_energy),[0]*
			len(set(frontier_energy))))
		for i,e in enumerate(frontier_energy) : 
			ctr[e] += 1
			pl.text(e[0],e[1]+0.1*ctr[e],str(i),fontsize=10)
			pl.text(e[0]+0.4,e[1]+0.1*ctr[e],frontier[i],fontsize=9)	
		pl.xlabel('Energy:'+feat1)
		pl.ylabel('Energy:'+feat2)
		pl.title('Energy Plot')
		xmin,xmax = pl.xlim()
		ymin,ymax = pl.ylim()
		pl.xlim(xmin,xmax)
		pl.ylim(ymin,ymax)
		pic_dir = '../docs/tex/pics/'
		pl.savefig(pic_dir+self.name+'.pdf')
		pl.savefig(pic_dir+self.name+'.png')

def trace(data, name, format='png', datarange=(None, None), suffix='', path='./', rows=1, columns=1, 
    num=1, last=True, fontmap = None, verbose=1):
    """
    Generates trace plot from an array of data.

    :Arguments:
        data: array or list
            Usually a trace from an MCMC sample.

        name: string
            The name of the trace.
            
        datarange: tuple or list
            Preferred y-range of trace (defaults to (None,None)).

        format (optional): string
            Graphic output format (defaults to png).

        suffix (optional): string
            Filename suffix.

        path (optional): string
            Specifies location for saving plots (defaults to local directory).
            
        fontmap (optional): dict
            Font map for plot.

    """

    if fontmap is None: fontmap = {1:10, 2:8, 3:6, 4:5, 5:4}

    # Stand-alone plot or subplot?
    standalone = rows==1 and columns==1 and num==1

    if standalone:
        if verbose>0:
            print_('Plotting', name)
        figure()

    subplot(rows, columns, num)
    pyplot(data.tolist())
    ylim(datarange)

    # Plot options
    title('\n\n   %s trace'%name, x=0., y=1., ha='left', va='top', fontsize='small')

    # Smaller tick labels
    tlabels = gca().get_xticklabels()
    setp(tlabels, 'fontsize', fontmap[rows/2])

    tlabels = gca().get_yticklabels()
    setp(tlabels, 'fontsize', fontmap[rows/2])

    if standalone:
        if not os.path.exists(path):
            os.mkdir(path)
        if not path.endswith('/'):
            path += '/'
        # Save to file
        savefig("%s%s%s.%s" % (path, name, suffix, format))

def param_set_averages_plot(results):
    averages_ocr = [
        a[1] for a in sorted(
            param_set_averages(results, metric='ocr').items(),
            key=lambda x: int(x[0].split('-')[1]))
    ]
    averages_q = [
        a[1] for a in sorted(
            param_set_averages(results, metric='q').items(),
            key=lambda x: int(x[0].split('-')[1]))
    ]
    averages_mse = [
        a[1] for a in sorted(
            param_set_averages(results, metric='mse').items(),
            key=lambda x: int(x[0].split('-')[1]))
    ]
    fig = plt.figure(figsize=(6, 4))
    # plt.tight_layout()
    plt.plot(averages_ocr, label='OCR', linewidth=2.0)
    plt.plot(averages_q, label='Q', linewidth=2.0)
    plt.plot(averages_mse, label='MSE', linewidth=2.0)
    plt.ylim([0, 1])
    plt.xlabel(u'Paslėptų neuronų skaičius')
    plt.ylabel(u'Vidurinė Q įverčio pokyčio reikšmė')
    plt.grid(True)
    plt.tight_layout()
    plt.legend(loc='lower right')
    plt.show()

def plotB3reg():
    w=loadStanFit('revE2B3BHreg.fit')
    printCI(w,'mmu')
    printCI(w,'mr')
    for b in range(2):
        subplot(1,2,b+1)
        plt.title('')
        px=np.array(np.linspace(-0.5,0.5,101),ndmin=2)
        a0=np.array(w['mmu'][:,b],ndmin=2).T
        a1=np.array(w['mr'][:,b],ndmin=2).T
        y=np.concatenate([sap(a0+a1*px,97.5,axis=0),sap(a0+a1*px[:,::-1],2.5,axis=0)])
        x=np.squeeze(np.concatenate([px,px[:,::-1]],axis=1))
        plt.plot(px[0,:],np.median(a0)+np.median(a1)*px[0,:],'red')
        #plt.plot([-1,1],[0.5,0.5],'grey')
        ax=plt.gca()
        ax.set_aspect(1)
        ax.add_patch(plt.Polygon(np.array([x,y]).T,alpha=0.2,fill=True,fc='red',ec='w'))
        y=np.concatenate([sap(a0+a1*px,75,axis=0),sap(a0+a1*px[:,::-1],25,axis=0)])
        ax.add_patch(plt.Polygon(np.array([x,y]).T,alpha=0.2,fill=True,fc='red',ec='w'))
        man=np.array([-0.4,-0.2,0,0.2,0.4])
        mus=[]
        for m in range(len(man)):
            mus.append(loadStanFit('revE2B3BH%d.fit'%m)['mmu'][:,b])
        mus=np.array(mus).T
        errorbar(mus,x=man)
        ax.set_xticks(man)
        plt.xlim([-0.5,0.5])
        plt.ylim([-0.4,0.8])
        #plt.xlabel('Manipulated Displacement')
        if b==0:
            plt.ylabel('Perceived Displacemet')
            plt.gca().set_yticklabels([])
        subplot_annotate()
    plt.text(-1.1,-0.6,'Pivot Displacement',fontsize=8);

def plot_sphere_x( s, fname ):
  """ put plot of ionization fractions from sphere `s` into fname """

  plt.figure()
  s.Edges.units = 'kpc'
  s.r_c.units = 'kpc'
  xx = s.r_c
  L = s.Edges[-1]

  plt.plot( xx, np.log10( s.xHe1 ),
            color='green', ls='-', label = r'$x_{\rm HeI}$' )
  plt.plot( xx, np.log10( s.xHe2 ),
            color='green', ls='--', label = r'$x_{\rm HeII}$' )
  plt.plot( xx, np.log10( s.xHe3 ),
            color='green', ls=':', label = r'$x_{\rm HeIII}$' )

  plt.plot( xx, np.log10( s.xH1 ),
            color='red', ls='-', label = r'$x_{\rm HI}$' )
  plt.plot( xx, np.log10( s.xH2 ),
            color='red', ls='--', label = r'$x_{\rm HII}$' )

  plt.xlim( -L/20, L+L/20 )
  plt.xlabel( 'r_c [kpc]' )

  plt.ylim( -4.5, 0.2 )
  plt.ylabel( 'log 10 ( x )' )

  plt.grid()
  plt.legend(loc='best', ncol=2)
  plt.tight_layout()
  plt.savefig( 'doc/img/x_' + fname )

    def InitializePlot(self, goal_config):
        self.fig = pl.figure()
        lower_limits, upper_limits = self.boundary_limits
        pl.xlim([lower_limits[0], upper_limits[0]])
        pl.ylim([lower_limits[1], upper_limits[1]])
        pl.plot(goal_config[0], goal_config[1], 'gx')

        # Show all obstacles in environment
        for b in self.robot.GetEnv().GetBodies():
            if b.GetName() == self.robot.GetName():
                continue
            bb = b.ComputeAABB()
            pl.plot([bb.pos()[0] - bb.extents()[0],
                     bb.pos()[0] + bb.extents()[0],
                     bb.pos()[0] + bb.extents()[0],
                     bb.pos()[0] - bb.extents()[0],
                     bb.pos()[0] - bb.extents()[0]],
                    [bb.pos()[1] - bb.extents()[1],
                     bb.pos()[1] - bb.extents()[1],
                     bb.pos()[1] + bb.extents()[1],
                     bb.pos()[1] + bb.extents()[1],
                     bb.pos()[1] - bb.extents()[1]], 'r')
                    
                     
        pl.ion()
        pl.show()

    def rmsdSpreadSubplot(multiplier=1.0, layout=(-1, -1)):
        rmsd_data   = dict( (e, rad_data[e]['innov'][quant])  for e in rad_data.iterkeys() )
        spread_data = dict( (e, rad_data[e]['spread'][quant]) for e in rad_data.iterkeys() )

        times = temp.getTimes()
        n_t = len(times)

        for exp, exp_name in exp_names.iteritems():
            pylab.plot(sawtooth(times, times)[:(n_t + 1)], rmsd_data[exp][:(n_t + 1)], color=colors[exp], linestyle='-')
            pylab.plot(times[(n_t / 2):], rmsd_data[exp][n_t::2], color=colors[exp], linestyle='-')
 
        for exp, exp_name in exp_names.iteritems():
            pylab.plot(sawtooth(times, times)[:(n_t + 1)], spread_data[exp][:(n_t + 1)], color=colors[exp], linestyle='--')
            pylab.plot(times[(n_t / 2):], spread_data[exp][n_t::2], color=colors[exp], linestyle='--')

        ylim = pylab.ylim()
        pylab.plot(times, -1 * np.ones((len(times),)), color='#999999', linestyle='-', label="RMS Innovation")
        pylab.plot(times, -1 * np.ones((len(times),)), color='#999999', linestyle='--', label="Spread")

        pylab.axhline(y=7, color='k', linestyle=':')
        pylab.axvline(x=14400, color='k', linestyle=':')

        pylab.ylabel("RMS Innovation/Spread (dBZ)", size='large')

        pylab.xlim(times[0], times[-1])
        pylab.ylim(ylim)

        pylab.legend(loc=4)

        pylab.xticks(times[::2], [ "" for t in times[::2] ])
        pylab.yticks(size='x-large')
        return

    def InitializePlot(self, goal_config):  # default
        self.fig = pl.figure()
        pl.xlim([self.lower_limits[0], self.upper_limits[0]])
        pl.ylim([self.lower_limits[1], self.upper_limits[1]])
        pl.plot(goal_config[0], goal_config[1], "gx")

        # Show all obstacles in environment
        for b in self.robot.GetEnv().GetBodies():
            if b.GetName() == self.robot.GetName():
                continue
            bb = b.ComputeAABB()
            pl.plot(
                [
                    bb.pos()[0] - bb.extents()[0],
                    bb.pos()[0] + bb.extents()[0],
                    bb.pos()[0] + bb.extents()[0],
                    bb.pos()[0] - bb.extents()[0],
                    bb.pos()[0] - bb.extents()[0],
                ],
                [
                    bb.pos()[1] - bb.extents()[1],
                    bb.pos()[1] - bb.extents()[1],
                    bb.pos()[1] + bb.extents()[1],
                    bb.pos()[1] + bb.extents()[1],
                    bb.pos()[1] - bb.extents()[1],
                ],
                "r",
            )

        pl.ion()
        pl.show()

def RosenbrockTest():
    nDim = 3
    numOfParticles = 20
    maxIteration = 200
    minX = array([-5.0]*nDim)
    maxX = array([5.0]*nDim)
    maxV = 0.2*(maxX - minX)
    minV = -1.0*maxV
    numOfTrial = 20
    gBest = array([0.0]*maxIteration)
    for i in xrange(numOfTrial):
        p1 = RPSO.PSOProblem(nDim, numOfParticles, maxIteration, minX, maxX, minV, maxV, RPSO.Rosenbrock)
        p1.run()
        gBest = gBest + p1.gBestArray[:maxIteration]
    gBest = gBest / numOfTrial
    pylab.title('$G_{best}$ over 20 trials')
    pylab.xlabel('The $N^{th}$ Iteratioin')
    pylab.ylabel('Average gBest over '+str(numOfTrial)+' runs (logscale)')
    pylab.grid(True)
#    pylab.yscale('log')
    ymin, ymax = -1.5, 2.5
    ystep = 0.5
    pylab.ylim(ymin, ymax)
    yticks = linspace(ymin, ymax, (ymax-ymin)/ystep+1)
    pylab.yticks(tuple(yticks),tuple(map(str,yticks)))
    pylab.plot(range(maxIteration), log10(gBest),'-', label='Global best')
    pylab.legend()
    pylab.show()

def tracks_movie(base, skip=1, frames=500, size=10):
    """
    A movie of each particle as a point
    """
    conf, track, pegs = load(base)

    fig = pl.figure(figsize=(size,size*conf['top']/conf['wall']))
    plot = None

    for t in xrange(1,max(frames, track.shape[1]/skip)):
        tmp = track[:,t*skip,:]
        if not ((tmp[:,0] > 0) & (tmp[:,1] > 0) & (tmp[:,0] < conf['wall']) & (tmp[:,1] < conf['top'])).any():
            continue

        if plot is None:
            plot = pl.plot(tmp[:,0], tmp[:,1], 'k,', alpha=1.0, ms=0.1)[0]
            pl.xticks([])
            pl.yticks([])
            pl.xlim(0,conf['wall'])
            pl.ylim(0,conf['top'])
            pl.tight_layout()
        else:
            plot.set_xdata(tmp[:,0])
            plot.set_ydata(tmp[:,1])
        pl.draw()
        pl.savefig(base+'-movie-%05d.png' % (t-1))

	def test_seq(self):
		sleep(1)
		
		self.cmd_logger(0)
		self.takeoffpub.publish(Empty())	;print 'takeoff' #takeoff
		sleep(12); print '4'; sleep(1); print '3'; sleep(1); print '2'; sleep(1); print '1'; sleep(1)
		
		self.cmd_logger(0)
		self.twist.linear.z = self.vzcmd
		self.cmd_logger(self.vzcmd)
		self.cmdpub.publish(self.twist)		;print 'vzcmd' #set vzcmd
		sleep(self.vzdur)			#wait for vzdur
		
		self.cmd_logger(self.vzcmd)
		self.clear_twist()
		self.cmd_logger(0)
		self.cmdpub.publish(self.twist)		;print 'clear vz' #clear vz
		sleep(4)
		
		self.cmd_logger(0)
		self.landpub.publish(Empty())		;print 'land' #land
		sleep(1)
		
		if not raw_input('show and save?') == 'n':
			pl.xlabel('time (s)')
			pl.ylim(0,2400)
			pl.plot(self.cmd_log['tm'],self.cmd_log['cmd'], 'b-s')
			pl.plot(self.nd_log['tm'],self.nd_log['vz'], 'g-+')
			pl.plot(self.nd_log['tm'],self.nd_log['al'], 'r-+')
			pl.grid(True)
			pl.show()

def plot_file_color(base, thin=True, start=0, size=14, save=False):
    conf, track, pegs = load(base)

    fig = pl.figure(figsize=(size,size*conf['top']/conf['wall']))

    track = track[start:]
    x = track[:,0];   y = track[:,1]
    t = np.linspace(0,1,x.shape[0])
    points = np.array([x,y]).transpose().reshape(-1,1,2)
    segs = np.concatenate([points[:-1],points[1:]],axis=1)
    lc = LineCollection(segs, linewidths=0.25, cmap=pl.cm.coolwarm)
    lc.set_array(t)
    pl.gca().add_collection(lc)

    #pl.scatter(x, y, c=np.arange(len(x)),linestyle='-',cmap=pl.cm.coolwarm)
    #pl.plot(track[-1000000:,0], track[-1000000:,1], '-', linewidth=0.0125, alpha=0.8)
    for peg in pegs:
        pl.gca().add_artist(pl.Circle(peg, conf['radius'], color='k', alpha=0.3))
    pl.xlim(0, conf['wall'])
    pl.ylim(0, conf['top'])
    pl.xticks([])
    pl.yticks([])
    pl.tight_layout()
    pl.show()
    if save:
        pl.savefig(base+".png", dpi=200)

def plot_prob_effector(sens, fpr, xmax=1, baserate=0.1):
    """Plots a line graph of P(effector|positive test) against
    the baserate of effectors in the input set to the classifier.
        
    The baserate argument draws an annotation arrow
    indicating P(pos|+ve) at that baserate
    """
    assert 0.1 <= xmax <= 1, "Max x axis value must be in range [0,1]"
    assert 0.01 <= baserate <= 1, "Baserate annotation must be in range [0,1]"
    baserates = pylab.arange(0, 1.05, xmax * 0.005)  
    probs = [p_correct_given_pos(sens, fpr, b) for b in baserates]
    pylab.plot(baserates, probs, 'r')
    pylab.title("P(eff|pos) vs baserate; sens: %.2f, fpr: %.2f" % (sens, fpr))
    pylab.ylabel("P(effector|positive)")
    pylab.xlabel("effector baserate")
    pylab.xlim(0, xmax)
    pylab.ylim(0, 1)
    # Add annotation arrow
    xpos, ypos = (baserate, p_correct_given_pos(sens, fpr, baserate))
    if baserate < xmax:
        if xpos > 0.7 * xmax:
            xtextpos = 0.05 * xmax
        else:
            xtextpos = xpos + (xmax-xpos)/5.
        if ypos > 0.5:
            ytextpos = ypos - 0.05
        else:
            ytextpos = ypos + 0.05
        pylab.annotate('baserate: %.2f, P(pos|+ve): %.3f' % (xpos, ypos), 
                       xy=(xpos, ypos), 
                       xytext=(xtextpos, ytextpos),
                       arrowprops=dict(facecolor='black', shrink=0.05))
    else:
        pylab.text(0.05 * xmax, 0.95, 'baserate: %.2f, P(pos|+ve): %.3f' % \
                   (xpos, ypos))

def pseudoSystem():

	#The corrolations discovered when answering this question shows the emergent effects of component evolution on CSE.
	#We can further study these correlations by creating an example component system consisting of many components where a a different component has a new version released every day.
	#By looking at users who Upgrade the system at different frequencies over 100 days, we present two graphs, Upgrade frequency to uttd and change.	
	
	
	l = 100
	uttdxy = []
	chxy = []
	for uf in range(1,20):
		uttd = range(uf)*(l*2/uf)
		uttd = uttd[1:l+1]
		uttdxy.append((uf,numpy.mean(uttd)))
		
		sh = [0]*(uf-1) + [uf]
		sh = sh*l
		sh = sh[:l]
		chxy.append((uf,sum(sh)))
		
	pylab.figure(20)
	x,y = zip(*sorted(uttdxy))
	pylab.plot(x,y)
	pylab.scatter(x,y)
	saveFigure("q1bpseudouttd")
	
	pylab.figure(21)
	
	x,y = zip(*sorted(chxy))
	pylab.plot(x,numpy.array(y))
	pylab.scatter(x,numpy.array(y))
	pylab.ylim([0,l+10])
	saveFigure("q1bpseudochange")

def plotFeatureImportance(featureImportance, title, originalImage=None, lim=0.06, colorate=None):
    """
    originalImage : the index of the original image. If None, ignore
    """
    indices = featureImportanceIndices(len(featureImportance), originalImage)
    pl.figure()
    pl.title(title)
    if colorate is not None:
        nbType = len(colorate)
        X = [[] for i in range(nbType)]
        Y = [[] for i in range(nbType)]
        for j, f in enumerate(featureImportance):
            X[j % nbType].append(j)
            Y[j % nbType].append(f)
        for i in range(nbType):
            pl.bar(X[i], Y[i], align="center", label=colorate[i][0], color=colorate[i][1])
        pl.legend()
    else:
        pl.bar(range(len(featureImportance)), featureImportance, align="center")
    #pl.xticks(pl.arange(len(indices)), indices, rotation=-90)
    pl.xlim([-1, len(indices)])
    pl.ylabel("Feature importance")
    pl.xlabel("Filter indices")
    pl.ylim(0, lim)
    pl.show()

def test_anns(results):
    data = read_metrics(results)
    split = int(len(data) * 0.8)
    num_input = len(data[0]) - 1
    random.shuffle(data)
    test_data = data[split:]

    # Read ANNs
    ann_dir = '/home/tomas/Dropbox/Git/ga_sandbox/projects/denoising/neural/trained_anns'
    trained_anns = []
    for filename in os.listdir(ann_dir):
        if filename.endswith('.net'):
            ann_path = os.path.join(ann_dir, filename)
            ann = libfann.neural_net()
            ann.create_from_file(ann_path)
            trained_anns.append(ann)

    points = []
    for row in test_data:
        actual_output = row[num_input]
        ann_mean_output = np.mean([
            ann.run(row[:num_input])
            for ann in trained_anns
        ])
        points.append([ann_mean_output, actual_output])
        print "actual: " + str(actual_output) + ", predicted: " + str(ann_mean_output)

    points = sorted(points, key=lambda p: p[1])
    fig = plt.figure(figsize=(6, 4))
    plt.plot([p[0] for p in points])
    plt.plot([p[1] for p in points])
    plt.ylim([0, 1.2])
    plt.show()

def correlation_matrix(data, size=8.0):
    """ Calculates and shows the correlation matrix of the pandas data frame
        'data' as a heat map.
        Only the correlations between numerical variables are calculated!
    """
    # calculate the correlation matrix
    corr = data.corr()
    #print corr
    lc = len(corr.columns)
    # set some settings for plottin'
    pl.pcolor(corr, vmin = -1, vmax = 1, edgecolor = "black")
    pl.colorbar()
    pl.xlim([-5,lc])
    pl.ylim([0,lc+5])
    pl.axis('off')
    # anotate the rows and columns with their corresponding variables
    ax = pl.gca()            
    for i in range(0,lc):
        ax.annotate(corr.columns[i], (-0.5, i+0.5), \
            size='large', horizontalalignment='right', verticalalignment='center')
        ax.annotate(corr.columns[i], (i+0.5, lc+0.5),\
            size='large', rotation='vertical',\
            horizontalalignment='center', verticalalignment='right')
    # change the size of the image
    fig = pl.figure(num=1)    
    fig.set_size_inches(size+(size/4), size)     
    
    pl.show()

    def plotMean(self, Channel, listPops, colors = [], ylabel = "Fluorescence (a.u.)"):
        """
        missing doc
        """

        if type(listPops) != type([]): listPops = [listPops]

        maxf = 0.
        minf = 1.e10
        maxx = 0.0
        minx = -1.0

        for pop in listPops:
            F = np.array( self.getFluorescence(pop, Channel).mean(axis=1) )
            plot.simplePlot( F, fillcolor=self.colors[pop], label = pop )

            if maxf < F.max().max() : maxf = F.max().max()
            if minf > F.min().min() : minf = F.min().min()

            maxx = max(maxx, F.shape[0])


        # setting labels and axes
        pl.xlabel('Time (h)')
        pl.ylabel(ylabel)
        pl.ylim(0.3*minf, 1.2*maxf)
        pl.xlim(minx, maxx)

        #pl.tight_layout()

        return

    def plot_roc(self, roc=None):
        """Plot ROC curves

        .. plot::
            :include-source:
            :width: 80%

            from dreamtools import rocs
            r = rocs.ROC()
            r.scores = [.9,.5,.6,.7,.1,.2,.6,.4,.7,.9, .2]
            r.classes = [1,0,1,0,0,1,1,0,0,1,1]
            r.plot_roc()

        """
        if roc == None:
            roc = self.get_roc()
        from pylab import plot, xlim, ylim ,grid, title, xlabel, ylabel
        x = roc['fpr']
        plot(x, roc['tpr'], '-o')
        plot([0,1], [0,1],'r')
        ylim([0, 1])
        xlim([0, 1])
        grid(True)
        title("ROC curve (AUC=%s)" % self.compute_auc(roc))
        xlabel("FPR")
        ylabel("TPR")

def makeimg(wav):
	global callpath
	global imgpath

	fs, frames = wavfile.read(os.path.join(callpath, wav))
	
	pylab.ion()

	# generate specgram
	pylab.figure(1)
	
	# generate specgram
	pylab.specgram(
		frames,
		NFFT=256, 
		Fs=22050, 
		detrend=pylab.detrend_none,
		window=numpy.hamming(256),
		noverlap=192,
		cmap=pylab.get_cmap('Greys'))
	
	x_width = len(frames)/fs
	
	pylab.ylim([0,11025])
	pylab.xlim([0,round(x_width,3)-0.006])
	
	img_path = os.path.join(imgpath, wav.replace(".wav",".png"))

	pylab.savefig(img_path)
	
	return img_path