def plotuttdpc():
	pylab.figure(1)

	auttdpc = uttdperc(always)
	modauttdpc = uttdperc(modalways)
	
	for name,pf,c in variables: 
		ivals = map(lambda x : uttdperc(x),pf)
		imean,istd,imeanpstd,imeanmstd = multimeanstd(ivals)
		pylab.fill_between(pallthedays, imeanpstd, imeanmstd, facecolor=c, alpha=0.3)
	
	for name,pf,c in variables: 
		ivals = map(lambda x : uttdperc(x),pf)
		imean,istd,imeanpstd,imeanmstd = multimeanstd(ivals)
		mdiff = numpy.mean(imean)
		pylab.plot(pallthedays,imean,color=c,label=("Mean (+-1std) UTTDpC of 30 \"%s\" users" % name))

	pylab.plot(pallthedays,auttdpc,color='black',label="UTTDpC of \"Always Upgrade\" user")
	pylab.plot(pallthedays,modauttdpc,color='red',label="UTTDpC of \"Progressive Always Upgrade\" user")
	
	print "Last uttd always",auttdpc[-1]
	print "Last uttd mod always",modauttdpc[-1]
	
	pylab.legend(loc="upper left")
	pylab.xlabel("Date")
	pylab.ylabel("Uptodate Distance per Component")
	pylab.title("Uptodate Distance per Component of Users")
	pylab.ylim([0,1])
	saveFigure("q4auttdperc")

def plotnew():
	fig = pylab.figure(20)
	
	
	
	for name,pf,c in variables: 
		ivals = map(lambda x : nntt(x),pf)
		imean,istd,imeanpstd,imeanmstd = multimeanstd(ivals)
		pylab.fill_between(pallthedays, imeanpstd, imeanmstd, facecolor=c, alpha=0.3)
	
	for name,pf,c in variables: 
		ivals = map(lambda x : nntt(x),pf)
		imean,istd,imeanpstd,imeanmstd = multimeanstd(ivals)
		mdiff = numpy.mean(imean)
		pylab.plot(pallthedays,imean,color=c,label=(name +"+-1std "))
	
	nntal = nntt(always)
	nntmod =nntt(modalways) 
	pylab.plot(pallthedays,nntal, color="black", label="Always Upgrade Mean change")
	pylab.plot(pallthedays,nntmod, color="red", label="Always Upgrade Mean change")
	
	print "Last new always",nntal[-1]
	print "Last new mod always",nntmod[-1]
	
	pylab.legend(loc="upper left")

	saveFigure("q4anew")

def plotchange():
	fig = pylab.figure(10)
	
	
	chtalw = chtt(always)
	chtmodalw= chtt(modalways)
	for name,pf,c in variables: 
		ivals = map(lambda x : chtt(x),pf)
		imean,istd,imeanpstd,imeanmstd = multimeanstd(ivals)
		pylab.fill_between(pallthedays, imeanpstd, imeanmstd, facecolor=c, alpha=0.3)
	
	for name,pf,c in variables: 
		ivals = map(lambda x : chtt(x),pf)
		imean,istd,imeanpstd,imeanmstd = multimeanstd(ivals)
		mdiff = numpy.mean(imean)
		pylab.plot(pallthedays,imean,color=c,label=("Mean (+-1std) Total Change of 30 \"%s\" users" % name))
	
	pylab.plot(pallthedays,chtalw, color="black", label="Total Change of \"Always Upgrade\" user")
	pylab.plot(pallthedays,chtmodalw, color="red", label="Total Change of \"Progressive Always Upgrade\" user")
	
	print "Last change always",chtalw[-1]
	print "Last change mod always",chtmodalw[-1]
	
	pylab.legend(loc="upper left")
	pylab.xlabel("Date")
	pylab.ylabel("Total Change")
	pylab.title("Total Change of Users")
	
	saveFigure("q4achange")

def Draw(func1, func2):
    # генирация точек графика
    xlist = mlab.frange(a, b, 0.01)
    ylist = [func1(x) for x in xlist]
    ylist2 = [func2(x) for x in xlist]

    # Генирирум ось
    y0 = [0 for x in xlist]
    pylab.plot(xlist, ylist)
    #pylab.plot(xlist, y0, label='line1', color='blue')
    pylab.plot(xlist, ylist2, label='$sin(x)/x)$', color='red')
    pylab.legend()

    # Включаем рисование сетки
    pylab.grid(True)

    pylab.fill_between(xlist, ylist, ylist2, color='green', alpha=0.25)
    # если мало разбиений, то переопереляем сетку под шаг
    if ((round((b - a) / h)) < 25):
        pylab.xticks([a + i * h for i in range(round((b - a) / h) + 1)])
    # рисуем корни, промерка того что корень не содержит ошибок
    for i in range(1, len(table)):
        if (table[i][4] != ':-('):
            pylab.scatter(table[i][3], table[i][4])

    # Рисуем фогрму с графиком
    pylab.show()

def shade_bands(edges, y_range=[-1e5,1e5],cmap='prism', **kwargs):
    '''
    Shades frequency bands.
    
    when plotting data over a set of frequency bands it is nice to 
    have each band visually seperated from the other. The kwarg `alpha`
    is useful.
    
    Parameters 
    --------------
    edges : array-like
        x-values seperating regions of a given shade
    y_range : tuple 
        y-values to shade in 
    cmap : str
        see matplotlib.cm  or matplotlib.colormaps for acceptable values
    \*\* : key word arguments
        passed to `matplotlib.fill_between`
        
    Examples 
    -----------
    >>> rf.shade_bands([325,500,750,1100], alpha=.2)
    '''
    cmap = plb.cm.get_cmap(cmap)
    for k in range(len(edges)-1):
        plb.fill_between(
            [edges[k],edges[k+1]], 
            y_range[0], y_range[1], 
            color = cmap(1.0*k/len(edges)),
            **kwargs)

 def plot(self):
     if not self.plot_state: return
     pop,best = self.plot_state
     with self.pylab_interface:
         import pylab
         pylab.clf()
         n,p = pop.shape
         iternum = numpy.arange(1,n+1)
         tail = int(0.25*n)
         pylab.hold(True)
         c = coordinated_colors(base=(0.4,0.8,0.2))
         if p==5:
             pylab.fill_between(iternum[tail:], pop[tail:,1], pop[tail:,3],
                                color=c['light'], label='_nolegend_')
             pylab.plot(iternum[tail:],pop[tail:,2],
                        label="80% range", color=c['base'])
             pylab.plot(iternum[tail:],pop[tail:,0],
                        label="_nolegend_", color=c['base'])
         else:
             pylab.plot(iternum,pop, label="population",
                        color=c['base'])
         pylab.plot(iternum[tail:], best[tail:], label="best",
                    color=c['dark'])
         pylab.xlabel('iteration number')
         pylab.ylabel('chisq')
         pylab.legend()
         #pylab.gca().set_yscale('log')
         pylab.hold(False)
         pylab.draw()

 def view_simple( self, stats, thetas ):
   # plotting params
   nbins       = 20
   alpha       = 0.5
   label_size  = 8
   linewidth   = 3
   linecolor   = "r"
   
   # extract from states
   #thetas = states_object.get_thetas()[burnin:,:]
   #stats  = states_object.get_statistics()[burnin:,:]
   #nsims  = states_object.get_sim_calls()[burnin:]
   
   # plot sample distribution of thetas, add vertical line for true theta, theta_star
   f = pp.figure()
   sp = f.add_subplot(111)
   pp.plot( self.fine_theta_range, self.posterior, linecolor+"-", lw = 1)
   ax = pp.axis()
   pp.hist( thetas, self.nbins_coarse, range=self.range,normed = True, alpha = alpha )
   
   pp.fill_between( self.fine_theta_range, self.posterior, color="m", alpha=0.5)
   
   pp.plot( self.posterior_bars_range, self.posterior_bars, 'ro')
   pp.vlines( thetas.mean(), ax[2], ax[3], color="b", linewidths=linewidth)
   #pp.vlines( self.theta_star, ax[2], ax[3], color=linecolor, linewidths=linewidth )
   pp.vlines( self.posterior_mode, ax[2], ax[3], color=linecolor, linewidths=linewidth )
   
   pp.xlabel( "theta" )
   pp.ylabel( "P(theta)" )
   pp.axis([self.range[0],self.range[1],ax[2],ax[3]])
   set_label_fonsize( sp, label_size )
   pp.show()

def plotBkMeasure(bk, ek, vk, figurePath):
    #print bk
    #print ek
    #print vk
    
    k = list(range(len(bk)))
    
    #for i,j in enumerate(bk):
    pylab.ioff()
    pylab.figure()
    pylab.plot(k, bk, '.', label='Bk')

    pylab.plot(k, ek, label='E(Bk)')
    #pylab.plot(k, ek+2*np.sqrt(vk), '-.r', label='limit range')
    #pylab.plot(k, ek-2*np.sqrt(vk), '-.r')
    #for i in range(len(ek)):
    pylab.fill_between(k, ek+4*np.sqrt(vk), ek-4*np.sqrt(vk), facecolor='red', interpolate=True )

    # figure setting
    pylab.xlim(2,k[-1])
    pylab.ylim(0,1.0)
    pylab.legend(loc='upper right')
    pylab.xlabel('Number of Clusters')
    pylab.ylabel('Bk')
    # pylab.title('Bk measure between two algorithm')

    # show result
    pylab.savefig(figurePath, format='svg')

def sampleplot_K(r,ylim=None,HDI_y=None):
    from pylab import plot,fill_between,gca,text
    
    x,y=histogram(r,plot=False)
    
    plot(x,y,'-o')
    
    fill_between(x,y,facecolor='blue', alpha=0.2)
    if ylim:
        gca().set_ylim(ylim)

    dx=x[1]-x[0]
    cs=np.cumsum(y)*dx
    
    HDI=np.percentile(r,[2.5,50,97.5])

    yl=gca().get_ylim()
    
    dy=0.05*yl[1]
    if HDI_y is None:
        HDI_y=yl[1]*.1
        
    
    text((HDI[0]+HDI[2])/2, HDI_y+dy,'95% HDI', ha='center', va='center',fontsize=12)
    plot(HDI,[HDI_y,HDI_y,HDI_y],'k.-',linewidth=1)
    for v in HDI:
        text(v, HDI_y-dy,'%.3f' % v, ha='center', va='center', 
             fontsize=12)
    xl=gca().get_xlim()
    
    text(.05*(xl[1]-xl[0])+xl[0], 0.9*yl[1],r'$\tilde{x}=%.3f$' % np.median(r), ha='left', va='center')

    def bootstrap(self, nBoot, nbins = 20):
        pops = np.zeros((nBoot, nbins))
        #medianpop = [[] for i in data.cat]
        pylab.figure(figsize = (20,14))
        for i in xrange(3):
            pylab.subplot(1,3,i+1)
            #if  i ==0:
                #pylab.title("Bootstrap on medians", fontsize = 20.)
            pop = self.angles[(self.categories == i)]# & (self.GFP > 2000)]
            for index in xrange(nBoot):
                newpop = np.random.choice(pop, size=len(pop), replace=True)
                #medianpop[i].append(np.median(newpop))
                newhist, binedges = np.histogram(newpop, bins = nbins)
                pops[index,:] = newhist/1./len(pop)
            #pylab.hist(medianpop[i], bins = nbins, label = "{2} median {0:.1f}, std {1:.1f}".format(np.median(medianpop[i]), np.std(medianpop[i]), data.cat[i]), color = data.colors[i], alpha =.2, normed = True)

            meanpop = np.sum(pops, axis = 0)/1./nBoot
            stdY = np.std(pops, axis = 0)
            print "width", binedges[1] - binedges[0]
            pylab.bar(binedges[:-1], meanpop, width = binedges[1] - binedges[0], label = "mean distribution", color = data.colors[i], alpha = 0.6)
            pylab.fill_between((binedges[:-1]+binedges[1:])/2., meanpop-stdY, meanpop+stdY, alpha = 0.3)
            pylab.legend()
            pylab.title(data.cat[i])
            pylab.xlabel("Angle(degree)", fontsize = 15)
            pylab.ylim([-.01, 0.23])

        pylab.savefig("/users/biocomp/frose/frose/Graphics/FINALRESULTS-diff-f3/distrib_nBootstrap{0}_bins{1}_GFPsup{2}_{3}.png".format(nBoot, nbins, 'all', randint(0,999)))

def plotWithVariance(x, y, variance, *args, **kwargs):
    """
    Plot data with variance indicated by shading within one sigma.
    """
    line = pylab.plot(x, y.flatten(), *args, **kwargs)[0]
    sigma = np.sqrt(variance)
    pylab.fill_between(x, y - sigma, y + sigma, color=line.get_color(), alpha=0.5)

def drawROC(points,zeTitle,zeFilename,visible,show_fig,save_fig=True,
            special_point=None,special_value=None,special_label=None):
    AUC=computeAUC(points)
    import pylab

    pylab.clf()
    pylab.grid(color='#aaaaaa', linestyle='-', linewidth=1,alpha=0.5)

    pylab.plot([x[0] for x in points], [y[1] for y in points], '-', linewidth=3,color="#000088",zorder=3)
    pylab.fill_between([x[0] for x in points], [y[1] for y in points],0,color='0.9')
    pylab.plot([0.0,1.0], [0.0, 1.0], '-',color="#AAAAAA")

    pylab.ylim((-0.01,1.01))
    pylab.xlim((-0.01,1.01))
    pylab.xticks(pylab.arange(0,1.1,.1))
    pylab.yticks(pylab.arange(0,1.1,.1))
    pylab.grid(True)

    ax=pylab.gca()
    r = pylab.Rectangle((0,0), 1, 1, edgecolor='#444444', facecolor='none',zorder=1)
    ax.add_patch(r)
    [spine.set_visible(False) for spine in ax.spines.values()]

    if len(points)<10:
      for i in range(1,len(points)-1):
        pylab.plot(points[i][0],points[i][1],'o',color="#000066",zorder=6)

    pylab.xlabel('False positive rate')
    pylab.ylabel('True positive rate')

    if special_point is not None:
        pylab.plot(special_point[0],special_point[1],'o',color="#DD9999",zorder=6)
        if special_value is not None:
            pylab.text(special_point[0]+0.01,special_point[1]-0.01, special_value,
                       {'color' : '#DD5555', 'fontsize' : 10},
                       horizontalalignment = 'left',
                       verticalalignment = 'top',
                       rotation = 0,
                       clip_on = False)
    if special_label is not None:
        if special_label!="":
            labels=[special_label]
            colors=['#DD9999']
            circles=[pylab.Circle((0, 0), 1, fc=colors[0])]
            legend_location = 'lower right'
            pylab.legend(circles, labels, loc=legend_location)

    pylab.text(0.5, 0.3,'AUC=%f'%AUC,
     horizontalalignment='center',
     verticalalignment='center',
     fontsize=18)

    pylab.title(zeTitle)

    if save_fig:
        pylab.savefig(zeFilename,dpi=300)
        print("\n result in "+zeFilename)

    if show_fig:
        pylab.show()

def visualize(generation_list):
    '''Generate pretty pictures using pylab and pygame'''
    best = []
    average = []
    stddev = []
    average_plus_stddev = []
    average_minus_stddev = []
    for pop in generation_list:
        best += [most_fit(pop).fitness]
        average += [avg_fitness(pop)]
        stddev += [fitness_stddev(pop)] 
        average_plus_stddev += [average[-1] + stddev[-1]]
        average_minus_stddev += [average[-1] - stddev[-1]]
    
    pylab.figure(1)
    pylab.fill_between(range(len(generation_list)), average_plus_stddev, average_minus_stddev, alpha=0.2, color='b', label="Standard deviation")
    pylab.plot(range(len(generation_list)), best, color='r', label='Best')
    pylab.plot(range(len(generation_list)), average, color='b', label='Average with std.dev.')
    pylab.title("Fitness plot - Beer-cog")
    pylab.xlabel("Generation")
    pylab.ylabel("Fitness")
    pylab.legend(loc="upper left")
    pylab.savefig("mincog_fitness.png")

    best_index = best.index(max(best))
    best_individual = most_fit(generation_list[-1])

    with open('last.txt','w') as f:
        f.write(str(best_individual.gtype))
    print best_individual.gtype
    
    game = min_cog_game.Game()
    game.play(best_individual.ptype, True)

def show_barlines(page):
    import pylab
    for i, barlines in enumerate(page.barlines):
        sd = page.staves.staff_dist[i]
        for j, barline_range in enumerate(barlines):
            barline_x = int(barline_range.mean())
            staff_y = page.staves.staff_y(i, barline_x)
            repeats = page.repeats[i][j]
            if repeats:
                # Draw thick bar
                pylab.fill_between([barline_x - sd/4,
                                    barline_x + sd/4],
                                   staff_y - sd*2,
                                   staff_y + sd*2,
                                   color='g')
                for letter, sign in (('L', -1), ('R', +1)):
                    if letter in repeats:
                        # Draw thin bar
                        bar_x = barline_x + sign * sd/2
                        pylab.plot([bar_x, bar_x],
                                   [staff_y - sd*2,
                                    staff_y + sd*2],
                                   color='g')
                        for y in (-1, +1):
                            circ = pylab.Circle((bar_x + sign*sd/2,
                                                 staff_y + y*sd/2),
                                                sd/4,
                                                color='g')
                            pylab.gcf().gca().add_artist(circ)
            else:
                pylab.plot([barline_x, barline_x],
                           [staff_y - sd*2,
                            staff_y + sd*2],
                           color='g')

def addqqplotinfo(qnull,M,xl='-log10(P) observed',yl='-log10(P) expected',xlim=None,ylim=None,alphalevel=0.05,legendlist=None,fixaxes=False):    
    distr='log10'
    pl.plot([0,qnull.max()], [0,qnull.max()],'k')
    pl.ylabel(xl)
    pl.xlabel(yl)
    if xlim is not None:
        pl.xlim(xlim)
    if ylim is not None:
        pl.ylim(ylim)        
    if alphalevel is not None:
        if distr == 'log10':
            betaUp, betaDown, theoreticalPvals = _qqplot_bar(M=M,alphalevel=alphalevel,distr=distr)
            lower = -sp.log10(theoreticalPvals-betaDown)
            upper = -sp.log10(theoreticalPvals+betaUp)
            pl.fill_between(-sp.log10(theoreticalPvals),lower,upper,color="grey",alpha=0.5)
            #pl.plot(-sp.log10(theoreticalPvals),lower,'g-.')
            #pl.plot(-sp.log10(theoreticalPvals),upper,'g-.')
    if legendlist is not None:
        leg = pl.legend(legendlist, loc=4, numpoints=1)
        # set the markersize for the legend
        for lo in leg.legendHandles:
            lo.set_markersize(10)

    if fixaxes:
        fix_axes()        

def GP_plotpred(xpred, x, y, cov_par, cov_func = None, cov_typ = 'SE',
             MF = None, MF_par = None, MF_args = None, MF_args_pred = None, \
             WhiteNoise = False, plot_color = None):
    '''
    Wrapper for GP_predict that takes care of merging the
    covariance and mean function parameters, and (optionally) plots
    the predictive distribution (as well as returning it)
    '''
    if MF != None:
        merged_par = scipy.append(cov_par, MF_par)
        n_MF_par = len(MF_par)
    else:
        merged_par = cov_par[:]
        n_MF_par = 0
    fpred, fpred_err = GP_predict(merged_par, xpred, x, y, \
                                  cov_func = cov_func, cov_typ = cov_typ, \
                                  MF = MF, n_MF_par = n_MF_par, \
                                  MF_args = MF_args, MF_args_pred = MF_args_pred, \
                                  WhiteNoise = WhiteNoise)
    xpl = scipy.array(xpred[:,0]).flatten()
    if plot_color != None:
        pylab.fill_between(xpl, fpred + 2 * fpred_err, fpred - 2 * fpred_err, \
                           color = plot_color, alpha = 0.1)
        pylab.fill_between(xpl, fpred + fpred_err, fpred - fpred_err, \
                           color = plot_color, alpha = 0.1)
        pylab.plot(xpl, fpred, '-', color = plot_color)
    return fpred, fpred_err

def plot_shaded_lines(my_xticks, y1, y2, error1, error2, ylab, xlab, filename):
    plt.figure(figsize=(6,6))
    from matplotlib import rcParams
    rcParams.update({'figure.autolayout': True})

    x = range(0, len(y1))
    plt.plot(x, y1, 'k-', color="blue",  label='men')
    plt.fill_between(x, y1-error1, y1+error1, facecolor='blue', alpha=.2)

    plt.plot(x, y2, 'k-', color="red",  label='women')
    plt.fill_between(x, y2-error2, y2+error2, facecolor='red', alpha=.2)

    #if isinstance(x, (int, long, float, complex)):
    #    plt.xlim(np.min(x), np.max(x))
    plt.gcf().subplots_adjust(bottom=0.3)

    plt.xticks(x, my_xticks)
    plt.xticks(rotation=70, fontsize=14)
    plt.yticks(fontsize=14)
    #plt.setp(ax.get_xticklabels(), rotation='vertical', fontsize=14)
    plt.ylabel(ylab, fontsize=14)
    plt.xlabel(xlab, fontsize=14)
    plt.legend()

    plt.savefig(filename)

def sun():
    dur = 1000
    cad = 29.4 / 60.0 / 24.0
    X = numpy.genfromtxt('%ssun/sun_composite_tsi_20130930.txt' % root).T
    date = X[0]
    t = X[1]
    irr = X[2]
    l = (irr > 0) * (t > 5000)
    t = t[l] - t[l].min()
    date_ref = date[l][0]
    irr = irr[l] / irr[l].max()
    pylab.figure(1)
    pylab.clf()
    pylab.plot(t, irr, 'k-')
    col = ['r','g','b','y','m']
    tstart = [1000, 2100, 2600, 3800, 5000]
    for i in numpy.arange(5):
        time = numpy.r_[tstart[i]:tstart[i]+dur:cad]
        if i == 0:
            y1 = numpy.zeros(len(time)) + 1
            y2 = y1 - 0.0045
        g = scipy.interpolate.interp1d(t, irr, bounds_error = False)
        dF = filter.boxcare(g(time), 10, fill = True)
        pylab.plot(time, dF, c = col[i])
        pylab.fill_between(time, y1, y2, color = col[i], alpha = 0.3)
        X = numpy.zeros((2,len(time)))
        X[0,:] = time - time[0]
        X[1,:] = dF
        numpy.savetxt('%ssun/sun_lightcurve_%02d.txt' % (root, i), X.T)
    pylab.xlim(t.min(), t.max())
    pylab.ylim(0.9955,1)
    pylab.xlabel('Days since %d' % date_ref)
    pylab.ylabel('Flux decrement')
    pylab.savefig('%ssun/sun_lightcurves.png' % root)
    return

def plotchange():
	fig = pylab.figure(10)
	
	
	for name,pf,c in variables: 
		ivals = map(lambda x : chtt(x),pf)
		imean,istd,imeanpstd,imeanmstd = multimeanstd(ivals)
		pylab.fill_between(pallthedays, imeanpstd, imeanmstd, facecolor=c, alpha=0.3)
	
	for name,pf,c in variables: 
		ivals = map(lambda x : chtt(x),pf)
		imean,istd,imeanpstd,imeanmstd = multimeanstd(ivals)
		mdiff = numpy.mean(imean)
		print "Final change",name,imean[-1]
		pylab.plot(pallthedays,imean,color=c,label=("Mean (+-1std) Total Change of 30 \"%s\" users" % name))
	
	pylab.legend(loc="upper left")
	pylab.xlabel("Date")
	pylab.ylabel("Total Change")
	pylab.title("Total Change of Users")
	pylab.ylim([0,1500])
	saveFigure("q1cchange")
	#This graph shows the range over which change can occur whn altering the probability to install.
	# given the results from this supports the idea that the probability to install corrolates to the amount of packages installed, this is straight forward.
	#An effect that did not occur, which was hypothesised might, was that over time the amount of installed packages decreases as common dependended packages are installed.
	#This does not occur, but may be because of the randomness that packages are selected.
	#It may still occur in reaility, as say a graphics designer will likely install graphics components that will require similar functionality.
	
	#Another interesting point from that can be seen in this graphi is the variability created by the soft failures, as discussed previously.
	# it can be seen that the std increased after a soft failure in install twice a week change curve.
	uivals = zip(alll,map(lambda x : rempd(x),alll))
	for name,ui in uivals:
		for date, remv in zip(allthedays,ui):
			if remv >= 100:
				print date,datetime.date.fromtimestamp(date),name

def demoplot(theta,args):
    colour=np.array([0,0,1.0])
    faded = 1-(1-colour)/2.0

    (X,y) = args
    (n, D) = np.shape(X)

    xrange = X.max() - X.min()
    Xtest = np.arange(X.min()-xrange/2,X.max()+xrange/2,(X.max()-X.min())/100)
    Xtest.shape = (len(Xtest),1)

    k = kernel2(X,X,theta,wantderiv=False)
    kstar = [kernel2(X,xs*np.ones((1,1)),theta,wantderiv=False,measnoise=False) for xs in Xtest]
    kstar = np.squeeze(kstar)
    kstarstar = [kernel2(xs*np.ones((1,1)),xs*np.ones((1,1)),theta,wantderiv=False,measnoise=False) for xs in Xtest]
    kstarstar = np.squeeze(kstarstar)

    L = np.linalg.cholesky(k)
    invk = np.linalg.solve(L.transpose(),np.linalg.solve(L,np.eye(np.shape(X)[0])))
    mean = np.dot(kstar,np.dot(invk,y))
    var = kstarstar - np.diag(np.dot(kstar,np.dot(invk,kstar.T)))
    #var = np.reshape(var,(100,1))

    pl.ion()
    fig = pl.figure()
    #ax1 = fig.add_subplot(211)
    #ax2 = fig.add_subplot(212,sharex=ax1,sharey=ax1)
    pl.plot(Xtest,mean,'-k')
    #pl.plot(xstar,mean+2*np.sqrt(var),'x-')
    #pl.plot(xstar,mean-2*np.sqrt(var),'x-')
    #print np.shape(xstar), np.shape(mean), np.shape(var)
    pl.fill_between(np.squeeze(Xtest),np.squeeze(mean-2*np.sqrt(var)),np.squeeze(mean+2*np.sqrt(var)),color='0.75')
    pl.plot(X,y,'ko')