python类xlim()的实例源码

def run_regression_1D_stoc():
    np.random.seed(42)

    print "create dataset ..."
    N = 200
    X = np.random.rand(N, 1)
    Y = np.sin(12 * X) + 0.5 * np.cos(25 * X) + np.random.randn(N, 1) * 0.2
    # plt.plot(X, Y, 'kx', mew=2)

    def plot(m):
        xx = np.linspace(-0.5, 1.5, 100)[:, None]
        mean, var = m.predict_f(xx)
        zu = m.sgp_layers[0].zu
        mean_u, var_u = m.predict_f(zu)
        plt.figure()
        plt.plot(X, Y, 'kx', mew=2)
        plt.plot(xx, mean, 'b', lw=2)
        plt.fill_between(
            xx[:, 0],
            mean[:, 0] - 2 * np.sqrt(var[:, 0]),
            mean[:, 0] + 2 * np.sqrt(var[:, 0]),
            color='blue', alpha=0.2)
        plt.errorbar(zu, mean_u, yerr=2 * np.sqrt(var_u), fmt='ro')
        plt.xlim(-0.1, 1.1)

    # inference
    print "create model and optimize ..."
    M = 20
    hidden_size = [2]
    model = aep.SDGPR(X, Y, M, hidden_size, lik='Gaussian')
    model.optimise(method='adam', alpha=1.0,
                   maxiter=50000, mb_size=M, adam_lr=0.001)
    plot(model)
    plt.show()
    plt.savefig('/tmp/aep_dgpr_1D_stoc.pdf')

def plot_q(model='cem', r_min=0.0, r_max=6371.0, dr=1.0):

    """
    Plot a radiallysymmetric Q model.

    plot_q(model='cem', r_min=0.0, r_max=6371.0, dr=1.0):

    r_min=minimum radius [km], r_max=maximum radius [km], dr=radius increment [km]

    Currently available models (model): cem, prem, ql6
    """

    r = np.arange(r_min, r_max+dr, dr)
    q = np.zeros(len(r))

    for k in range(len(r)):

        if model=='cem':
            q[k]=q_cem(r[k])
        elif model=='ql6':
            q[k]=q_ql6(r[k])
        elif model=='prem':
            q[k]=q_prem(r[k])


    plt.plot(r,q,'k')
    plt.xlim((0.0,r_max))
    plt.xlabel('radius [km]')
    plt.ylabel('Q')
    plt.show()


###################################################################################################
#- CEM, EUMOD
###################################################################################################

def plot_pr(auc_score, name, phase, precision, recall, label=None):
    pylab.clf()
    pylab.figure(num=None, figsize=(5, 4))
    pylab.grid(True)
    pylab.fill_between(recall, precision, alpha=0.5)
    pylab.plot(recall, precision, lw=1)
    pylab.xlim([0.0, 1.0])
    pylab.ylim([0.0, 1.0])
    pylab.xlabel('Recall')
    pylab.ylabel('Precision')
    pylab.title('P/R curve (AUC=%0.2f) / %s' % (auc_score, label))
    filename = name.replace(" ", "_")
    pylab.savefig(os.path.join(CHART_DIR, "pr_%s_%s.png" %
                  (filename, phase)), bbox_inches="tight")

def plot_roc(auc_score, name, fpr, tpr):
    pylab.figure(num=None, figsize=(6, 5))
    pylab.plot([0, 1], [0, 1], 'k--')
    pylab.xlim([0.0, 1.0])
    pylab.ylim([0.0, 1.0])
    pylab.xlabel('False Positive Rate')
    pylab.ylabel('True Positive Rate')
    pylab.title('Receiver operating characteristic (AUC=%0.2f)\n%s' % (
        auc_score, name))
    pylab.legend(loc="lower right")
    pylab.grid(True, linestyle='-', color='0.75')
    pylab.fill_between(tpr, fpr, alpha=0.5)
    pylab.plot(fpr, tpr, lw=1)
    pylab.savefig(
        os.path.join(CHART_DIR, "roc_" + name.replace(" ", "_") + ".png"))

def plot_pr(auc_score, name, precision, recall, label=None):
    pylab.figure(num=None, figsize=(6, 5))
    pylab.xlim([0.0, 1.0])
    pylab.ylim([0.0, 1.0])
    pylab.xlabel('Recall')
    pylab.ylabel('Precision')
    pylab.title('P/R (AUC=%0.2f) / %s' % (auc_score, label))
    pylab.fill_between(recall, precision, alpha=0.5)
    pylab.grid(True, linestyle='-', color='0.75')
    pylab.plot(recall, precision, lw=1)
    filename = name.replace(" ", "_")
    pylab.savefig(os.path.join(CHART_DIR, "pr_" + filename + ".png"))

def plot_roc(auc_score, name, fpr, tpr):
    pylab.figure(num=None, figsize=(6, 5))
    pylab.plot([0, 1], [0, 1], 'k--')
    pylab.xlim([0.0, 1.0])
    pylab.ylim([0.0, 1.0])
    pylab.xlabel('False Positive Rate')
    pylab.ylabel('True Positive Rate')
    pylab.title('Receiver operating characteristic (AUC=%0.2f)\n%s' % (
        auc_score, name))
    pylab.legend(loc="lower right")
    pylab.grid(True, linestyle='-', color='0.75')
    pylab.fill_between(tpr, fpr, alpha=0.5)
    pylab.plot(fpr, tpr, lw=1)
    pylab.savefig(os.path.join(CHART_DIR, "roc_" + name.replace(" ", "_")+ ".png"))

def plot_pr(auc_score, name, precision, recall, label=None):
    pylab.figure(num=None, figsize=(6, 5))
    pylab.xlim([0.0, 1.0])
    pylab.ylim([0.0, 1.0])
    pylab.xlabel('Recall')
    pylab.ylabel('Precision')
    pylab.title('P/R (AUC=%0.2f) / %s' % (auc_score, label))
    pylab.fill_between(recall, precision, alpha=0.5)
    pylab.grid(True, linestyle='-', color='0.75')
    pylab.plot(recall, precision, lw=1)
    filename = name.replace(" ", "_")
    pylab.savefig(os.path.join(CHART_DIR, "pr_" + filename + ".png"))

def plot_roc(auc_score, name, fpr, tpr):
    pylab.figure(num=None, figsize=(6, 5))
    pylab.plot([0, 1], [0, 1], 'k--')
    pylab.xlim([0.0, 1.0])
    pylab.ylim([0.0, 1.0])
    pylab.xlabel('False Positive Rate')
    pylab.ylabel('True Positive Rate')
    pylab.title('Receiver operating characteristic (AUC=%0.2f)\n%s' % (
        auc_score, name))
    pylab.legend(loc="lower right")
    pylab.grid(True, linestyle='-', color='0.75')
    pylab.fill_between(tpr, fpr, alpha=0.5)
    pylab.plot(fpr, tpr, lw=1)
    pylab.savefig(os.path.join(CHART_DIR, "roc_" + name.replace(" ", "_")+ ".png"))

def plot_pr(auc_score, name, precision, recall, label=None):
    pylab.figure(num=None, figsize=(6, 5))
    pylab.xlim([0.0, 1.0])
    pylab.ylim([0.0, 1.0])
    pylab.xlabel('Recall')
    pylab.ylabel('Precision')
    pylab.title('P/R (AUC=%0.2f) / %s' % (auc_score, label))
    pylab.fill_between(recall, precision, alpha=0.5)
    pylab.grid(True, linestyle='-', color='0.75')
    pylab.plot(recall, precision, lw=1)
    filename = name.replace(" ", "_")
    pylab.savefig(os.path.join(CHART_DIR, "pr_" + filename + ".png"))

def plot_entropy_vs_rVals(self):
        if not doViz:
            self.skipTest("Required module matplotlib unavailable.")
        H = np.sum(calcRlogR(self.R), axis=1)
        Hnew_exact = np.sum(calcRlogR(self.Rnew_Exact), axis=1)
        Hnew_approx = np.sum(calcRlogR(self.Rnew_Approx), axis=1)

        rVals = self.rVals

        np.set_printoptions(precision=4, suppress=True)
        print ''
        print '--- rVals'
        print rVals[:3], rVals[-3:]

        print '--- R original'
        print self.R[:3]
        print self.R[-3:, :]

        print '--- R proposal'
        print self.Rnew_Exact[:3]
        print self.Rnew_Exact[-3:, :]

        pylab.plot(rVals, H, 'k-', label='H original')
        pylab.plot(rVals, Hnew_exact, 'b-', label='H proposal exact')
        pylab.plot(rVals, Hnew_approx, 'r-', label='H proposal approx')
        pylab.legend(loc='best')
        pylab.xlim([rVals.min() - .01, rVals.max() + .01])
        ybuf = 0.05 * H.max()
        pylab.ylim([ybuf, H.max() + ybuf])
        pylab.show(block=True)

def plotSpeedupFigure(AllInfo, maxWorker=1, **kwargs):
    pylab.figure(2)
    xs = AllInfo['nWorker']
    ts_mono = AllInfo['t_monolithic']

    xgrid = np.linspace(0, maxWorker + 0.1, 100)
    pylab.plot(xgrid, xgrid, 'y--', label='ideal parallel')

    for method in getMethodNames(**kwargs):
        speedupRatio = ts_mono / AllInfo['t_' + method]
        pylab.plot(xs, speedupRatio, 'o-',
                   label=method,
                   color=ColorMap[method],
                   markeredgecolor=ColorMap[method])

    pylab.xlim([-0.2, maxWorker + 0.5])
    pylab.ylim([0, maxWorker + 0.5])
    pylab.legend(loc='upper left')
    pylab.xlabel('Number of Workers')
    pylab.ylabel('Speedup over Monolithic')

    if kwargs['savefig']:
        title = 'BenchmarkPlot_%s_%s_minDur=%.2f_Speedup.eps'\
            % (platform.node(), kwargs['task'], kwargs['minSliceDuration'])
        pylab.savefig(title,
                      format='eps',
                      bbox_inches='tight',
                      pad_inches=0)

def plotBoundVsK(KVals=np.arange(1,50),
                 alpha=0.5,
                 gamma=10,
                 labels=None,
                 betaFunc='prior'):
    if labels is None:
        txtlabel = str(alpha)
        labels = [None, None]
    else:
        txtlabel = 'alpha\n' + str(alpha)
    exactVals = np.zeros(len(KVals))
    boundVals = np.zeros(len(KVals))
    for ii, K in enumerate(KVals):
        betaVec = 1.0/(1.0 + gamma) * np.ones(K+1)
        for k in range(1, K):
            betaVec[k] = betaVec[k] * (1 - np.sum(betaVec[:k]))
        betaVec[-1] = 1 - np.sum(betaVec[:-1])
        print betaVec
        assert np.allclose(betaVec.sum(), 1.0)
        exactVals[ii] = cDir_exact(alpha, betaVec)
        boundVals[ii] = cDir_surrogate(alpha, betaVec)
    assert np.all(exactVals >= boundVals)
    pylab.plot(KVals, exactVals,
               'k-', linewidth=LINEWIDTH, label=labels[0])
    pylab.plot(KVals, boundVals,
               'r--', linewidth=LINEWIDTH, label=labels[1])
    index = -1

    pylab.text(KVals[index]+.25, boundVals[index],
        txtlabel, fontsize=LEGENDSIZE-8)
    pylab.xlim([0, np.max(KVals)+7.5])
    pylab.gca().set_xticks([0, 10, 20, 30, 40, 50])
    pylab.xlabel("K", fontsize=FONTSIZE)
    pylab.ylabel("cDir function", fontsize=FONTSIZE)
    pylab.tick_params(axis='both', which='major', labelsize=TICKSIZE)

def makeFigure(hmmKappa=0):
    Data, trueResp = makeDataAndTrueResp()

    kemptyVals = np.asarray([0, 1, 2, 3.])
    ELBOVals = np.zeros_like(kemptyVals, dtype=np.float)

    # Iterate over the number of empty states (0, 1, 2, ...)
    for ii, kempty in enumerate(kemptyVals):
        resp = makeNewRespWithEmptyStates(trueResp, kempty)
        ELBOVals[ii] = resp2ELBO_HDPHMM(Data, resp, hmmKappa=hmmKappa)

    # Make largest value the one with kempty=0, to make plot look good
    ELBOVals -= ELBOVals[0]

    # Plot the results
    from matplotlib import pylab
    figH = pylab.figure(figsize=(6, 4))
    plotargs = dict(markersize=10, linewidth=3)
    pylab.plot(kemptyVals, ELBOVals, 'o--', label='HDP surrogate',
               color='b', markeredgecolor='b',
               **plotargs)
    pylab.xlabel('num. empty topics', fontsize=20)
    pylab.ylabel('change in ELBO', fontsize=20)
    B = 0.25
    pylab.xlim([-B, kemptyVals[-1] + B])
    pylab.xticks(kemptyVals)

    axH = pylab.gca()
    axH.tick_params(axis='both', which='major', labelsize=15)
    legH = pylab.legend(loc='upper left', prop={'size': 15})

    figH.subplots_adjust(bottom=0.16, left=0.2)
    pylab.show(block=True)

def plot_convergence_data(errors, label):
    error_line, = plt.plot(errors, label=label)
    # plt.xlim([0, 100])
    plt.xlim([100, 1000])
    plt.ylim([0, 100])
    plt.xlabel("Generation")
    plt.ylabel("Error")
    plt.legend(loc=0, prop={'size': 8})

def plot_prof_1(self, species, keystring, xlim1, xlim2, ylim1,
                    ylim2, symbol=None, show=False):
        '''
        Plot one species for cycle between xlim1 and xlim2 Only works
        with instances of se and mesa _profile.

        Parameters
        ----------
        species : list
            Which species to plot.
        keystring : string or integer
            Label that appears in the plot or in the case of se, a
            cycle.
        xlim1, xlim2 : integer or float
            Mass coordinate range.
        ylim1, ylim2 : integer or float
            Mass fraction coordinate range.
        symbol : string, optional
            Which symbol you want to use.  If None symbol is set to '-'.
            The default is None.
        show : boolean, optional
            Show the ploted graph.  The default is False.

        '''
        plotType=self._classTest()
        if plotType=='se':
            #tot_mass=self.se.get(keystring,'total_mass')
            tot_mass=self.se.get('mini')
            age=self.se.get(keystring,'age')
            mass=self.se.get(keystring,'mass')
            Xspecies=self.se.get(keystring,'iso_massf',species)

            mod=keystring
        elif plotType=='mesa_profile':
            tot_mass=self.header_attr['star_mass']
            age=self.header_attr['star_age']
            mass=self.get('mass')
            mod=self.header_attr['model_number']
            Xspecies=self.get(species)
        else:
            print('This method is not supported for '+str(self.__class__))
            return

        if symbol == None:
            symbol = '-'

        x,y=self._logarithm(Xspecies,mass,True,False,10)
        #print x
        pl.plot(y,x,symbol,label=str(species))
        pl.xlim(xlim1,xlim2)
        pl.ylim(ylim1,ylim2)
        pl.legend()

        pl.xlabel('$Mass$ $coordinate$', fontsize=20)
        pl.ylabel('$X_{i}$', fontsize=20)
        #pl.title('Mass='+str(tot_mass)+', Time='+str(age)+' years, cycle='+str(mod))
        pl.title('Mass='+str(tot_mass)+', cycle='+str(mod))
        if show:
            pl.show()

def run_regression_1D():
    np.random.seed(42)

    print "create dataset ..."
    N = 50
    rng = np.random.RandomState(42)
    X = np.sort(2 * rng.rand(N, 1) - 1, axis=0)
    Y = np.array([np.pi * np.sin(10 * X).ravel(),
                  np.pi * np.cos(10 * X).ravel()]).T
    Y += (0.5 - rng.rand(*Y.shape))
    Y = Y / np.std(Y, axis=0)

    def plot(model, alpha, fname):
        xx = np.linspace(-1.2, 1.2, 200)[:, None]
        if isinstance(model, IndepSGPR):
            mf, vf = model.predict_f(xx, alpha)
        else:
            # mf, vf = model.predict_f(xx, alpha, use_mean_only=False)
            mf, vf = model.predict_f(xx, alpha, use_mean_only=True)

        colors = ['r', 'b']
        plt.figure()
        for i in range(model.Dout):
            plt.subplot(model.Dout, 1, i + 1)
            plt.plot(X, Y[:, i], 'x', color=colors[i], mew=2)
            zu = model.models[i].zu
            mean_u, var_u = model.models[i].predict_f(zu, alpha)
            plt.plot(xx, mf[:, i], '-', color=colors[i], lw=2)
            plt.fill_between(
                xx[:, 0],
                mf[:, i] - 2 * np.sqrt(vf[:, i]),
                mf[:, i] + 2 * np.sqrt(vf[:, i]),
                color=colors[i], alpha=0.3)
            # plt.errorbar(zu[:, 0], mean_u, yerr=2*np.sqrt(var_u), fmt='ro')
            plt.xlim(-1.2, 1.2)
        plt.savefig(fname)

    # inference
    print "create independent output model and optimize ..."
    M = N
    alpha = 0.01
    indep_model = IndepSGPR(X, Y, M)
    indep_model.train(alpha=alpha)
    plot(indep_model, alpha, '/tmp/reg_indep_multioutput.pdf')

    print "create correlated output model and optimize ..."
    M = N
    ar_model = AutoSGPR(X, Y, M)
    ar_model.train(alpha=alpha)
    plot(ar_model, alpha, '/tmp/reg_autoreg_multioutput.pdf')

def run_step_1D():
    np.random.seed(42)

    def step(x):
        y = x.copy()
        y[y < 0.0] = 0.0
        y[y > 0.0] = 1.0
        return y + 0.02 * np.random.randn(x.shape[0], 1)

    print "create dataset ..."
    N = 100
    X = np.random.rand(N, 1) * 3 - 1.5
    Y = step(X) - 0.5
    # plt.plot(X, Y, 'kx', mew=2)

    def plot(m):
        xx = np.linspace(-3, 3, 100)[:, None]
        mean, var = m.predict_f(xx)
        zu = m.sgp_layers[0].zu
        mean_u, var_u = m.predict_f(zu)
        plt.figure()
        plt.plot(X, Y, 'kx', mew=2)
        plt.plot(xx, mean, 'b', lw=2)
        plt.fill_between(
            xx[:, 0],
            mean[:, 0] - 2 * np.sqrt(var[:, 0]),
            mean[:, 0] + 2 * np.sqrt(var[:, 0]),
            color='blue', alpha=0.2)
        plt.errorbar(zu, mean_u, yerr=2 * np.sqrt(var_u), fmt='ro')

        no_samples = 20
        xx = np.linspace(-3, 3, 500)[:, None]
        f_samples = m.sample_f(xx, no_samples)
        for i in range(no_samples):
            plt.plot(xx, f_samples[:, :, i], linewidth=0.5, alpha=0.5)

        plt.xlim(-3, 3)

    # inference
    print "create model and optimize ..."
    M = 20
    hidden_size = [2]
    model = aep.SDGPR(X, Y, M, hidden_size, lik='Gaussian')
    # model.optimise(method='L-BFGS-B', alpha=1, maxiter=1000)
    model.optimise(method='adam', adam_lr=0.05, alpha=1, maxiter=2000)
    plot(model)
    plt.show()

def plot_model_no_control(model, plot_title='', name_suffix=''):
    # plot function
    mx, vx = model.get_posterior_x()
    mins = np.min(mx, axis=0) - 0.5
    maxs = np.max(mx, axis=0) + 0.5
    nGrid = 50
    xspaced = np.linspace(mins[0], maxs[0], nGrid)
    yspaced = np.linspace(mins[1], maxs[1], nGrid)
    xx, yy = np.meshgrid(xspaced, yspaced)
    Xplot = np.vstack((xx.flatten(), yy.flatten())).T
    mf, vf = model.predict_f(Xplot)
    fig = plt.figure()
    plt.imshow((mf[:, 0]).reshape(*xx.shape),
               vmin=mf.min(), vmax=mf.max(), origin='lower',
               extent=[mins[0], maxs[0], mins[1], maxs[1]], aspect='auto')
    plt.colorbar()
    plt.contour(
        xx, yy, (mf[:, 0]).reshape(*xx.shape),
        colors='k', linewidths=2, zorder=100)
    zu = model.dyn_layer.zu
    plt.plot(zu[:, 0], zu[:, 1], 'wo', mew=0, ms=4)
    for i in range(mx.shape[0] - 1):
        plt.plot(mx[i:i + 2, 0], mx[i:i + 2, 1],
                 '-bo', ms=3, linewidth=2, zorder=101)
    plt.xlabel(r'$x_{t, 1}$')
    plt.ylabel(r'$x_{t, 2}$')
    plt.xlim([mins[0], maxs[0]])
    plt.ylim([mins[1], maxs[1]])
    plt.title(plot_title)
    plt.savefig('/tmp/hh_gpssm_dim_0' + name_suffix + '.pdf')

    fig = plt.figure()
    plt.imshow((mf[:, 1]).reshape(*xx.shape),
               vmin=mf.min(), vmax=mf.max(), origin='lower',
               extent=[mins[0], maxs[0], mins[1], maxs[1]], aspect='auto')
    plt.colorbar()
    plt.contour(
        xx, yy, (mf[:, 1]).reshape(*xx.shape),
        colors='k', linewidths=2, zorder=100)
    zu = model.dyn_layer.zu
    plt.plot(zu[:, 0], zu[:, 1], 'wo', mew=0, ms=4)
    for i in range(mx.shape[0] - 1):
        plt.plot(mx[i:i + 2, 0], mx[i:i + 2, 1],
                 '-bo', ms=3, linewidth=2, zorder=101)
    plt.xlabel(r'$x_{t, 1}$')
    plt.ylabel(r'$x_{t, 2}$')
    plt.xlim([mins[0], maxs[0]])
    plt.ylim([mins[1], maxs[1]])
    plt.title(plot_title)
    plt.savefig('/tmp/hh_gpssm_dim_1' + name_suffix + '.pdf')

def run_regression_1D_pep_inference():
    np.random.seed(42)

    print "create dataset ..."
    N = 200
    X = np.random.rand(N, 1)
    Y = np.sin(12 * X) + 0.5 * np.cos(25 * X) + np.random.randn(N, 1) * 0.2
    # plt.plot(X, Y, 'kx', mew=2)

    def plot(m):
        xx = np.linspace(-0.5, 1.5, 100)[:, None]
        mean, var = m.predict_f(xx)
        zu = m.sgp_layer.zu
        mean_u, var_u = m.predict_f(zu)
        plt.figure()
        plt.plot(X, Y, 'kx', mew=2)
        plt.plot(xx, mean, 'b', lw=2)
        plt.fill_between(
            xx[:, 0],
            mean[:, 0] - 2 * np.sqrt(var[:, 0]),
            mean[:, 0] + 2 * np.sqrt(var[:, 0]),
            color='blue', alpha=0.2)
        plt.errorbar(zu, mean_u, yerr=2 * np.sqrt(var_u), fmt='ro')
        plt.xlim(-0.1, 1.1)

    # inference
    print "create aep model and optimize ..."
    M = 20
    alpha = 0.5
    model_aep = aep.SGPR(X, Y, M, lik='Gaussian')
    model_aep.optimise(method='L-BFGS-B', alpha=alpha, maxiter=2000)
    # plot(model_aep)
    # plt.show()
    # plt.savefig('/tmp/gpr_aep_reg.pdf')

    start_time = time.time()
    model = pep.SGPR_rank_one(X, Y, M, lik='Gaussian')
    model.update_hypers(model_aep.get_hypers())
    # model.update_hypers(model.init_hypers())
    model.run_pep(np.arange(N), 10, alpha=alpha, parallel=False, compute_energy=True)
    end_time = time.time()
    print "sequential updates: %.4f" % (end_time - start_time)
    # plot(model)
    # plt.savefig('/tmp/gpr_pep_reg_seq.pdf')


    start_time = time.time()
    model = pep.SGPR_rank_one(X, Y, M, lik='Gaussian')
    model.update_hypers(model_aep.get_hypers())
    # model.update_hypers(model.init_hypers(Y))
    model.run_pep(np.arange(N), 10, alpha=alpha, parallel=True, compute_energy=False)
    end_time = time.time()
    print "parallel updates: %.4f" % (end_time - start_time)
    # plot(model)
    # plt.savefig('/tmp/gpr_pep_reg_par.pdf')

def run_regression_1D():
    np.random.seed(42)

    print "create dataset ..."
    N = 200
    X = np.random.rand(N, 1)
    Y = np.sin(12 * X) + 0.5 * np.cos(25 * X) + np.random.randn(N, 1) * 0.2
    # plt.plot(X, Y, 'kx', mew=2)

    def plot(m):
        xx = np.linspace(-0.5, 1.5, 100)[:, None]
        mean, var = m.predict_f(xx)
        zu = m.sgp_layer.zu
        mean_u, var_u = m.predict_f(zu)
        plt.figure()
        plt.plot(X, Y, 'kx', mew=2)
        plt.plot(xx, mean, 'b', lw=2)
        plt.fill_between(
            xx[:, 0],
            mean[:, 0] - 2 * np.sqrt(var[:, 0]),
            mean[:, 0] + 2 * np.sqrt(var[:, 0]),
            color='blue', alpha=0.2)
        plt.errorbar(zu, mean_u, yerr=2 * np.sqrt(var_u), fmt='ro')
        plt.xlim(-0.1, 1.1)

    # inference
    print "create model and optimize ..."
    M = 20
    alpha = 0.1
    model_aep = aep.SGPR(X, Y, M, lik='Gaussian')
    model_aep.optimise(method='L-BFGS-B', alpha=alpha, maxiter=2000)
    plot(model_aep)
    plt.savefig('/tmp/gpr_aep_reg.pdf')

    start_time = time.time()
    model = pep.SGPR(X, Y, M, lik='Gaussian')
    model.update_hypers(model_aep.get_hypers())
    # model.update_hypers(model.init_hypers())
    model.inference(alpha=alpha, no_epochs=10)
    end_time = time.time()
    print "sequential updates: %.4f" % (end_time - start_time)
    plot(model)
    plt.savefig('/tmp/gpr_pep_reg_seq.pdf')


    start_time = time.time()
    model = pep.SGPR(X, Y, M, lik='Gaussian')
    model.update_hypers(model_aep.get_hypers())
    # model.update_hypers(model.init_hypers())
    model.inference(alpha=alpha, no_epochs=10, parallel=True)
    end_time = time.time()
    print "parallel updates: %.4f" % (end_time - start_time)
    plot(model)
    plt.savefig('/tmp/gpr_pep_reg_par.pdf')

    # plt.show()