diff --git a/swp2.py b/swp2.py index f9703ee..313dc69 100644 --- a/swp2.py +++ b/swp2.py @@ -7,7 +7,7 @@ from matplotlib.backends.backend_pdf import PdfPages from matplotlib.ticker import MaxNLocator import io from mpl_toolkits.axes_grid1.inset_locator import inset_axes - +from matplotlib.ticker import MultipleLocator def log_facto(k): k = int(k) if k > 1e6: @@ -259,10 +259,11 @@ def plot_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc # intialize figure my_dpi = 300 fnt_size = 18 + # plt.rcParams['font.size'] = fnt_size fig, ax1 = plt.subplots(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) else: fnt_size = 12 - plt.rcParams['font.size'] = fnt_size + # plt.rcParams['font.size'] = fnt_size ax1 = ax[0,0] #ax1.set_xlim(1e-3, 1) ax1.set_yscale('log') @@ -293,26 +294,30 @@ def plot_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc best10_plots = sorted(plots_likelihoods)[:10] top_plot_lik = str(best10_plots[0]) plot_handles = [] + # plt.rcParams['font.size'] = fnt_size p0, = ax1.plot(top_plots[top_plot_lik][0], top_plots[top_plot_lik][1], 'o', linestyle = "-", - alpha=1, lw=2, label = str(top_plots[top_plot_lik][2])+' epoch | Lik='+top_plot_lik) + alpha=1, lw=2, label = str(top_plots[top_plot_lik][2])+' brks | Lik='+top_plot_lik) plot_handles.append(p0) for k, plot_Lk in enumerate(best10_plots[1:]): plot_Lk = str(plot_Lk) + # plt.rcParams['font.size'] = fnt_size p, = ax1.plot(top_plots[plot_Lk][0], top_plots[plot_Lk][1], 'o', linestyle = "--", - alpha=1/(k+1), lw=1.5, label = str(top_plots[plot_Lk][2])+' epoch | Lik='+plot_Lk) + alpha=1/(k+1), lw=1.5, label = str(top_plots[plot_Lk][2])+' brks | Lik='+plot_Lk) plot_handles.append(p) if theta_scale: - ax1.set_xlabel("Coal. time") - ax1.set_ylabel("Pop. size scaled by N0") + ax1.set_xlabel("Coal. time", fontsize=fnt_size) + ax1.set_ylabel("Pop. size scaled by N0", fontsize=fnt_size) # recent_scale_lower_bound = 0.01 # recent_scale_upper_bound = 0.1 # ax1.axvline(x=recent_scale_lower_bound) # ax1.axvline(x=recent_scale_upper_bound) else: # years - plt.set_xlabel("Time (years)") - plt.set_ylabel("Individuals (N)") + plt.set_xlabel("Time (years)", fontsize=fnt_size) + plt.set_ylabel("Individuals (N)", fontsize=fnt_size) ax1.set_xlim(1e-5, 1) + # plt.rcParams['font.size'] = fnt_size + # print(fnt_size, "rcParam font.size=", plt.rcParams['font.size']) ax1.legend(handles = plot_handles, loc='best', fontsize = fnt_size*0.5) ax1.set_title(title) if ax is None: @@ -335,24 +340,25 @@ def plot_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc # basic plot likelihood if ax is None: fig, ax2 = plt.subplots(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) - plt.rcParams['font.size'] = '18' + # plt.rcParams['font.size'] = fnt_size else: - plt.rcParams['font.size'] = fnt_size + #plt.rcParams['font.size'] = fnt_size ax2 = ax[2,0] ax2.plot(np.array(brkpt_lik)[:, 0], np.array(brkpt_lik)[:, 1].astype(float), 'o', linestyle = "dotted", lw=2) ax2.axhline(y=-Ln, linestyle = "-.", color = "red", label = "$-\log\mathcal{L}$ = "+str(round(-Ln, 2))) ax2.set_yscale('log') ax2.set_xlabel("# breakpoints", fontsize=fnt_size) - ax2.set_ylabel("$-\log\mathcal{L}$") - ax2.legend(loc='best', fontsize = fnt_size*0.8) + ax2.set_ylabel("$-\log\mathcal{L}$", fontsize=fnt_size) + ax2.legend(loc='best', fontsize = fnt_size*0.5) ax2.set_title(title+" Likelihood gain from # breakpoints") if ax is None: plt.savefig(title+'_Breakpts_Likelihood.pdf') # AIC if ax is None: fig, ax3 = plt.subplots(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) - plt.rcParams['font.size'] = '18' + # plt.rcParams['font.size'] = '18' else: + #plt.rcParams['font.size'] = fnt_size ax3 = ax[2,1] AIC = 2*(len(brkpt_lik)+1)+2*np.array(brkpt_lik)[:, 1].astype(float) ax3.plot(np.array(brkpt_lik)[:, 0], AIC, 'o', linestyle = "dotted", lw=2) @@ -362,12 +368,14 @@ def plot_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc ax3.set_yscale('log') ax3.set_xlabel("# breakpoints", fontsize=fnt_size) ax3.set_ylabel("AIC") - ax3.legend(loc='best', fontsize = fnt_size*0.8) + ax3.legend(loc='best', fontsize = fnt_size*0.5) ax3.set_title(title+" AIC") if ax is None: plt.savefig(title+'_Breakpts_Likelihood_AIC.pdf') + print("S", S) return ax -def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, breaks_max = 10, ax = None): + +def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, breaks_max = 10, ax = None, n_ticks = 10): """ Use theta values as is to do basic plots. """ @@ -390,16 +398,14 @@ def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, # intialize figure 1 my_dpi = 300 fnt_size = 18 + # plt.rcParams['font.size'] = fnt_size fig, ax1 = plt.subplots(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) - # Add some extra space for the second axis at the bottom - fig.subplots_adjust(bottom=0.15) else: fnt_size = 12 + # plt.rcParams['font.size'] = fnt_size ax1 = ax[0, 1] plt.subplots_adjust(wspace=0.3, hspace=0.3) - twin = ax1.twiny() - plots = [] for epoch, theta in epochs.items(): groups = np.array(list(theta.values()), dtype=object)[:, 1].tolist() @@ -421,35 +427,32 @@ def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, prop = prop[::-1] # print(prop, "\n", sum(prop)) # normalise to N0 (N0 of epoch1) + x_ticks = ax1.get_xticks() for i in range(len(y)): y[i] = y[i]/N0 # plot x_plot, y_plot = plot_straight_x_y(x, y) #plt.plot(x, y, 'o', linestyle="dotted", alpha=0.75, lw=2, label = str(epoch)+' brks') - p, = ax1.plot(x_plot, y_plot, 'o', linestyle="-.", alpha=0.75, lw=2, label = str(epoch)+' brks') + p, = ax1.plot(x_plot, y_plot, 'o', linestyle="-", alpha=0.75, lw=2, label = str(epoch)+' brks') # add plot to the list of all plots to superimpose plots.append(p) - # virtual line to get the second x axis for proportions - p0, = twin.plot(prop, y, alpha = 0, label="Proportion") - # Move twinned axis ticks and label from top to bottom - twin.xaxis.set_ticks_position("bottom") - twin.xaxis.set_label_position("bottom") - # Offset the twin axis below the host - if ax is None: - # arrange differently the second x axis if the plot is plain - twin.spines["bottom"].set_position(("axes", -0.15)) - else: - # in a combined plot, more space between the fig and the axis - twin.spines["bottom"].set_position(("axes", -0.35)) + #print(prop, "\n", sum(prop)) #ax.legend(handles=[p0]+plots) - ax1.set_xlabel("# bin") + ax1.set_xlabel("# bin", fontsize=fnt_size) # Set the x-axis locator to reduce the number of ticks to 10 - ax1.xaxis.set_major_locator(MaxNLocator(nbins=10)) - twin.xaxis.set_major_locator(MaxNLocator(nbins=10)) - ax1.set_ylabel("theta") - twin.set_ylabel("Proportion") - ax1.set_title("Title") + ax1.set_ylabel("theta", fontsize=fnt_size) + ax1.set_title("Title", fontsize=fnt_size) ax1.legend(handles=plots, loc='best', fontsize = fnt_size*0.5) + ax1.set_xticks(x_ticks) + if len(prop) >= 18: + ax1.locator_params(nbins=n_ticks) + # new scale of ticks if too many values + cumul = 0 + prop_cumul = [] + for val in prop: + prop_cumul.append(val+cumul) + cumul = val+cumul + ax1.set_xticklabels([f'{x[k]}\n{val:.2f}' for k, val in enumerate(prop_cumul)]) if ax is None: plt.savefig(title+'_raw'+str(k)+'.pdf') # fig 2 & 3 @@ -457,6 +460,7 @@ def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, fig2, ax2 = plt.subplots(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) fig3, ax3 = plt.subplots(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) else: + # plt.rcParams['font.size'] = fnt_size # place of plots on the grid ax2 = ax[1,0] ax3 = ax[1,1] @@ -487,43 +491,48 @@ def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, x_2 = [0]+x_2 y = [y[0]]+y x2_plot, y2_plot = plot_straight_x_y(x_2, y) - p2, = ax2.plot(x2_plot, y2_plot, 'o', linestyle="-.", alpha=0.75, lw=2, label = str(epoch)+' brks') + p2, = ax2.plot(x2_plot, y2_plot, 'o', linestyle="-", alpha=0.75, lw=2, label = str(epoch)+' brks') lines_fig2.append(p2) # Plotting (fig 3) which is the same but log scale for x - p3, = ax3.plot(x2_plot, y2_plot, 'o', linestyle="-.", alpha=0.75, lw=2, label = str(epoch)+' brks') + p3, = ax3.plot(x2_plot, y2_plot, 'o', linestyle="-", alpha=0.75, lw=2, label = str(epoch)+' brks') lines_fig3.append(p3) - ax2.set_xlabel("Relative scale") - ax2.set_ylabel("theta") - ax2.set_title("Title") + ax2.set_xlabel("Relative scale", fontsize=fnt_size) + ax2.set_ylabel("theta", fontsize=fnt_size) + ax2.set_title("Title", fontsize=fnt_size) ax2.legend(handles=lines_fig2, loc='best', fontsize = fnt_size*0.5) if ax is None: plt.savefig(title+'_plot2_'+str(k)+'.pdf') ax3.set_xscale('log') - ax3.set_xlabel("log Relative scale") - ax3.set_ylabel("theta") - ax3.set_title("Title") + ax3.set_yscale('log') + ax3.set_xlabel("log Relative scale", fontsize=fnt_size) + ax3.set_ylabel("theta", fontsize=fnt_size) + ax3.set_title("Title", fontsize=fnt_size) ax3.legend(handles=lines_fig3, loc='best', fontsize = fnt_size*0.5) if ax is None: plt.savefig(title+'_plot3_'+str(k)+'_log.pdf') + plt.clf() # return plots return ax def combined_plot(folder_path, mu, tgen, breaks, title = "Title", theta_scale = True): my_dpi = 300 # Add some extra space for the second axis at the bottom + #plt.rcParams['font.size'] = 18 fig, axs = plt.subplots(3, 2, figsize=(5000/my_dpi, 2970/my_dpi), dpi=my_dpi) + #plt.rcParams['font.size'] = 12 ax = plot_all_epochs_thetafolder(folder_path, mu, tgen, title, theta_scale, ax = axs) ax = plot_test_theta(folder_path, mu, tgen, title, theta_scale, breaks_max = breaks, ax = axs) # Adjust layout to prevent clipping of titles plt.tight_layout() # Adjust absolute space between the top and bottom rows - plt.subplots_adjust(hspace=0.7) # Adjust this value based on your requirement + #plt.subplots_adjust(hspace=0.7) # Adjust this value based on your requirement # Save the entire grid as a single figure plt.savefig(title+'_combined.pdf') - plt.close() + plt.clf() # second call for individual plots plot_all_epochs_thetafolder(folder_path, mu, tgen, title, theta_scale, ax = None) plot_test_theta(folder_path, mu, tgen, title, theta_scale, breaks_max = breaks, ax = None) + plt.clf() if __name__ == "__main__":