python类subplots()的实例源码

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.

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()

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

def show_clusters_over_time(
        task_output_path=None,
        query_laps=[0, 1, 2, 10, 20, None],
        nrows=2):
    '''
    '''
    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)
        cur_ax_handle = ax_handle_list.flatten()[plot_id]
        bnpy.viz.PlotComps.plotCompsFromHModel(
            cur_model, dataset=dataset, ax_handle=cur_ax_handle)
        cur_ax_handle.set_xlim([-4.5, 4.5])
        cur_ax_handle.set_xlabel("lap: %d" % lap_val)
    pylab.tight_layout()


###############################################################################
#
# Run with *merge* moves only, from K=5 initial clusters
# --------------------------------------------------------
#
# Unfortunately, no pairwise merge is accepted.
# The model is stuck using 5 clusters when one cluster would do.

def show_top_words_over_time(
        task_output_path=None,
        vocabList=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)
        # Plot the current model
        cur_ax_list = ax_handles_RC[row_id].flatten().tolist()
        bnpy.viz.PrintTopics.plotCompsFromHModel(
            cur_model,
            vocabList=vocabList,
            fontsize=9,
            Ktop=7,
            ax_list=cur_ax_list)
        cur_ax_list[0].set_ylabel("lap: %d" % lap_val)
    pylab.subplots_adjust(
        wspace=0.04, hspace=0.1, 
        left=0.01, right=0.99, top=0.99, bottom=0.1)
    pylab.tight_layout()


###############################################################################
#
# Show the topics over time

def show_top_words_over_time(
        task_output_path=None,
        vocabList=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)
        # Plot the current model
        cur_ax_list = ax_handles_RC[row_id].flatten().tolist()
        bnpy.viz.PrintTopics.plotCompsFromHModel(
            cur_model,
            vocabList=vocabList,
            fontsize=9,
            Ktop=7,
            ax_list=cur_ax_list)
        cur_ax_list[0].set_ylabel("lap: %d" % lap_val)
    pylab.subplots_adjust(
        wspace=0.04, hspace=0.1, 
        left=0.01, right=0.99, top=0.99, bottom=0.1)
    pylab.tight_layout()


###############################################################################
#
# Show the topics over time

def show_clusters_over_time(
        task_output_path=None,
        query_laps=[0, 1, 2, 10, 20, None],
        nrows=2):
    '''
    '''
    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)
        cur_ax_handle = ax_handle_list.flatten()[plot_id]
        bnpy.viz.PlotComps.plotCompsFromHModel(
            cur_model, dataset=dataset, ax_handle=cur_ax_handle)
        cur_ax_handle.set_title("lap: %d" % lap_val)
        cur_ax_handle.set_xlabel(dataset.column_names[0])
        cur_ax_handle.set_ylabel(dataset.column_names[1])
        cur_ax_handle.set_xlim(data_ax_h.get_xlim())
        cur_ax_handle.set_ylim(data_ax_h.get_ylim())
    pylab.tight_layout()



###############################################################################
#
# *DiagGauss* observation model, without moves
# --------------------------------------------
#
# Start with too many clusters (K=25)

def show_clusters_over_time(
        task_output_path=None,
        query_laps=[0, 1, 2, 10, 20, None],
        nrows=2):
    ''' Show 2D elliptical contours overlaid on raw data.
    '''
    ncols = int(np.ceil(len(query_laps) // float(nrows)))
    fig_handle, ax_handle_list = pylab.subplots(
        figsize=(SMALL_FIG_SIZE[0] * ncols, SMALL_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)
        cur_ax_handle = ax_handle_list.flatten()[plot_id]
        bnpy.viz.PlotComps.plotCompsFromHModel(
            cur_model, dataset=dataset, ax_handle=cur_ax_handle)
        cur_ax_handle.set_title("lap: %d" % lap_val)
        cur_ax_handle.set_xlabel(dataset.column_names[0])
        cur_ax_handle.set_ylabel(dataset.column_names[1])
        cur_ax_handle.set_xlim(data_ax_h.get_xlim())
        cur_ax_handle.set_ylim(data_ax_h.get_ylim())
    pylab.tight_layout()


###############################################################################
#
# *DiagGauss* observation model
# -----------------------------
#
# Assume diagonal covariances.
#
# Start with too many clusters (K=20)

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()

###############################################################################
# Training from K=3 with births
# -----------------------------
# 
# Initialization: 3 topics, using random initial guess

def show_bars_over_time(
        task_output_path=None,
        query_laps=[0, 1, 2, 5, None],
        ncols=10):
    ''' Show square-image visualization of estimated topics over time.

    Post Condition
    --------------
    New matplotlib figure with visualization (one row per lap).
    '''
    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()

###############################################################################
#
# Examine the bars 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()

###############################################################################
# Train LDA topic model
# ---------------------
# 
# Using 10 clusters and the 'randexamples' initialization procedure.

def _viz_Gauss_before_after(
        curModel=None, propModel=None,
        curSS=None, propSS=None,
        Plan=None,
        propLscore=None, curLscore=None,
        Data_b=None, Data_t=None, 
        **kwargs):
    pylab.subplots(
        nrows=1, ncols=2, figsize=(8, 4), num=1)
    h1 = pylab.subplot(1, 2, 1)
    h1.clear()
    GaussViz.plotGauss2DFromHModel(
        curModel, compsToHighlight=Plan['btargetCompID'], figH=h1)
    if curLscore is not None:
        pylab.title('%.4f' % (curLscore))

    h2 = pylab.subplot(1, 2, 2, sharex=h1, sharey=h1)
    h2.clear()
    newCompIDs = np.arange(curModel.obsModel.K, propModel.obsModel.K)
    GaussViz.plotGauss2DFromHModel(
        propModel, compsToHighlight=newCompIDs, figH=h2, Data=Data_t)
    if propLscore is not None:
        pylab.title('%.4f' % (propLscore))

        Lgain = propLscore - curLscore
        if Lgain > 0:
            pylab.xlabel('ACCEPT +%.2f' % (Lgain))
        else:
            pylab.xlabel('REJECT %.2f' % (Lgain))
    pylab.draw()
    pylab.subplots_adjust(hspace=0.1, top=0.9, bottom=0.15,
                          left=0.15, right=0.95)

def plot_colat_slice(self, component, colat, valmin, valmax, iteration=0, verbose=True):

        #- Some initialisations. ------------------------------------------------------------------

        colat = np.pi * colat / 180.0

        n_procs = self.setup["procs"]["px"] * self.setup["procs"]["py"] * self.setup["procs"]["pz"]

        vmax = float("-inf")
        vmin = float("inf")

        fig, ax = plt.subplots()

        #- Loop over processor boxes and check if colat falls within the volume. ------------------
        for p in range(n_procs):

            if (colat >= self.theta[p,:].min()) & (colat <= self.theta[p,:].max()):

                #- Read this field and make lats & lons. ------------------------------------------
                field = self.read_single_box(component,p,iteration)

                r, lon = np.meshgrid(self.z[p,:], self.phi[p,:])

                x = r * np.cos(lon)
                y = r * np.sin(lon)

                #- Find the colat index and plot for this one box. --------------------------------
                idx=min(np.where(min(np.abs(self.theta[p,:]-colat))==np.abs(self.theta[p,:]-colat))[0])

                colat_effective = self.theta[p,idx]*180.0/np.pi

                #- Find min and max values. -------------------------------------------------------

                vmax = max(vmax, field[idx,:,:].max())
                vmin = min(vmin, field[idx,:,:].min())

                #- Make a nice colourmap and plot. ------------------------------------------------
                my_colormap=cm.make_colormap({0.0:[0.1,0.0,0.0], 0.2:[0.8,0.0,0.0], 0.3:[1.0,0.7,0.0],0.48:[0.92,0.92,0.92], 0.5:[0.92,0.92,0.92], 0.52:[0.92,0.92,0.92], 0.7:[0.0,0.6,0.7], 0.8:[0.0,0.0,0.8], 1.0:[0.0,0.0,0.1]})

                cax = ax.pcolor(x, y, field[idx,:,:], cmap=my_colormap, vmin=valmin,vmax=valmax)


        #- Add colobar and title. ------------------------------------------------------------------
        cbar = fig.colorbar(cax)
        if component in UNIT_DICT:
            cb.set_label(UNIT_DICT[component], fontsize="x-large", rotation=0)

        plt.suptitle("Vertical slice of %s at %i degree colatitude" % (component, colat_effective), size="large")

        plt.axis('equal')
        plt.show()


    #==============================================================================================
    #- Plot depth slice.
    #==============================================================================================

def main(args_dict):
    # Extract configuration from command line arguments
    MK = args_dict['MK']
    Kmin = args_dict['Kmin']

    # Load data
    data = json_load('data/astar_rbr_MK%d.json' % (MK))
    lnZ = data['lnZ']
    MAPs = np.array(data['MAPs'])
    print('Loaded %d MAP samples from A* sampling' % (len(MAPs)))

    # Estimate MSE of lnZ estimators from Gumbel and Exponential tricks
    MSEs_Gumb = []
    MSEs_Expo = []
    Ms = xrange(1, MK / Kmin)
    for M in Ms:
        # Computation with M samples, repeated K >= Kmin times with a new set every time
        K = MK / M
        myMAPs = np.reshape(MAPs[:(K*M)], (K, M))

        # Compute unbiased estimators of ln(Z)
        lnZ_Gumb = np.mean(myMAPs, axis=1)
        lnZ_Expo = EULER - np.log(np.mean(np.exp(- myMAPs), axis=1)) - (np.log(M) - digamma(M))

        # Save MSE estimates
        MSEs_Gumb.append(np.mean((lnZ_Gumb - lnZ) ** 2))
        MSEs_Expo.append(np.mean((lnZ_Expo - lnZ) ** 2))

    # Set up plot
    matplotlib_configure_as_notebook()
    fig, ax = plt.subplots(1, 1, facecolor='w', figsize=(4.25, 3.25))
    ax.set_xscale('log')
    ax.set_xlabel('desired MSE (lower to the right)')
    ax.set_ylabel('required number of samples $M$')
    ax.grid(b=True, which='both', linestyle='dotted', lw=0.5, color='black', alpha=0.3)

    # Plot MSEs
    ax.plot(MSEs_Gumb, Ms, color=tableau20(0), label='Gumbel')
    ax.plot(MSEs_Expo, Ms, color=tableau20(2), label='Exponential')

    # Finalize plot
    ax.set_xlim((1e-2, 2))
    ax.invert_xaxis()
    lgd = ax.legend(loc='upper left')
    save_plot(fig, 'figures/fig3a', (lgd,))

def gaussian_twoD_testing():
    """ Implement and check the estimator for a 2D gaussian fit. """

    data = np.empty((121,1))
    amplitude=np.random.normal(3e5,1e5)
    center_x=91+np.random.normal(0,0.8)
    center_y=14+np.random.normal(0,0.8)
    sigma_x=np.random.normal(0.7,0.2)
    sigma_y=np.random.normal(0.7,0.2)
    offset=0
    x = np.linspace(90,92,11)
    y = np.linspace(13,15,12)
    xx, yy = np.meshgrid(x, y)

    axes=(xx.flatten(), yy.flatten())

    theta_here=10./360.*(2*np.pi)

#            data=qudi_fitting.twoD_gaussian_function((xx,yy),*(amplitude,center_x,center_y,sigma_x,sigma_y,theta_here,offset))
    gmod,params = qudi_fitting.make_twoDgaussian_model()

    data = gmod.eval(x=axes, amplitude=amplitude, center_x=center_x,
                     center_y=center_y, sigma_x=sigma_x, sigma_y=sigma_y,
                     theta=theta_here, offset=offset)
    data += 50000*np.random.random_sample(np.shape(data))

    gmod, params = qudi_fitting.make_twoDgaussian_model()

    para=Parameters()
#    para.add('theta',vary=False)
#    para.add('center_x',expr='0.5*center_y')
#    para.add('sigma_x',min=0.2*((92.-90.)/11.), max=   10*(x[-1]-y[0]) )
#    para.add('sigma_y',min=0.2*((15.-13.)/12.), max=   10*(y[-1]-y[0]))
#    para.add('center_x',value=40,min=50,max=100)

    result = qudi_fitting.make_twoDgaussian_fit(xy_axes=axes, data=data)#,add_parameters=para)

#
#            FIXME: What does "Tolerance seems to be too small." mean in message?
#            print(result.message)
    plt.close('all')

    fig, ax = plt.subplots(1, 1)
    ax.hold(True)
    ax.imshow(result.data.reshape(len(y),len(x)),
              cmap=plt.cm.jet, origin='bottom', extent=(x.min(), x.max(),
                                       y.min(), y.max()),interpolation="nearest")
    ax.contour(x, y, result.best_fit.reshape(len(y),len(x)), 8
                , colors='w')
    plt.show()
#    plt.close('all')

    print(result.fit_report())

#            print('Message:',result.message)

def show_many_random_initial_models(
        obsPriorArgsDict,
        initArgsDict,
        nrows=1, ncols=6):
    ''' Create plot of many different random initializations
    '''
    fig_handle, ax_handle_list = pylab.subplots(
        figsize=(SMALL_FIG_SIZE[0] * ncols, SMALL_FIG_SIZE[1] * nrows),
        nrows=nrows, ncols=ncols, sharex=True, sharey=True)
    for trial_id in range(nrows * ncols):
        cur_model = bnpy.make_initialized_model(
            dataset,
            allocModelName='FiniteMixtureModel',
            obsModelName='Gauss',
            algName='VB',
            allocPriorArgsDict=dict(gamma=10.0),
            obsPriorArgsDict=obsPriorArgsDict,
            initArgsDict=initArgsDict,
            seed=int(trial_id),
            )
        # Plot the current model
        cur_ax_handle = ax_handle_list.flatten()[trial_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])
    pylab.tight_layout()



###############################################################################
# initname: 'randexamples'
# ------------------------
# This procedure selects K examples uniformly at random.
# Each cluster is then initialized from one selected example,
# using a standard global step update. 
#
# **Example 1**:
# Initialize with 8 clusters, with prior biased towards small covariances
#
# .. math::
#
#   \E_{\mbox{prior}}[ \Sigma_k ] = 0.01 I_D