Python matplotlib.pylab 模块，subplots() 实例源码

def plot_trajectory(name):
    STEPS = 600
    DELTA = 1 if name != 'linear' else 0.1
    trajectory = create_trajectory(name, STEPS)

    x = [trajectory.get_position_at(i * DELTA).x for i in range(STEPS)]
    y = [trajectory.get_position_at(i * DELTA).y for i in range(STEPS)]

    trajectory_fig, trajectory_plot = plt.subplots(1, 1)
    trajectory_plot.plot(x, y, label='trajectory', lw=3)
    trajectory_plot.set_title(name.title() + ' Trajectory', fontsize=20)
    trajectory_plot.set_xlabel(r'$x{\rm[m]}$', fontsize=18)
    trajectory_plot.set_ylabel(r'$y{\rm[m]}$', fontsize=18)
    trajectory_plot.legend(loc=0)
    trajectory_plot.grid()
    plt.show()

def grid_search_gamma(rbf_svm, X, y):
    ## grid search - gamma only
    # use a full grid over all parameters
    param_grid = {'gamma': np.logspace(-15, 4, num = 5000, base = 2.0)}
    grid_search = GridSearchCV(rbf_svm, param_grid = param_grid, scoring = 'roc_auc',
                               cv = 10, pre_dispatch = '2*n_jobs', n_jobs = -1)
    # re-fit on the whole training data
    grid_search.fit(X, y)
    grid_search_scores = [score[1] for score in grid_search.grid_scores_]
    print('Best parameters : {}'.format(grid_search.best_params_))
    print('Best score : {}'.format(grid_search.best_score_))

    # set canvas
    fig, ax = plt.subplots(1, 1)
    # ax.scatter(X[:, 0], X[:, 1], c = y)
    ax.plot(param_grid['gamma'], grid_search_scores)
    ax.set_title('AUC = f(gamma, C = 1.0)', fontsize = 'large')
    ax.set_xlabel('gamma', fontsize = 'medium')
    ax.set_ylabel('AUC', fontsize = 'medium')

    return fig

def plotValueFunction(self, valueFunction, prefix):
        '''3d plot of a value function.'''
        fig, ax = plt.subplots(subplot_kw = dict(projection = '3d'))
        X, Y = np.meshgrid(np.arange(self.numCols), np.arange(self.numRows))
        Z = valueFunction.reshape(self.numRows, self.numCols)

        for i in xrange(len(X)):
            for j in xrange(len(X[i])/2):
                tmp = X[i][j]
                X[i][j] = X[i][len(X[i]) - j - 1]
                X[i][len(X[i]) - j - 1] = tmp

        my_col = cm.jet(np.random.rand(Z.shape[0],Z.shape[1]))

        ax.plot_surface(X, Y, Z, rstride = 1, cstride = 1,
            cmap = plt.get_cmap('jet'))
        plt.gca().view_init(elev=30, azim=30)
        plt.savefig(self.outputPath + prefix + 'value_function.png')
        plt.close()

def plot_classes(y, cord, names, test_error, message=""):
    plt.close("all")
    cord = np.array(cord)
    colors = ('b', 'g', 'r', 'c', 'm', 'y', 'k')
    un = np.unique(y)
    fig, ax = plt.subplots()
    for u, col in zip(un, colors):
        ind = np.argwhere(y == u)
        x = cord[ind, :]
        x = x.reshape(x.shape[0], cord.shape[1])
        ax.scatter(x[:, 0], x[:, 1], label="class:" + str(u),
                   color=col)

    plt.legend(loc='upper right', fancybox=True, shadow=True, prop={'size': 8})
    fig.suptitle(
        "Output prediction. Test error:" + str(test_error*100) + "%. " +
        message, fontsize=8)
    return fig

def histogram(data,variables):

    sns.set_context("notebook", font_scale=1.5, rc={"lines.linewidth": 0})

    sns.set_style('white')

    var_length = len(variables)

    fig, axes = plt.subplots(1, var_length, figsize=(19, 5))

    for i in range(var_length):

        axes[i].hist(data[variables[i]],lw=0,color="indianred",bins=8);
        axes[i].tick_params(axis='both', which='major', pad=15)
        axes[i].set_xlabel(variables[i])
        axes[i].set_yticklabels("");

    sns.despine(left=True)

def show_bars_over_time(
        task_output_path=None,
        query_laps=[0, 1, 2, 5, None],
        ncols=10):
    '''
    '''
    nrows = len(query_laps)
    fig_handle, ax_handles_RC = pylab.subplots(
        figsize=(SMALL_FIG_SIZE[0] * ncols, SMALL_FIG_SIZE[1] * nrows),
        nrows=nrows, ncols=ncols, sharex=True, sharey=True)
    for row_id, lap_val in enumerate(query_laps):
        cur_model, lap_val = bnpy.load_model_at_lap(task_output_path, lap_val)
        cur_topics_KV = cur_model.obsModel.getTopics()
        # Plot the current model
        cur_ax_list = ax_handles_RC[row_id].flatten().tolist()
        bnpy.viz.BarsViz.show_square_images(
            cur_topics_KV,
            vmin=0.0, vmax=0.06,
            ax_list=cur_ax_list)
        cur_ax_list[0].set_ylabel("lap: %d" % lap_val)
    pylab.tight_layout()

###############################################################################
#
# Show the clusters over time

def show_bars_over_time(
        task_output_path=None,
        query_laps=[0, 1, 2, 5, None],
        ncols=10):
    '''
    '''
    nrows = len(query_laps)
    fig_handle, ax_handles_RC = pylab.subplots(
        figsize=(SMALL_FIG_SIZE[0] * ncols, SMALL_FIG_SIZE[1] * nrows),
        nrows=nrows, ncols=ncols, sharex=True, sharey=True)
    for row_id, lap_val in enumerate(query_laps):
        cur_model, lap_val = bnpy.load_model_at_lap(task_output_path, lap_val)
        cur_topics_KV = cur_model.obsModel.getTopics()
        # Plot the current model
        cur_ax_list = ax_handles_RC[row_id].flatten().tolist()
        bnpy.viz.BarsViz.show_square_images(
            cur_topics_KV,
            vmin=0.0, vmax=0.06,
            ax_list=cur_ax_list)
        cur_ax_list[0].set_ylabel("lap: %d" % lap_val)
    pylab.tight_layout()

###############################################################################
#
# Show the clusters over time

def plot_true_and_augmented_data(sample,noised_sample,label,n_examples):
    output_dir = os.path.split(FLAGS.output)[0]
    # Save augmented data
    plt.clf()
    fig, ax = plt.subplots(3,1)
    for t in range(noised_sample.shape[1]):
        ax[t].plot(noised_sample[:,t])
        ax[t].set_xlabel('time (samples)')
        ax[t].set_ylabel('amplitude')
    ax[0].set_title('window {:03d}, cluster_id: {}'.format(n_examples,label))
    plt.savefig(os.path.join(output_dir, "augmented_data",
                            'augmented_{:03d}.pdf'.format(n_examples)))
    plt.close()

    # Save true data
    plt.clf()
    fig, ax = plt.subplots(3,1)
    for t in range(sample.shape[1]):
        ax[t].plot(sample[:,t])
        ax[t].set_xlabel('time (samples)')
        ax[t].set_ylabel('amplitude')
    ax[0].set_title('window {:03d}, cluster_id: {}'.format(n_examples,label))
    plt.savefig(os.path.join(output_dir, "true_data",
                            'true__{:03d}.pdf'.format(n_examples)))
    plt.close()

def plot(ts):
        if not plt:
            print ""
        fig, ax = plt.subplots()
        lined = dict()

        ax.set_title('Click on legend line to toggle line on/off')
        lines = [ax.plot(ts[col], label=col) for col in ts.columns]
        leg = ax.legend(loc='best')

        for legline, origline in zip(leg.get_lines(), lines):
            legline.set_picker(5)  # 5 pts tolerance
            lined[legline] = origline[0]

        def onpick(event):
            # on the pick event, find the orig line corresponding to the
            # legend proxy line, and toggle the visibility
            legline = event.artist
            origline = lined[legline]
            vis = not origline.get_visible()
            origline.set_visible(vis)
            # Change the alpha on the line in the legend so we can see what lines
            # have been toggled
            if vis:
                legline.set_alpha(1.0)
            else:
                legline.set_alpha(0.2)
            fig.canvas.draw()

        fig.canvas.mpl_connect('pick_event', onpick)
        plt.show(False)

def imageSegmentor(imageFilePath, matFilePath):      


    mat = readMatFile(matFilePath);                                             # read mat file                                 
    image = getImage(imageFilePath);                                            # input the image
    typeOfFruit = getTypeOfFruit(image);                                        # on basis of counting or temperature value there are 2 types of fruit

    plt.imshow(image);

    _fft = getFFT(image);
    _mag = getMag(_fft);
    _ang = getAngleInDegrees(_fft);    

    edges = edgeDetector(image);                                                # detects the edges of the image
    _segmentation = segmentation(image, typeOfFruit);                           # segments different parts of image    
    filteredImage = filterImageFromSegmentation(image, _segmentation);          # filter the object part of image

    outputMatrix = imageMapping(filteredImage, mat['IR']);                      # map the value part of the image and else 0

    prefix =  re.split('IR_|.pgm', imageFilePath)[0];
    postfix = re.split('IR_|.pgm', imageFilePath)[1];    
    nameOfFile = prefix + "csv_" 
    nameOfFile = nameOfFile + postfix;
    print(nameOfFile);    
    writeToCSV(outputMatrix, nameOfFile);                                      # write it to the CSV file
    writeFF2CSV(outputMatrix, _mag, _ang, nameOfFile);  

    fig, ((fig1, fig2), (fig3, fig4)) = plt.subplots(2, 2, figsize = (10, 8));  # subplot the different plots
    fig1.imshow(image, cmap = plt.cm.gray);                                     # colormap used here is gray    
    fig2.imshow(image, cmap = plt.cm.gray); 
    fig3.imshow(edges, cmap = plt.cm.gray);
    fig4.imshow(filteredImage, cmap = plt.cm.gray);

    return

# header file

def save(self, out_path):
        '''Saves a figure for the monitor

        Args:
            out_path: str
        '''

        plt.clf()
        np.set_printoptions(precision=4)
        font = {
            'size': 7
        }
        matplotlib.rc('font', **font)
        y = 2
        x = ((len(self.d) - 1) // y) + 1
        fig, axes = plt.subplots(y, x)
        fig.set_size_inches(20, 8)

        for j, (k, v) in enumerate(self.d.iteritems()):
            ax = axes[j // x, j % x]
            ax.plot(v, label=k)
            if k in self.d_valid.keys():
                ax.plot(self.d_valid[k], label=k + '(valid)')
            ax.set_title(k)
            ax.legend()

        plt.tight_layout()
        plt.savefig(out_path, facecolor=(1, 1, 1))
        plt.close()

def generate_box_plot(dataset, methods, position_rmses, orientation_rmses):

  num_methods = len(methods)
  x_ticks = np.linspace(0., 1., num_methods)

  width = 0.3 / num_methods
  spacing = 0.3 / num_methods
  fig, ax1 = plt.subplots()
  ax1.set_ylabel('RMSE position [m]', color='b')
  ax1.tick_params('y', colors='b')
  fig.suptitle(
      "Hand-Eye Calibration Method Error {}".format(dataset), fontsize='24')
  bp_position = ax1.boxplot(position_rmses, 0, '',
                            positions=x_ticks - spacing, widths=width)
  plt.setp(bp_position['boxes'], color='blue', linewidth=line_width)
  plt.setp(bp_position['whiskers'], color='blue', linewidth=line_width)
  plt.setp(bp_position['fliers'], color='blue',
           marker='+', linewidth=line_width)
  plt.setp(bp_position['caps'], color='blue', linewidth=line_width)
  plt.setp(bp_position['medians'], color='blue', linewidth=line_width)
  ax2 = ax1.twinx()
  ax2.set_ylabel('RMSE Orientation [$^\circ$]', color='g')
  ax2.tick_params('y', colors='g')
  bp_orientation = ax2.boxplot(
      orientation_rmses, 0, '', positions=x_ticks + spacing, widths=width)
  plt.setp(bp_orientation['boxes'], color='green', linewidth=line_width)
  plt.setp(bp_orientation['whiskers'], color='green', linewidth=line_width)
  plt.setp(bp_orientation['fliers'], color='green',
           marker='+')
  plt.setp(bp_orientation['caps'], color='green', linewidth=line_width)
  plt.setp(bp_orientation['medians'], color='green', linewidth=line_width)

  plt.xticks(x_ticks, methods)
  plt.xlim(x_ticks[0] - 2.5 * spacing, x_ticks[-1] + 2.5 * spacing)

  plt.show()

def generate_time_plot(methods, datasets, runtimes_per_method, colors):
  num_methods = len(methods)
  num_datasets = len(datasets)
  x_ticks = np.linspace(0., 1., num_methods)

  width = 0.6 / num_methods / num_datasets
  spacing = 0.4 / num_methods / num_datasets
  fig, ax1 = plt.subplots()
  ax1.set_ylabel('Time [s]', color='b')
  ax1.tick_params('y', colors='b')
  ax1.set_yscale('log')
  fig.suptitle("Hand-Eye Calibration Method Timings", fontsize='24')
  handles = []
  for i, dataset in enumerate(datasets):
    runtimes = [runtimes_per_method[dataset][method] for method in methods]
    bp = ax1.boxplot(
        runtimes, 0, '',
        positions=(x_ticks + (i - num_datasets / 2. + 0.5) *
                   spacing * 2),
        widths=width)
    plt.setp(bp['boxes'], color=colors[i], linewidth=line_width)
    plt.setp(bp['whiskers'], color=colors[i], linewidth=line_width)
    plt.setp(bp['fliers'], color=colors[i],
             marker='+', linewidth=line_width)
    plt.setp(bp['medians'], color=colors[i],
             marker='+', linewidth=line_width)
    plt.setp(bp['caps'], color=colors[i], linewidth=line_width)
    handles.append(mpatches.Patch(color=colors[i], label=dataset))
  plt.legend(handles=handles, loc=2)

  plt.xticks(x_ticks, methods)
  plt.xlim(x_ticks[0] - 2.5 * spacing * num_datasets,
           x_ticks[-1] + 2.5 * spacing * num_datasets)

  plt.show()

def plot_distortion(training_data_instances):
    # dimension of a training data instance
    d = training_data_instances.shape[1]
    # first m instances considered
    m = 20

    fig, axes = plt.subplots(1, 1)
    fig.suptitle("Distortion of random projection", fontsize = "x-large")

    for k in [50, 100, 500]:
        ## generate random projection matrix
        random_projection_matrix =  generate_random_projection_matrix(k, d)
        ## random projection
        m_instances = training_data_instances[0:m]
        projected_m_instances = np.dot(m_instances, np.transpose(random_projection_matrix))
        # print random_projected_matrix[0], random_projected_matrix.shape
        ## evaluate distortion - line chart
        m_instances_distortions = []
        for i in range(m):
            for j in range(i + 1, m):
                m_instances_distortions.append(euclidean(projected_m_instances[i], projected_m_instances[j]) / euclidean(m_instances[i], m_instances[j]))
        m_instances_distortions = np.array(m_instances_distortions)
        mean, std = np.mean(m_instances_distortions), np.std(m_instances_distortions)
        # line chart
        axes.plot(m_instances_distortions, label = "k=" + str(k))
        axes.plot([0, m_instances_distortions.size], [mean, mean], label = "k=" + str(k) + ", mean = " + str(round(mean, 4)))

        print "k = ", k, "distortion =", mean, "+-", std
    axes.set_xlabel("pairs of instances", fontsize = "large")
    axes.set_ylabel("distortion", fontsize = "large")
    axes.legend(loc = "center right", fontsize = "medium")
    plt.show()

def plot_gplvm_vfe_gaussian_stochastic():
    N_train = 2000
    M = 50
    D = 2
    Q = 3
    y_train = np.random.randn(N_train, Q)
    model = vfe.SGPLVM(y_train, D, M, lik='Gaussian')

    # init hypers, inducing points and q(u) params
    params = model.init_hypers(y_train)
    logZ, grad_all = model.objective_function(params, N_train)
    mbs = np.logspace(-2, 0, 10)
    reps = 20
    times = np.zeros(len(mbs))
    objs = np.zeros((len(mbs), reps))
    for i, mb in enumerate(mbs):
        no_points = int(N_train * mb)
        start_time = time.time()
        for k in range(reps):
            objs[i, k] = model.objective_function(
                params, no_points, alpha=alpha)[0]
        times[i] = time.time() - start_time

    f, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
    ax1.plot(mbs, times, 'x-')
    ax1.set_xlabel("Minibatch proportion")
    ax1.set_ylabel("Time taken")
    ax1.set_xscale("log", nonposx='clip')

    ax2.plot(mbs, objs, 'kx')
    ax2.axhline(logZ, color='b')
    ax2.set_xlabel("Minibatch proportion")
    ax2.set_ylabel("ELBO estimates")
    ax2.set_xscale("log", nonposx='clip')
    plt.savefig('/tmp/gaussian_stochastic_vfe_gplvm.pdf')

def plot_gplvm_vfe_probit_stochastic():
    N_train = 2000
    M = 50
    D = 2
    Q = 3
    y_train = 2 * np.random.randint(0, 2, size=(N_train, Q)) - 1
    model = vfe.SGPLVM(y_train, D, M, lik='Probit')

    # init hypers, inducing points and q(u) params
    params = model.init_hypers(y_train)
    logZ, grad_all = model.objective_function(params, N_train)
    mbs = np.logspace(-2, 0, 10)
    reps = 20
    times = np.zeros(len(mbs))
    objs = np.zeros((len(mbs), reps))
    for i, mb in enumerate(mbs):
        no_points = int(N_train * mb)
        start_time = time.time()
        for k in range(reps):
            objs[i, k] = model.objective_function(
                params, no_points, alpha=alpha)[0]
        times[i] = time.time() - start_time

    f, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
    ax1.plot(mbs, times, 'x-')
    ax1.set_xlabel("Minibatch proportion")
    ax1.set_ylabel("Time taken")
    ax1.set_xscale("log", nonposx='clip')

    ax2.plot(mbs, objs, 'kx')
    ax2.axhline(logZ, color='b')
    ax2.set_xlabel("Minibatch proportion")
    ax2.set_ylabel("ELBO estimates")
    ax2.set_xscale("log", nonposx='clip')
    plt.savefig('/tmp/probit_stochastic_vfe_gplvm.pdf')

def plot_gpr_vfe_gaussian_stochastic():
    N_train = 2000
    M = 50
    D = 2
    Q = 3
    y_train = np.random.randn(N_train, Q)
    x_train = np.random.randn(N_train, D)
    model = vfe.SGPR(x_train, y_train, M, lik='Gaussian')

    # init hypers, inducing points and q(u) params
    params = model.init_hypers(y_train)
    logZ, grad_all = model.objective_function(params, N_train)
    mbs = np.logspace(-2, 0, 10)
    reps = 20
    times = np.zeros(len(mbs))
    objs = np.zeros((len(mbs), reps))
    for i, mb in enumerate(mbs):
        no_points = int(N_train * mb)
        start_time = time.time()
        for k in range(reps):
            objs[i, k] = model.objective_function(
                params, no_points, alpha=alpha)[0]
        times[i] = time.time() - start_time

    f, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
    ax1.plot(mbs, times, 'x-')
    ax1.set_xlabel("Minibatch proportion")
    ax1.set_ylabel("Time taken")
    ax1.set_xscale("log", nonposx='clip')

    ax2.plot(mbs, objs, 'kx')
    ax2.axhline(logZ, color='b')
    ax2.set_xlabel("Minibatch proportion")
    ax2.set_ylabel("ELBO estimates")
    ax2.set_xscale("log", nonposx='clip')
    plt.savefig('/tmp/gaussian_stochastic_vfe_gpr.pdf')

def plot_gpr_vfe_probit_stochastic():
    N_train = 2000
    alpha = 0.5
    M = 50
    D = 2
    Q = 3
    x_train = np.random.randn(N_train, D)
    y_train = 2 * np.random.randint(0, 2, size=(N_train, Q)) - 1
    model = vfe.SGPR(x_train, y_train, M, lik='Probit')

    # init hypers, inducing points and q(u) params
    params = model.init_hypers(y_train)
    logZ, grad_all = model.objective_function(params, N_train)
    mbs = np.logspace(-2, 0, 10)
    reps = 20
    times = np.zeros(len(mbs))
    objs = np.zeros((len(mbs), reps))
    for i, mb in enumerate(mbs):
        no_points = int(N_train * mb)
        start_time = time.time()
        for k in range(reps):
            objs[i, k] = model.objective_function(
                params, no_points, alpha=alpha)[0]
        times[i] = time.time() - start_time

    f, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
    ax1.plot(mbs, times, 'x-')
    ax1.set_xlabel("Minibatch proportion")
    ax1.set_ylabel("Time taken")
    ax1.set_xscale("log", nonposx='clip')

    ax2.plot(mbs, objs, 'kx')
    ax2.axhline(logZ, color='b')
    ax2.set_xlabel("Minibatch proportion")
    ax2.set_ylabel("ELBO estimates")
    ax2.set_xscale("log", nonposx='clip')
    plt.savefig('/tmp/probit_stochastic_vfe_gpr.pdf')

def plot_dgpr_vfe_gaussian_stochastic():
    N_train = 2000
    M = 50
    D = 2
    Q = 3
    y_train = np.random.randn(N_train, Q)
    x_train = np.random.randn(N_train, D)
    hidden_size = [3, 2]
    model = vfe.SDGPR(x_train, y_train, M, hidden_size, lik='Gaussian')

    # init hypers, inducing points and q(u) params
    params = model.init_hypers(y_train)
    logZ, grad_all = model.objective_function(params, N_train)
    mbs = np.logspace(-2, 0, 10)
    reps = 20
    times = np.zeros(len(mbs))
    objs = np.zeros((len(mbs), reps))
    for i, mb in enumerate(mbs):
        no_points = int(N_train * mb)
        start_time = time.time()
        for k in range(reps):
            objs[i, k] = model.objective_function(
                params, no_points, alpha=1.0)[0]
        times[i] = time.time() - start_time

    f, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
    ax1.plot(mbs, times, 'x-')
    ax1.set_xlabel("Minibatch proportion")
    ax1.set_ylabel("Time taken")
    ax1.set_xscale("log", nonposx='clip')

    ax2.plot(mbs, objs, 'kx')
    ax2.axhline(logZ, color='b')
    ax2.set_xlabel("Minibatch proportion")
    ax2.set_ylabel("ELBO estimates")
    ax2.set_xscale("log", nonposx='clip')
    plt.savefig('/tmp/gaussian_stochastic_vfe_dgpr.pdf')

def plot_gplvm_aep_gaussian_stochastic():
    N_train = 2000
    alpha = 0.5
    M = 50
    D = 2
    Q = 3
    y_train = np.random.randn(N_train, Q)
    model = aep.SGPLVM(y_train, D, M, lik='Gaussian')

    # init hypers, inducing points and q(u) params
    params = model.init_hypers(y_train)
    logZ, grad_all = model.objective_function(params, N_train, alpha=alpha)
    mbs = np.logspace(-2, 0, 10)
    reps = 20
    times = np.zeros(len(mbs))
    objs = np.zeros((len(mbs), reps))
    for i, mb in enumerate(mbs):
        no_points = int(N_train * mb)
        start_time = time.time()
        for k in range(reps):
            objs[i, k] = model.objective_function(
                params, no_points, alpha=alpha)[0]
        times[i] = time.time() - start_time

    f, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
    ax1.plot(mbs, times, 'x-')
    ax1.set_xlabel("Minibatch proportion")
    ax1.set_ylabel("Time taken")
    ax1.set_xscale("log", nonposx='clip')

    ax2.plot(mbs, objs, 'kx')
    ax2.axhline(logZ, color='b')
    ax2.set_xlabel("Minibatch proportion")
    ax2.set_ylabel("ELBO estimates")
    ax2.set_xscale("log", nonposx='clip')
    plt.savefig('/tmp/gaussian_stochastic_aep_gplvm.pdf')

def plot_gplvm_aep_probit_stochastic():
    N_train = 2000
    alpha = 0.5
    M = 50
    D = 2
    Q = 3
    y_train = 2 * np.random.randint(0, 2, size=(N_train, Q)) - 1
    model = aep.SGPLVM(y_train, D, M, lik='Probit')

    # init hypers, inducing points and q(u) params
    params = model.init_hypers(y_train)
    logZ, grad_all = model.objective_function(params, N_train, alpha=alpha)
    mbs = np.logspace(-2, 0, 10)
    reps = 20
    times = np.zeros(len(mbs))
    objs = np.zeros((len(mbs), reps))
    for i, mb in enumerate(mbs):
        no_points = int(N_train * mb)
        start_time = time.time()
        for k in range(reps):
            objs[i, k] = model.objective_function(
                params, no_points, alpha=alpha)[0]
        times[i] = time.time() - start_time

    f, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
    ax1.plot(mbs, times, 'x-')
    ax1.set_xlabel("Minibatch proportion")
    ax1.set_ylabel("Time taken")
    ax1.set_xscale("log", nonposx='clip')

    ax2.plot(mbs, objs, 'kx')
    ax2.axhline(logZ, color='b')
    ax2.set_xlabel("Minibatch proportion")
    ax2.set_ylabel("ELBO estimates")
    ax2.set_xscale("log", nonposx='clip')
    plt.savefig('/tmp/probit_stochastic_aep_gplvm.pdf')

def plot_gpr_aep_gaussian_stochastic():
    N_train = 2000
    alpha = 0.5
    M = 50
    D = 2
    Q = 3
    y_train = np.random.randn(N_train, Q)
    x_train = np.random.randn(N_train, D)
    model = aep.SGPR(x_train, y_train, M, lik='Gaussian')

    # init hypers, inducing points and q(u) params
    params = model.init_hypers(y_train)
    logZ, grad_all = model.objective_function(params, N_train, alpha=alpha)
    mbs = np.logspace(-2, 0, 10)
    reps = 20
    times = np.zeros(len(mbs))
    objs = np.zeros((len(mbs), reps))
    for i, mb in enumerate(mbs):
        no_points = int(N_train * mb)
        start_time = time.time()
        for k in range(reps):
            objs[i, k] = model.objective_function(
                params, no_points, alpha=alpha)[0]
        times[i] = time.time() - start_time

    f, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
    ax1.plot(mbs, times, 'x-')
    ax1.set_xlabel("Minibatch proportion")
    ax1.set_ylabel("Time taken")
    ax1.set_xscale("log", nonposx='clip')

    ax2.plot(mbs, objs, 'kx')
    ax2.axhline(logZ, color='b')
    ax2.set_xlabel("Minibatch proportion")
    ax2.set_ylabel("ELBO estimates")
    ax2.set_xscale("log", nonposx='clip')
    plt.savefig('/tmp/gaussian_stochastic_aep_gpr.pdf')

def plot_gpr_aep_probit_stochastic():
    N_train = 2000
    alpha = 0.5
    M = 50
    D = 2
    Q = 3
    x_train = np.random.randn(N_train, D)
    y_train = 2 * np.random.randint(0, 2, size=(N_train, Q)) - 1
    model = aep.SGPR(x_train, y_train, M, lik='Probit')

    # init hypers, inducing points and q(u) params
    params = model.init_hypers(y_train)
    logZ, grad_all = model.objective_function(params, N_train, alpha=alpha)
    mbs = np.logspace(-2, 0, 10)
    reps = 20
    times = np.zeros(len(mbs))
    objs = np.zeros((len(mbs), reps))
    for i, mb in enumerate(mbs):
        no_points = int(N_train * mb)
        start_time = time.time()
        for k in range(reps):
            objs[i, k] = model.objective_function(
                params, no_points, alpha=alpha)[0]
        times[i] = time.time() - start_time

    f, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
    ax1.plot(mbs, times, 'x-')
    ax1.set_xlabel("Minibatch proportion")
    ax1.set_ylabel("Time taken")
    ax1.set_xscale("log", nonposx='clip')

    ax2.plot(mbs, objs, 'kx')
    ax2.axhline(logZ, color='b')
    ax2.set_xlabel("Minibatch proportion")
    ax2.set_ylabel("ELBO estimates")
    ax2.set_xscale("log", nonposx='clip')
    plt.savefig('/tmp/probit_stochastic_aep_gpr.pdf')

def plot_dgpr_aep_gaussian_stochastic():
    N_train = 2000
    M = 50
    D = 2
    Q = 3
    y_train = np.random.randn(N_train, Q)
    x_train = np.random.randn(N_train, D)
    hidden_size = [3, 2]
    model = aep.SDGPR(x_train, y_train, M, hidden_size, lik='Gaussian')

    # init hypers, inducing points and q(u) params
    params = model.init_hypers(y_train)
    logZ, grad_all = model.objective_function(params, N_train, alpha=1.0)
    mbs = np.logspace(-2, 0, 10)
    reps = 20
    times = np.zeros(len(mbs))
    objs = np.zeros((len(mbs), reps))
    for i, mb in enumerate(mbs):
        no_points = int(N_train * mb)
        start_time = time.time()
        for k in range(reps):
            objs[i, k] = model.objective_function(
                params, no_points, alpha=1.0)[0]
        times[i] = time.time() - start_time

    f, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
    ax1.plot(mbs, times, 'x-')
    ax1.set_xlabel("Minibatch proportion")
    ax1.set_ylabel("Time taken")
    ax1.set_xscale("log", nonposx='clip')

    ax2.plot(mbs, objs, 'kx')
    ax2.axhline(logZ, color='b')
    ax2.set_xlabel("Minibatch proportion")
    ax2.set_ylabel("ELBO estimates")
    ax2.set_xscale("log", nonposx='clip')
    plt.savefig('/tmp/gaussian_stochastic_aep_dgpr.pdf')

def main(args_dict):
    # Extract configuration from command line arguments
    MK = np.array(args_dict['MK'])
    M = 100
    K = MK / M
    print('M = %d; K = %d' % (M, K))
    x_type = args_dict['x_type']
    deltas = args_dict['deltas']
    do_confidence = args_dict['confidence']

    # Load data from JSON files generated by (non-public) Matlab code
    jsons = [json_load('data/bandits_normal_delta%s_MK%d.json' % (delta, MK)) for delta in deltas]
    lnZs = np.array([json['lnZ'] for json in jsons])
    MAPs = np.array([json['MAPs_ttest'] for json in jsons])

    # Estimate estimator MSEs for the various tricks (as specified by alphas)
    alphas = np.linspace(-0.2, 1.5, 100)
    MSEs, MSEs_stdev = MAPs_to_estimator_MSE_vs_alpha(1, MAPs, lnZs, alphas, K)

    # Set up plot
    matplotlib_configure_as_notebook()
    fig, ax = plt.subplots(1, 1, facecolor='w', figsize=(4.25, 3.25))
    ax.set_xlabel('trick parameter $\\alpha$')
    ax.set_ylabel('MSE of estimator of $\ln Z$')

    # Plot the MSEs
    labels = ['$\\delta = %g$' % (delta) for delta in deltas]
    colors = [plt.cm.plasma((np.log10(delta) - (-3)) / (0 - (-3))) for delta in deltas]
    plot_MSEs_to_axis(ax, alphas, MSEs, MSEs_stdev, do_confidence, labels, colors)

    # Finalize plot
    for vertical in [0.0, 1.0]:
        ax.axvline(vertical, color='black', linestyle='dashed', alpha=.7)
    ax.annotate('Gumbel trick', xy=(0.0, 0.0052), rotation=90, horizontalalignment='right', verticalalignment='bottom')
    ax.annotate('Exponential trick', xy=(1.0, 0.0052), rotation=90, horizontalalignment='right', verticalalignment='bottom')
    lgd = ax.legend(loc='upper center')
    ax.set_ylim((5*1e-3, 5*1e-2))
    save_plot(fig, 'figures/fig3b', bbox_extra_artists=(lgd,))

def main(args_dict):
    # Set up plot
    matplotlib_configure_as_notebook()
    fig, ax = plt.subplots(1, 2, facecolor='w', figsize=(9.25, 3.25))

    # Estimating Z
    Ms = np.arange(3, args_dict['M']+1)

    ax[0].set_xlabel('number of samples $M$')
    ax[0].set_ylabel('MSE of $\hat{Z}$, in units of $Z^2$')
    ax[0].set_xlim((np.min(Ms), np.max(Ms)))
    ax[0].set_xscale('log')
    ax[0].set_yscale('log')
    ax[0].grid(b=True, which='major', linestyle='dotted', lw=.5, color='black', alpha=0.5)

    ax[0].plot(Ms, Z_Gumbel_MSE(Ms), linestyle='-', color=tableau20(0), label='Gumbel: MSE')
    ax[0].plot(Ms, Z_Gumbel_var(Ms), linestyle='dashed', color=tableau20(0), label='Gumbel: var')
    ax[0].plot(Ms, Z_Exponential_MSE(Ms), linestyle='-', color=tableau20(2), label='Exponential: MSE')
    ax[0].plot(Ms, Z_Exponential_var(Ms), linestyle='dashed', color=tableau20(2), label='Exponential: var')

    # Estimating ln Z
    Ms = np.arange(1, args_dict['M']+1)

    ax[1].set_xlabel('number of samples $M$')
    ax[1].set_ylabel('MSE of $\widehat{\ln Z}$, in units of $1$')
    ax[1].set_xlim((np.min(Ms), np.max(Ms)))
    ax[1].set_xscale('log')
    ax[1].set_yscale('log')
    ax[1].grid(b=True, which='major', linestyle='dotted', lw=0.5, color='black', alpha=0.5)

    ax[1].plot(Ms, lnZ_Gumbel_MSE(Ms), linestyle='-', color=tableau20(0), label='Gumbel: MSE')
    ax[1].plot(Ms, lnZ_Exponential_MSE(Ms), linestyle='-', color=tableau20(2), label='Exponential: MSE')
    ax[1].plot(Ms, lnZ_Exponential_var(Ms), linestyle='dashed', color=tableau20(2), label='Exponential: var')

    # Finalize plot
    lgd0 = ax[0].legend()
    lgd1 = ax[1].legend()
    plt.tight_layout()
    save_plot(fig, 'figures/fig1', (lgd0, lgd1,))

def plotBasisFunctions(self, eigenvalues, eigenvectors):
        '''3d plot of the basis function. Right now I am plotting eigenvectors,
           so each coordinate of the eigenvector correspond to the value to be
           plotted for the correspondent state.''' 
        for i in xrange(len(eigenvalues)):  
            fig, ax = plt.subplots(subplot_kw = dict(projection = '3d'))
            X, Y = np.meshgrid(np.arange(self.numRows), np.arange(self.numCols))
            Z = eigenvectors[:,i].reshape(self.numCols, self.numRows)

            for ii in xrange(len(X)):
                for j in xrange(len(X[ii])/2):
                    tmp = X[ii][j]
                    X[ii][j] = X[ii][len(X[ii]) - j - 1]
                    X[ii][len(X[ii]) - j - 1] = tmp

            my_col = cm.jet(np.random.rand(Z.shape[0],Z.shape[1]))

            ax.plot_surface(X, Y, Z, rstride = 1, cstride = 1,
                cmap = plt.get_cmap('jet'))
            plt.gca().view_init(elev=30, azim=30)
            plt.savefig(self.outputPath + str(i) + '_eig' + '.png')
            plt.close()


        plt.plot(eigenvalues, 'o')
        plt.savefig(self.outputPath + 'eigenvalues.png')

def plot_debug_grad(debug, tag, fold_exp, trg):
    plt.close("all")
    # f = plt.figure(figsize=(15, 10.8), dpi=300)
    nbr_rows = int(len(debug["grad_sup"][0])/2)
    f, axs = plt.subplots(nbr_rows, 2, sharex=True, sharey=False,
                          figsize=(15, 12.8), dpi=300)

    if trg == "sup":
        grad = np.array(debug["grad_sup"])
    elif trg == "hint":
        grad = np.array(debug["grad_hint"])
    print grad.shape, trg
    j = 0
    for i in range(0, nbr_rows*2, 2):
        w_vl = grad[:, i]
        b_vl = grad[:, i+1]
        axs[j, 0].plot(w_vl, label=trg)
        axs[j, 0].set_title("w"+str(j))
        axs[j, 1].plot(b_vl, label=trg)
        axs[j, 1].set_title("b"+str(j))
        axs[j, 0].grid(True)
        axs[j, 1].grid(True)
        j += 1
    f.suptitle("Grad sup/hint:" + tag, fontsize=8)
    plt.legend()
    f.savefig(fold_exp+"/grad_" + trg + ".png", bbox_inches='tight')
    plt.close("all")
    del f

def plot_debug_ratio_grad(debug, fold_exp, r="h/s"):
    plt.close("all")
    # f = plt.figure(figsize=(15, 10.8), dpi=300)
    nbr_rows = int(len(debug["grad_sup"][0])/2)
    f, axs = plt.subplots(nbr_rows, 2, sharex=True, sharey=False,
                          figsize=(15, 12.8), dpi=300)

    grads = np.array(debug["grad_sup"])
    gradh = np.array(debug["grad_hint"])
    if gradh.size != grads.size:
        print "Can't calculate the ratio. It looks like you divided the " +\
            "hint batch..."
        return 0
    print gradh.shape, grads.shape
    j = 0
    for i in range(0, nbr_rows*2, 2):
        w_vls = grads[:, i]
        b_vls = grads[:, i+1]
        w_vl_h = gradh[:, i]
        b_vlh = gradh[:, i+1]
        if r == "h/s":
            ratio_w = np.divide(w_vl_h, w_vls)
            ratio_b = np.divide(b_vlh, b_vls)
        elif r == "s/h":
            ratio_w = np.divide(w_vls, w_vl_h)
            ratio_b = np.divide(b_vls, b_vlh)
        else:
            raise ValueError("Either h/s or s/h.")
        axs[j, 0].plot(ratio_w, label=r)
        axs[j, 0].set_title("w"+str(j))
        axs[j, 1].plot(ratio_b, label=r)
        axs[j, 1].set_title("b"+str(j))
        axs[j, 0].grid(True)
        axs[j, 1].grid(True)
        j += 1
    f.suptitle("Ratio gradient: " + r, fontsize=8)
    plt.legend()
    f.savefig(fold_exp+"/ratio_grad_" + r.replace("/", "-") + ".png",
              bbox_inches='tight')
    plt.close("all")
    del f

def plot_validation_cost(train_error, val_error, class_rate=None, savefilename=None):
    epochs = range(len(train_error))
    fig, ax1 = plt.subplots()
    ax1.plot(epochs, train_error, label='train error')
    ax1.plot(epochs, val_error, label='validation error')
    ax1.set_xlabel('epoch')
    ax1.set_ylabel('cost')
    plt.title('Validation Cost')
    lines = ax1.get_lines()
    # Shrink current axis's height by 10% on the bottom
    box = ax1.get_position()
    ax1.set_position([box.x0, box.y0 + box.height * 0.1,
                      box.width, box.height * 0.9])
    if class_rate is not None:
        ax2 = plt.twinx(ax1)
        ax2.plot(epochs, class_rate, label='classification rate', color='r')
        ax2.set_ylabel('classification rate')
        lines.extend(ax2.get_lines())
        ax2.set_position([box.x0, box.y0 + box.height * 0.1,
                          box.width, box.height * 0.9])

    labels = [l.get_label() for l in lines]
    # Put a legend below current axis
    ax1.legend(lines, labels, loc='upper center', bbox_to_anchor=(0.5, -0.05),
               fancybox=False, shadow=False, ncol=5)
    # ax1.legend(lines, labels, loc='lower right')
    if savefilename:
        plt.savefig(savefilename)
    plt.show()

def show_clusters_over_time(
        task_output_path=None,
        query_laps=[0, 1, 2, 5, 10, None],
        nrows=2):
    ''' Read model snapshots from provided folder and make visualizations

    Post Condition
    --------------
    New matplotlib plot with some nice pictures.
    '''
    ncols = int(np.ceil(len(query_laps) // float(nrows)))
    fig_handle, ax_handle_list = pylab.subplots(
        figsize=(FIG_SIZE[0] * ncols, FIG_SIZE[1] * nrows),
        nrows=nrows, ncols=ncols, sharex=True, sharey=True)
    for plot_id, lap_val in enumerate(query_laps):
        cur_model, lap_val = bnpy.load_model_at_lap(task_output_path, lap_val)
        # Plot the current model
        cur_ax_handle = ax_handle_list.flatten()[plot_id]
        bnpy.viz.PlotComps.plotCompsFromHModel(
            cur_model, Data=dataset, ax_handle=cur_ax_handle)
        cur_ax_handle.set_xticks([-2, -1, 0, 1, 2])
        cur_ax_handle.set_yticks([-2, -1, 0, 1, 2])
        cur_ax_handle.set_xlabel("lap: %d" % lap_val)
    pylab.tight_layout()


###############################################################################
# Training from K=1 cluster
# -------------------------
# 
# Using 1 initial cluster, with birth and merge proposal moves.