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.
python类subplots()的实例源码
plot-03-demo=vb+proposals-model=dp_mix+gauss.py 文件源码
项目:bnpy
作者: bnpy
项目源码
文件源码
阅读 28
收藏 0
点赞 0
评论 0
plot-02-demo=vb_single_run-model=dp_mix+gauss.py 文件源码
项目:bnpy
作者: bnpy
项目源码
文件源码
阅读 28
收藏 0
点赞 0
评论 0
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.
run-04-demo=topic_vb+proposals-model=hdp_topic+mult.py 文件源码
项目:bnpy
作者: bnpy
项目源码
文件源码
阅读 38
收藏 0
点赞 0
评论 0
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
run-03-demo=topic_vb_single_run-model=hdp_topic+mult.py 文件源码
项目:bnpy
作者: bnpy
项目源码
文件源码
阅读 38
收藏 0
点赞 0
评论 0
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
plot-02-demo=merge_and_delete-model=dp_mix+gauss.py 文件源码
项目:bnpy
作者: bnpy
项目源码
文件源码
阅读 23
收藏 0
点赞 0
评论 0
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)
run-05-demo=topic_model_vb+births-model=hdp_topic+mult.py 文件源码
项目:bnpy
作者: bnpy
项目源码
文件源码
阅读 30
收藏 0
点赞 0
评论 0
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
run-04-demo=topic_model_vb+merges-model=hdp_topic+mult.py 文件源码
项目:bnpy
作者: bnpy
项目源码
文件源码
阅读 35
收藏 0
点赞 0
评论 0
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
run-03-demo=topic_model_vb_single_run-model=hdp_topic+mult.py 文件源码
项目:bnpy
作者: bnpy
项目源码
文件源码
阅读 32
收藏 0
点赞 0
评论 0
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