From 33a120ed7e2bc0da9f0d74d4e837461dda9edb59 Mon Sep 17 00:00:00 2001 From: richard Date: Tue, 25 Aug 2026 10:26:30 +0200 Subject: [PATCH] =?UTF-8?q?version=20utilis=C3=A9e=20pour=20la=20th=C3=A8s?= =?UTF-8?q?e=20/=20PCIA?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- dependences/pybam.py | 12 +-- sfs_tools.py | 25 +++++- swp2.py | 180 ++++++++++++++++++++++++++++++++----------- 3 files changed, 161 insertions(+), 56 deletions(-) diff --git a/dependences/pybam.py b/dependences/pybam.py index 9f02c02..0f132a1 100644 --- a/dependences/pybam.py +++ b/dependences/pybam.py @@ -114,7 +114,7 @@ Github: http://github.com/JohnLonginotto/pybam my_bam = pybam.read('/my/data.bam',decompressor='internal') [ Parse Words (hah) ]''' -wat += '\n'+''.join([('\n===============================================================================================\n\n ' if code is 'file_alignments_read' or code is 'sam' else ' ')+(code+' ').ljust(25,'-')+description+'\n' for code,description in sorted(parse_codes.items())]) + '\n' +wat += '\n'+''.join([('\n===============================================================================================\n\n ' if code == 'file_alignments_read' or code == 'sam' else ' ')+(code+' ').ljust(25,'-')+description+'\n' for code,description in sorted(parse_codes.items())]) + '\n' class read(): ''' @@ -155,7 +155,7 @@ class read(): if fields is not False: print(fields) - if type(fields) is not list or len(fields) is 0: + if type(fields) is not list or len(fields) == 0: raise PybamError('\n\nFields for the static parser must be provided as a non-empty list. You gave a ' + str(type(fields)) + '\n') else: for field in fields: @@ -167,7 +167,7 @@ class read(): if decompressor: if type(decompressor) is str: - if decompressor is not 'internal' and '{}' not in decompressor: raise PybamError('\n\nWhen a custom decompressor is used and the input file is a string, the decompressor string must contain at least one occurence of "{}" to be substituted with a filepath by pybam.\n') + if decompressor != 'internal' and '{}' not in decompressor: raise PybamError('\n\nWhen a custom decompressor is used and the input file is a string, the decompressor string must contain at least one occurence of "{}" to be substituted with a filepath by pybam.\n') else: raise PybamError('\n\nUser-supplied decompressor must be a string that when run on the command line decompresses a named file (or stdin), to stdout:\ne.g. "lzma --decompress --stdout {}" if pybam is provided a path as input file, where {} is substituted for that path.\nor just "lzma --decompress --stdout" if pybam is provided a file object instead of a file path, as data from that file object will be piped via stdin to the decompression program.\n') ## First we make a generator that will return chunks of uncompressed data, regardless of how we choose to decompress: @@ -203,7 +203,7 @@ class read(): elif magic == b"\x1f\x8b\x08\x04": # The user has passed us compressed gzip/bgzip data, which is typical for a BAM file # use custom decompressor if provided: - if decompressor is not False and decompressor is not 'internal': + if decompressor != False and decompressor != 'internal': if type(f) is str: self._subprocess = subprocess.Popen( decompressor.replace('{}',f), shell=True, stdout=subprocess.PIPE, stderr=DEVNULL) else: self._subprocess = subprocess.Popen('{ printf "'+magic+'"; cat; } | ' + decompressor, stdin=self._file, shell=True, stdout=subprocess.PIPE, stderr=DEVNULL) self.file_decompressor = decompressor @@ -230,7 +230,7 @@ class read(): use = 'gzip' except OSError: use = 'internal' - if use is not 'internal' and decompressor is not 'internal': + if use != 'internal' and decompressor != 'internal': if type(f) is str: self._subprocess = subprocess.Popen([ use , '--decompress','--stdout', f ], stdout=subprocess.PIPE, stderr=DEVNULL) else: self._subprocess = subprocess.Popen('{ printf "'+magic+'"; cat; } | ' + use + ' --decompress --stdout', stdin=f, shell=True, stdout=subprocess.PIPE, stderr=DEVNULL) time.sleep(1) @@ -307,7 +307,7 @@ class read(): yield b''.join(internal_cache) return - elif decompressor is not False and decompressor is not 'internal': + elif decompressor != False and decompressor != 'internal': # It wouldn't be safe to just print to the shell four random bytes from the beginning of a file, so instead it's # written to a temp file and cat'd. The idea here being that we trust the decompressor string as it was written by # someone with access to python, so it has system access anyway. The file/data, however, should not be trusted. diff --git a/sfs_tools.py b/sfs_tools.py index 481f115..26a5515 100755 --- a/sfs_tools.py +++ b/sfs_tools.py @@ -241,7 +241,9 @@ def sfs_from_parsed_vcf(n, vcf_dict, folded = True, diploid = True, phased = Fal return SFS_values, count_pluriall -def barplot_sfs(sfs, xlab, ylab, folded=True, title = "Barplot", transformed = False, normalized = False, ploidy = 2, output = None): + +def barplot_sfs(sfs, xlab, ylab, folded=True, title = "Barplot", transformed = False, + normalized = False, ploidy = 2, output = None, step = 8): sfs_val = [] n = len(sfs.values()) sum_sites = sum(list(sfs.values())) @@ -300,8 +302,8 @@ def barplot_sfs(sfs, xlab, ylab, folded=True, title = "Barplot", transformed = if transformed: print("Transformed SFS ( n =",n_title, ") :", sfs_val) - #plt.axhline(y=1/n, color='r', linestyle='-') - plt.bar([x+0.2 for x in list(sfs.keys())], [1/n]*n, fill=False, hatch="///", linestyle='-', width = 0.4, label= "H0 Theoric constant") + plt.axhline(y=1/n, color='r', linestyle='-', label= "H0 Theoric constant") + #plt.bar([x+0.2 for x in list(sfs.keys())], [1/n]*n, fill=False, hatch="///", linestyle='-', width = 0.4, label= "H0 Theoric constant") else: if normalized: @@ -311,12 +313,27 @@ def barplot_sfs(sfs, xlab, ylab, folded=True, title = "Barplot", transformed = print(expected_y) plt.bar([x+0.2 for x in list(sfs.keys())], expected_y, fill=False, hatch="///", linestyle='-', width = 0.4, label= "H0 Theoric constant") print(sum(expected_y)) + + # graphics parameters if output is not None: # if write in a file, don't open the window dynamically plot = False else: plot = True - customgraphics.barplot(x = [x-0.2 for x in X_axis], width=0.4, y= sfs_val, xlab = xlab, ylab = ylab, title = title, label = "H1 Observed spectrum", xticks =list(sfs.keys()), plot = plot ) + # customgraphics.barplot(x = [x-0.2 for x in X_axis], width=0.4, y= sfs_val, xlab = xlab, ylab = ylab, title = title, label = "H1 Observed spectrum", xticks =list(sfs.keys()), plot = plot ) + customgraphics.barplot(x = X_axis, width=0.4, y= sfs_val, xlab = xlab, ylab = ylab, title = title, label = "H1 Observed spectrum", xticks =list(sfs.keys()), plot = plot ) + + + + + # graphics parameters + if n > 12: + # plt.locator_params(axis='x', nbins=10) + step = n//step + x_ticks = [k for k in range(1,n,step)]+[n] + print(x_ticks) + plt.xticks(x_ticks) + if output: plt.savefig(f"{output}/{original_title}_SFS.pdf") else: diff --git a/swp2.py b/swp2.py index 7c6893b..c8e20ff 100644 --- a/swp2.py +++ b/swp2.py @@ -78,18 +78,20 @@ def parse_stwp_theta_file(stwp_theta_file, breaks, mu, tgen, relative_theta_scal j = int(group.split(",")[-1]) t[i] = 0 #t = + #T += y[i]*2 / (x[i]*(x[i]-1)) + if len(group.split(',')) == 1: k = i if relative_theta_scale: - t[i] += ((theta_L[group_nb] ) / (k*(k-1))) + t[i] += ((theta_L[group_nb]*2 ) / (k*(k-1))) else: - t[i] += ((theta_L[group_nb] ) / (k*(k-1)) * tgen) / mu + t[i] += ((theta_L[group_nb]*2 ) / (k*(k-1)) ) / mu else: for k in range(j, i-1, -1 ): if relative_theta_scale: - t[i] += ((theta_L[group_nb] ) / (k*(k-1))) + t[i] += ((theta_L[group_nb]*2 ) / (k*(k-1))) else: - t[i] += ((theta_L[group_nb] ) / (k*(k-1)) * tgen) / mu + t[i] += ((theta_L[group_nb]*2 ) / (k*(k-1)) ) / mu # we add the cumulative times at the end t[i] += sum_t sum_t = t[i] @@ -126,7 +128,7 @@ def plot_straight_x_y(x,y): def plot_all_epochs_thetafolder(full_dict, mu, tgen, title = "Title", theta_scale = True, ax = None, input = None, output = None): - my_dpi = 500 + my_dpi = 800 L = full_dict["L"] if ax is None: # intialize figure @@ -143,14 +145,14 @@ def plot_all_epochs_thetafolder(full_dict, mu, tgen, title = "Title", plot_handles = [] best_plot = full_dict['all_epochs']['best'] p0, = ax1.plot(best_plot[0], best_plot[1], linestyle = "-", - alpha=1, lw=2, label = str(best_plot[2])+' brks | Lik='+best_plot[3]) + alpha=1, lw=2, label = str(best_plot[2])+' brks | Lik='+str(round(float(best_plot[3]), 4))) plot_handles.append(p0) #ax1.grid(True,which="both", linestyle='--', alpha = 0.3) for k, plot_Lk in enumerate(full_dict['all_epochs']['plots']): plot_Lk = str(full_dict['all_epochs']['plots'][k][3]) # plt.rcParams['font.size'] = fnt_size p, = ax1.plot(full_dict['all_epochs']['plots'][k][0], full_dict['all_epochs']['plots'][k][1], linestyle = "-", - alpha=1/(k+1), lw=1.5, label = str(full_dict['all_epochs']['plots'][k][2])+' brks | Lik='+plot_Lk) + alpha=1/(k+1), lw=1.5, label = str(full_dict['all_epochs']['plots'][k][2])+' brks | Lik='+str(round(float(plot_Lk), 4))) plot_handles.append(p) if theta_scale: ax1.set_xlabel("Coal. time", fontsize=fnt_size) @@ -172,11 +174,12 @@ def plot_all_epochs_thetafolder(full_dict, mu, tgen, title = "Title", # ax1.set_xticklabels([f'{k}\n{k/(mu)}\n{k/(mu)*tgen}' for k in x_ticks], fontsize = fnt_size*0.8) # 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.legend(handles = plot_handles, fontsize = fnt_size*0.5, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.) ax1.set_title(title) breaks = len(full_dict['all_epochs']['plots']) if ax is None: - plt.savefig(title+'_best_'+str(breaks+1)+'_epochs.pdf') + fig.tight_layout() + plt.savefig(title+'_best_'+str(breaks+1)+'_epochs.pdf', bbox_inches='tight') # plot likelihood against nb of breakpoints if ax is None: fig, ax2 = plt.subplots(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) @@ -192,12 +195,13 @@ def plot_all_epochs_thetafolder(full_dict, mu, tgen, title = "Title", ax2.scatter(full_dict['Ln_Brks'][0], full_dict['Ln_Brks'][1], s=50, c=colors, marker='o', zorder=2) ax2.axhline(y=full_dict['best_Ln'], linestyle = "-.", color = "red", label = "$-\log\mathcal{L}$ = "+str(round(full_dict['best_Ln'], 2))) ax2.set_yscale('log') - ax2.set_xlabel("# breakpoints", fontsize=fnt_size) + ax2.set_xlabel("# breakpoints", 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') + fig.tight_layout() + plt.savefig(title+'_Breakpts_Likelihood.pdf', bbox_inches='tight') # AIC if ax is None: fig, ax3 = plt.subplots(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) @@ -212,11 +216,12 @@ def plot_all_epochs_thetafolder(full_dict, mu, tgen, title = "Title", ax3.axhline(y=full_dict['best_AIC'], linestyle = "-.", color = "red", label = "Min. AIC = "+str(round(full_dict['best_AIC'], 2))) ax3.set_yscale('log') - ax3.set_xlabel("# breakpoints", fontsize=fnt_size) - ax3.set_ylabel("AIC") + ax3.set_xlabel("# breakpoints", fontsize=fnt_size*0.8) + ax3.set_ylabel("AIC", fontsize=fnt_size*0.8) ax3.legend(loc='best', fontsize = fnt_size*0.5) ax3.set_title(title+" AIC") if ax is None: + fig.tight_layout() plt.savefig(title+'_Breakpts_Likelihood_AIC.pdf') else: # return plots @@ -307,7 +312,7 @@ def save_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc for i in range(2, my_n): an +=1.0/i - print("an=", an, "theta_w", S/an, "theta_w_p_site", (S/an)/L) + #print("an=", an, "theta_w", S/an, "theta_w_p_site", (S/an)/L) # compute Ln Ln = log_facto(S+S0) - log_facto(S0) + np.log(float(S0)/(S+S0)) * S0 for xi in range(0, len(SFS_stored)): @@ -319,10 +324,13 @@ def save_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc AIC = [] for brk in np.array(brkpt_lik)[:, 0]: brk = int(brk) - AIC.append((2*brk+1)+2*np.array(brkpt_lik)[brk, 1].astype(float)) + AIC.append((2*(brk+1))+(2*np.array(brkpt_lik)[brk, 1].astype(float))) AIC_Brks = [list(np.array(brkpt_lik)[:, 0]), AIC] # AIC = 2*k - 2ln(L) ; where k is the number of parameters, here brks+1 AIC_ln = 2*(len(brkpt_lik)+1) - 2*Ln + # AIC_new = [] + # for brk in range(len(brkpt_lik)): + # AIC_new.append(2*brk+1 - 2*Ln) best_AIC = AIC_ln selected_brks_nb = AIC.index(min(AIC)) # to return : plots ; Ln_Brks ; AIC_Brks ; best_Ln ; best_AIC @@ -463,7 +471,7 @@ def save_k_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, x[i] = int(x[i]) # compute the times as: theta_k / (k*(k-1)) for i in range(0, len(x)): - T += y[i]*2 / (x[i]*(x[i]-1)) + T += y[i]*4 / (x[i]*(x[i]-1)) x_2.append(T) # Save plotting (fig 2) # x_2 = [0]+x_2 @@ -501,8 +509,10 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax = nb_epochs = len(plot_lines) # fig 2 & 3 if ax is None: - my_dpi = 500 + my_dpi = 800 fnt_size = 18 + # make sure that everything is closed before... + plt.close() 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: @@ -519,12 +529,46 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax = swp2_lines[0][k] = swp2_lines[0][k]/tgen for k in range(len(swp2_lines[1])): swp2_lines[1][k] = swp2_lines[1][k] - # x2_plot, y2_plot = plot_straight_x_y(swp2_lines[0],swp2_lines[1]) - x2_plot, y2_plot = swp2_lines[0], swp2_lines[1] - p2, = ax2.plot(x2_plot, y2_plot, linestyle="-", alpha=0.75, lw=2, label = 'swp2', color="black") + #x2_plot, y2_plot = swp2_lines[0], swp2_lines[1] + x2_plot, y2_plot = plot_straight_x_y(swp2_lines[0],swp2_lines[1]) + + title_swp2 = "swp2" + + swp2_x = swp2_lines[0] + swp2_y = swp2_lines[1] + # SIM ONLY + # sim1 + # title_swp2 = "sim1" + # x2_plot = [0, 100, 100, 4600, 4600, 10000] + # y2_plot = [5000, 5000, 25000, 25000, 125000, 125000] + + # sim2 + # title_swp2 = "sim2" + # x2_plot = [0, 100, 100, 280, 280, 10000] + # y2_plot = [5000, 5000, 1000, 1000, 200, 200] + + # # sim3 + # title_swp2 = "sim3" + # x2_plot = [0, 100, 100, 280, 280, 10000] + # y2_plot = [5000, 5000, 1000, 1000, 5000, 5000] + + + # # sim4 + # title_swp2 = "sim4" + # x2_plot = [0, 100, 100, 4600, 4600, 10000] + # y2_plot = [5000, 5000, 25000, 25000, 5000, 5000] + + + + # # to comment if not using sim!! + # swp2_lines[0] = x2_plot + # swp2_lines[1] = y2_plot + ## END OF SIM ONLY + + p2, = ax2.plot(x2_plot, y2_plot, linestyle="-", alpha=0.75, lw=2, label = title_swp2, color="black") lines_fig2.append(p2) # Plotting (fig 3) which is the same but log scale for x - p3, = ax3.plot(x2_plot, y2_plot, linestyle="-", alpha=0.75, lw=2, label = 'swp2', color="black") + p3, = ax3.plot(x2_plot, y2_plot, linestyle="-", alpha=0.75, lw=2, label = title_swp2, color="black") lines_fig3.append(p3) min_x = 1 min_y = 1 @@ -535,6 +579,14 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax = x2_plot, y2_plot = plot_straight_x_y(x,y) if subset is not None: if breaks in subset: + with open(title+"_theo.csv", "w") as filin: + filin.write("xtheo,ytheo\n") + for i in range(len(swp2_x)): + filin.write(str(swp2_x[i]) + "," + str(swp2_y[i]) + "\n") + with open(title+"_obs.csv", "w") as filin: + filin.write("xobs,yobs\n") + for i in range(len(x)): + filin.write(str(x2_plot[i]) + "," + str(y2_plot[i]) + "\n") masking_alpha = 0.75 autoscale = True min_x = min(min_x, min(x2_plot)) @@ -554,20 +606,31 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax = Ne_max_below_limit = y[min(x.index(t_max_below_limit)+1, len(y)-1)] Ne_min_below_limit = y[x.index(t_min_below_limit)] if recent_change: - print(f"\n{breaks} breaks ; This is below the recent limit of {recent_limit_years} years:\n", - f"t_min (most recent time point under the limit) : {t_min_below_limit/mu*tgen:.1f} t_max (most ancient time point under the limit) : {t_max_below_limit/mu*tgen:.1f}", - f"\nNe_min (effective size at t_min) : {Ne_min_below_limit/(4*mu):.1f} Ne_max (effective size at t_max) : {Ne_max_below_limit/(4*mu):.1f}", - f"\nNe_min/Ne_max = {(Ne_min_below_limit/(4*mu)) / (Ne_max_below_limit/(4*mu)):.1f}", - f"\nEvolution: {((Ne_min_below_limit/(4*mu)) - (Ne_max_below_limit/(4*mu)))/((Ne_max_below_limit/(4*mu)))*100:.1f}%") - else: - print(f"Recent event under {recent_limit_years} years: NA") - # need to compute the last change and when it occured - tmin = x[1] - tmin_plus_1 = x[2] - Ne_min = y[1] - Ne_min_plus_1 = y[2] - print(f"Last was {tmin/mu*tgen:.1f} years ago. And was of {((Ne_min/(4*mu)) - (Ne_min_plus_1/(4*mu)))/(Ne_min_plus_1/(4*mu))*100:.1f}%") + # print(t_max_below_limit, y[min(x.index(t_max_below_limit)+1, len(y)-1)]) + recent_events_x = x[0:min(x.index(t_max_below_limit)+1, len(y)-1)+1] + recent_events_y = y[0:min(x.index(t_max_below_limit)+1, len(y)-1)+1] + size_changes = [] + for i in range(len(recent_events_y)-1): + size = recent_events_y[i] + # size of t+1 + size_1 = recent_events_y[i+1] + size_changes.append(((size) - (size_1))/((size_1))*100) + print(f"\n{breaks} breaks ; This is below the recent limit of {recent_limit_years} years:\n", + f"t_min (most recent time point under the limit) : {t_min_below_limit*tgen:.1f} t_max (most ancient time point under the limit) : {t_max_below_limit*tgen:.1f}", + f"\nNe_min (effective size at t_min) : {Ne_min_below_limit:.1f} Ne_max (effective size at t_max) : {Ne_max_below_limit:.1f}", + f"\nNe_min/Ne_max = {(Ne_min_below_limit) / (Ne_max_below_limit):.1f}", + f"\nEvolution: {((Ne_min_below_limit) - (Ne_max_below_limit))/((Ne_max_below_limit))*100:.1f}%", + f"\nDetail: {len(size_changes)} events since the last {recent_limit_years} years (each event size change in %) : {[round(k,2) for k in size_changes]}") + + else: + print(f"{breaks} breaks ; Last recent event under {recent_limit_years} years: NA") + # get the last change and when it occured in all cases (under the recent limit or not) + tmin = x[1] + tmin_plus_1 = x[2] + Ne_min = y[1] + Ne_min_plus_1 = y[2] + print(f"{breaks} breaks ; Last event was {tmin*tgen:.1f} years ago. And size change was of {((Ne_min) - (Ne_min_plus_1))/(Ne_min_plus_1)*100:.1f}%") else: masking_alpha = 0 autoscale = False @@ -594,12 +657,15 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax = if ax is None: # if not ax, then use the plt syntax, not ax... plt.xlabel(xlabel, fontsize=fnt_size) - plt.ylabel(ylabel, fontsize=fnt_size) + plt.ylabel(ylabel, fontsize=fnt_size*0.8) plt.gca().set_xlim(0, recent_limit * 3) if recent_change: plt.ylim(Ne_min_below_limit/3, Ne_max_below_limit *3) else: - plt.ylim(y2_plot[0]/3, y2_plot[0]) + plt.ylim(y[0]/3, y[0]+y[0]*0.5) + plt.gca().set_xlim(0, x[0] * 3) + + # plt.ylim(0, max(max_y+(max_y*0.05), max(swp2_lines[1])+(max(swp2_lines[1])*0.05))) #plt.xlim(0, recent_limit * 3) #xlim_val = plt.gca().get_xlim() @@ -607,7 +673,7 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax = # plt.xlim(min(min_x,min(swp2_lines[0])), max(max(swp2_lines[0]), max_x)) # x_ticks = list(plt.gca().get_xticks()) # plt.gca().set_xticks(x_ticks) - # plt.xticks(x_ticks) + plt.xticks(x_ticks) # plt.gca().set_xlim(xlim_val) # plt.gca().set_xticklabels([f'{k:.0e}\n{k/(mu):.0e}\n{k/(mu)*tgen:.0e}' for k in x_ticks], fontsize = fnt_size*0.5) plt.gca().set_xticklabels([f'{k:.1f}\n{k*tgen:.1f}' for k in x_ticks], fontsize = fnt_size*0.5) @@ -623,20 +689,21 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax = # plt.title(title, fontsize=fnt_size) # plt.legend(handles=lines_fig2, loc='best', fontsize = fnt_size*0.5) # # plt.text(-0.13, -0.135, 'Coal. time\nGen. time\nYears', ha='left', va='bottom', transform=ax3.transAxes) - plt.text(-0.13, -0.135, 'Gen. time\nYears', ha='left', va='bottom', transform=ax3.transAxes) + plt.text(-0.15, -0.155, 'Gen.\nYears', ha='left', va='bottom', transform=ax3.transAxes) plt.subplots_adjust(bottom=0.2) # Adjust the value as needed + plt.tight_layout() plt.savefig(title+'_plotB_'+str(nb_epochs)+'_epochs.pdf') # close fig2 to save memory plt.close(fig2) else: # when ax subplotting is used ax2.set_xlabel(xlabel, fontsize=fnt_size) - ax2.set_ylabel(ylabel, fontsize=fnt_size) + ax2.set_ylabel(ylabel, fontsize=fnt_size*0.8) ax2.set_title(title, fontsize=fnt_size) ax2.legend(handles=lines_fig2, loc='best', fontsize = fnt_size*0.5) ax3.set_xlabel(xlabel, fontsize=fnt_size) - ax3.set_ylabel(ylabel, fontsize=fnt_size) + ax3.set_ylabel(ylabel, fontsize=fnt_size*0.8) ax3.set_title(title, fontsize=fnt_size) ax3.legend(handles=lines_fig3, loc='best', fontsize = fnt_size*0.5) ax3.set_xscale('log') @@ -644,10 +711,13 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax = # Scale the x-axis # x_ticks = list(ax3.get_xticks()) # ax3.set_xticks(x_ticks) - # x_ticks = [i for i in range(0.1,max(max_x, max(swp2_lines[0]))), ] + # x_ticks = [i for i in range(0.1,max(max_x, max(swp2_lines[0]))) ] + # ax3.set_xticks(x_ticks) - ax3.set_xlim(0.1, max(max_x, max(swp2_lines[0]))) + ax3.set_xlim(1, max(max_x, max(swp2_lines[0]))) x_ticks = ax3.get_xticks() + ax3.set_xticks(x_ticks) + ax3.set_xlim(1, max(max_x, max(swp2_lines[0]))) # ax3.set_xlim(min(min(x_ticks), min(swp2_lines[0])), max(max_x, max(swp2_lines[0]))) # ax3.set_xlim(1, max(max_x, max(swp2_lines[0]))) # ax3.set_xticklabels([f'{k:.0e}\n{k/(mu):.0e}\n{k/(mu)*tgen:.0e}' for k in x_ticks], fontsize = fnt_size*0.5) @@ -664,11 +734,14 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax = # ax3.set_yticklabels([f'{k/(4*mu):.0e}' for k in y_ticks], fontsize = fnt_size*0.5) # plt.text(-0.13, -0.135, 'Coal. time\nGen. time\nYears', ha='left', va='bottom', transform=ax3.transAxes) # plt.text(-0.13, -0.135, 'Gen. time\nYears', ha='left', va='bottom', transform=ax3.transAxes) - plt.text(-0.13, -0.085, 'Gen. time\nYears', ha='left', va='bottom', transform=ax3.transAxes) + # plt.text(-0.13, -0.085, 'Gen. time\nYears', ha='left', va='bottom', transform=ax3.transAxes) plt.subplots_adjust(bottom=0.2) # Adjust the value as needed + plt.subplots_adjust(top=0.9) # Adjust the value as needed + if ax is None: # nb of plot_lines represent the number of epochs stored (len(plot_lines) = #breaks+1) + # plt.tight_layout() plt.savefig(title+'_plotC_'+str(nb_epochs)+'_epochs_log.pdf') # close fig3 to save memory plt.close(fig3) @@ -700,7 +773,9 @@ def plot_raw_stairs(plot_lines, prop, title, ax = None, n_ticks = 10, rescale = x,y = plot x_plot, y_plot = plot_straight_x_y(x,y) p, = ax1.plot(x_plot, y_plot, 'o', linestyle="-", alpha=0.75, lw=2, label = str(breaks)+' brks') - print("breaks=", breaks, "theta0", y[0]) + print("breaks=", breaks) + # for xi in range(len(x)): + # print("times", xi, "theta", y[xi]) # add plot to the list of all plots to superimpose plots.append(p) x_ticks = x @@ -767,18 +842,31 @@ def combined_plot(folder_path, mu, tgen, breaks, title = "Title", theta_scale = # Start of Parsing real swp2 output folder_splitted = folder_path.split("/") + swp2_blueprint = "/".join(folder_splitted[:-2])+".blueprint" + # parsing blueprint and get some parameters from the run: + with open(swp2_blueprint) as filin: + for line in filin: + if line.startswith("year_per_generation"): + blueprint_tgen = float(line.strip().split(":")[-1].strip()) + + #print("BLUEPRINT",swp2_blueprint, "tgen", blueprint_tgen ) swp2_summary = "/".join(folder_splitted[:-2])+'/'+folder_splitted[-3]+".final.summary" swp2_vals = parse_stairwayplot_output_summary(stwplt_out = swp2_summary) swp2_x, swp2_y = swp2_vals[0], swp2_vals[1] + if blueprint_tgen != tgen: + # need to rescale swp2 output to the new tgen + for i in range(len(swp2_x)): + swp2_x[i] = swp2_x[i] / blueprint_tgen * tgen remove_back_and_forth_points(swp2_x, swp2_y) # End of Parsing real swp2 output plot_raw_stairs(plot_lines = loaded_data['raw_stairs'], prop = loaded_data['prop'], title = title, ax = None, max_breaks = breaks) - plot_scaled_theta(plot_lines = loaded_data['scaled_stairs'], mu = mu, tgen = tgen, subset=[loaded_data['best_epoch_by_AIC']]+selected_breaks, + #plot_scaled_theta(plot_lines = loaded_data['scaled_stairs'], mu = mu, tgen = tgen, subset=[loaded_data['best_epoch_by_AIC']]+selected_breaks, + plot_scaled_theta(plot_lines = loaded_data['scaled_stairs'], mu = mu, tgen = tgen, subset=selected_breaks, # plot_scaled_theta(plot_lines = loaded_data['scaled_stairs'], subset=list(range(0,3))+[loaded_data['best_epoch_by_AIC']]+selected_breaks, prop = loaded_data['prop'], title = title, swp2_lines = [swp2_x, swp2_y], ax = None) plot_all_epochs_thetafolder(loaded_data, mu, tgen, title, theta_scale, ax = None) - + plt.close() # plt.close(fig1) # plt.close(fig2)