From 2253ce5f6775c9ad59fdd926a4096042c3965eee Mon Sep 17 00:00:00 2001 From: tforest Date: Tue, 28 Nov 2023 19:36:49 +0100 Subject: [PATCH] theoritical logLn swp2 --- swp2.py | 150 ++++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 118 insertions(+), 32 deletions(-) diff --git a/swp2.py b/swp2.py index 2f7938d..c321d00 100644 --- a/swp2.py +++ b/swp2.py @@ -1,12 +1,38 @@ import matplotlib.pyplot as plt import os import numpy as np +import math +from scipy.special import gammaln +from matplotlib.backends.backend_pdf import PdfPages + +def log_facto(k): + k = int(k) + if k > 1e6: + return k * np.log(k) - k + np.log(2*math.pi*k)/2 + val = 0 + for i in range(2, k+1): + val += np.log(i) + return val + +def log_facto_1(k): + startf = 1 # start of factorial sequence + stopf = int(k+1) # end of of factorial sequence + + q = gammaln(range(startf+1, stopf+1)) # n! = G(n+1) + + return q[-1] def return_x_y_from_stwp_theta_file(stwp_theta_file, breaks, mu, tgen, relative_theta_scale = False): with open(stwp_theta_file, "r") as swp_file: # Read the first line line = swp_file.readline() L = float(line.split()[2]) + rands = swp_file.readline() + line = swp_file.readline() + # skip empty lines before SFS + while line == "\n": + line = swp_file.readline() + sfs = np.array(line.split()).astype(float) # Process lines until the end of the file while line: # check at each line @@ -27,7 +53,7 @@ def return_x_y_from_stwp_theta_file(stwp_theta_file, breaks, mu, tgen, relative_ #### END of parsing # quit this file if the number of dimensions is incorrect if dim < breaks+1: - return 0,0,0 + return 0,0,0,0,0 # get n, the last bin of the last group # revert the list of groups as the most recent times correspond # to the closest and last leafs of the coal. tree. @@ -90,13 +116,19 @@ def return_x_y_from_stwp_theta_file(stwp_theta_file, breaks, mu, tgen, relative_ # # # divide by N0 # # y[i] = y[i]/N0 # # x[i] = x[i]/N0 - return x,y,likelihood + return x,y,likelihood,sfs,L def return_x_y_from_stwp_theta_file_as_is(stwp_theta_file, breaks, mu, tgen, relative_theta_scale = False): with open(stwp_theta_file, "r") as swp_file: # Read the first line line = swp_file.readline() L = float(line.split()[2]) + rands = swp_file.readline() + line = swp_file.readline() + # skip empty lines before SFS + while line == "\n": + line = swp_file.readline() + sfs = np.array(line.split()).astype(float) # Process lines until the end of the file while line: # check at each line @@ -117,7 +149,7 @@ def return_x_y_from_stwp_theta_file_as_is(stwp_theta_file, breaks, mu, tgen, rel #### END of parsing # quit this file if the number of dimensions is incorrect if dim < breaks+1: - return 0,0,0 + return 0,0 # get n, the last bin of the last group # revert the list of groups as the most recent times correspond # to the closest and last leafs of the coal. tree. @@ -130,7 +162,7 @@ def return_x_y_from_stwp_theta_file_as_is(stwp_theta_file, breaks, mu, tgen, rel groups[i] = groups[i].split(',') #print(groups[i], len(groups[i])) thetas[i] = [float(theta_site[i]), groups[i], likelihood] - return thetas + return thetas, sfs def plot_k_epochs_thetafolder(folder_path, mu, tgen, breaks = 2, title = "Title", theta_scale = True): scenari = {} @@ -138,9 +170,9 @@ def plot_k_epochs_thetafolder(folder_path, mu, tgen, breaks = 2, title = "Title" for file_name in os.listdir(folder_path): if os.path.isfile(os.path.join(folder_path, file_name)): # Perform actions on each file - x,y,likelihood = return_x_y_from_stwp_theta_file(folder_path+file_name, breaks = breaks, + x,y,likelihood,sfs,L = return_x_y_from_stwp_theta_file(folder_path+file_name, breaks = breaks, tgen = tgen, - mu = mu, relative_theta_scale = theta_scale) + mu = mu, relative_theta_scale = theta_scale) if x == 0 or y == 0: continue cpt +=1 @@ -189,15 +221,17 @@ def plot_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc breaks = 0 cpt +=1 if os.path.isfile(os.path.join(folder_path, file_name)): - x, y, likelihood = return_x_y_from_stwp_theta_file(folder_path+file_name, breaks = breaks, + x, y, likelihood, sfs, L = return_x_y_from_stwp_theta_file(folder_path+file_name, breaks = breaks, tgen = tgen, mu = mu, relative_theta_scale = theta_scale) + SFS_stored = sfs + L_stored = L while not (x == 0 and y == 0): if breaks not in epochs.keys(): epochs[breaks] = {} epochs[breaks][likelihood] = x,y breaks += 1 - x,y,likelihood = return_x_y_from_stwp_theta_file(folder_path+file_name, breaks = breaks, + x,y,likelihood,sfs,L = return_x_y_from_stwp_theta_file(folder_path+file_name, breaks = breaks, tgen = tgen, mu = mu, relative_theta_scale = theta_scale) print("\n*******\n"+title+"\n--------\n"+"mu="+str(mu)+"\ntgen="+str(tgen)+"\nbreaks="+str(breaks)+"\n*******\n") @@ -208,7 +242,7 @@ def plot_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc plt.figure(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) plt.xlim(1e-3, 1) #plt.ylim(0, 10) - #plt.yscale('log') + plt.yscale('log') plt.xscale('log') plt.grid(True,which="both", linestyle='--', alpha = 0.3) brkpt_lik = [] @@ -228,7 +262,15 @@ def plot_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc # divide by N0 y[i] = y[i]/N0 x[i] = x[i]/N0 - plt.plot(x, y, '-', alpha=0.75, lw=2, label = str(epoch)+' BrkPt | Lik='+greatest_likelihood) + sum_theta_i = 0 + print(epoch, x, y) + for i in range(2, len(y)-1): + sum_theta_i=y[i] / (i-1) + prop = [] + for k in range(2, len(y)-1): + prop.append(y[k+1] / (k - 1) / sum_theta_i) + #print(epoch, prop) + plt.plot(x, y, 'o', linestyle = "-", alpha=0.75, lw=2, label = str(epoch)+' BrkPt | Lik='+greatest_likelihood) if theta_scale: plt.xlabel("Coal. time") plt.ylabel("Pop. size scaled by N0") @@ -246,35 +288,67 @@ def plot_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc plt.title(title) plt.savefig(title+'_b'+str(breaks)+'.pdf') # plot likelihood against nb of breakpoints + + # best possible likelihood from SFS + # Segregating sites + S = sum(SFS_stored) + # number of monomorphic sites + L = L_stored + S0 = L-S + print("SFS", SFS_stored) + print("S", S, "L", L, "S0=", S0) + # 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)): + p_i = SFS_stored[xi] / float(S+S0) + Ln += np.log(p_i) * SFS_stored[xi] - log_facto(SFS_stored[xi]) + res = Ln + print(res) + # basic plot likelihood plt.figure(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) plt.rcParams['font.size'] = '18' - AIC = 2*(len(brkpt_lik)+1)+2*np.array(brkpt_lik)[:, 1].astype(float) - plt.plot(np.array(brkpt_lik)[:, 0], AIC, 'o', linestyle = "dotted", lw=2) - plt.axhline(y=106) + plt.plot(np.array(brkpt_lik)[:, 0], np.array(brkpt_lik)[:, 1].astype(float), 'o', linestyle = "dotted", lw=2) + # plt.ylim(0,100) + # plt.axhline(y=res) plt.yscale('log') plt.xlabel("# breakpoints", fontsize=20) plt.ylabel("$-\log\mathcal{L}$") #plt.legend(loc='upper right') plt.title(title) plt.savefig(title+'_Breakpts_Likelihood.pdf') + # AIC + plt.figure(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) + plt.rcParams['font.size'] = '18' + AIC = 2*(len(brkpt_lik)+1)+2*np.array(brkpt_lik)[:, 1].astype(float) + plt.plot(np.array(brkpt_lik)[:, 0], AIC, 'o', linestyle = "dotted", lw=2) + # plt.axhline(y=106) + plt.yscale('log') + plt.xlabel("# breakpoints", fontsize=20) + plt.ylabel("AIC") + #plt.legend(loc='upper right') + plt.title(title) + plt.savefig(title+'_Breakpts_Likelihood_AIC.pdf') -def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, breaks_max = 6): +def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, breaks_max = 5): + """ + Use theta values as is to do basic plots. + """ cpt = 0 epochs = {} for file_name in os.listdir(folder_path): cpt +=1 if os.path.isfile(os.path.join(folder_path, file_name)): for k in range(breaks_max): - thetas = return_x_y_from_stwp_theta_file_as_is(folder_path+file_name, breaks = k, + thetas,sfs = return_x_y_from_stwp_theta_file_as_is(folder_path+file_name, breaks = k, tgen = tgen, mu = mu, relative_theta_scale = theta_scale) - if thetas[0] == 0: + if thetas == 0: continue epochs[k] = thetas print("\n*******\n"+title+"\n--------\n"+"mu="+str(mu)+"\ntgen="+str(tgen)+"\nbreaks="+str(k)+"\n*******\n") print(cpt, "theta file(s) have been scanned.") - # intialize figure + # intialize figure 1 my_dpi = 300 plt.figure(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) for epoch, theta in epochs.items(): @@ -294,7 +368,7 @@ def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, plt.ylabel("theta") plt.legend(loc='upper right') plt.savefig(title+'_test'+str(k)+'.pdf') - # fig 2 + # fig 2 & 3 plt.figure(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) for epoch, theta in epochs.items(): groups = np.array(list(theta.values()), dtype=object)[:, 1].tolist() @@ -308,29 +382,41 @@ def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, N0 = y[0] for i in range(len(y)): y[i] = y[i]/N0 - # x_2 = [] T = 0 - # k allant de de 14 à 2 for i in range(len(x)): x[i] = int(x[i]) - #print(x[2]) + # compute the times as: theta_k / (k*(k-1)) for i in range(0, len(x)): - k = x[i] - #print(k, y[k-2]) - - #theta_k = y[k] / (k*(k-1)) T += y[i] / (x[i]*(x[i]-1)) x_2.append(T) + # Plotting (fig 2) + plt.plot(x_2, y, 'o', linestyle="dotted", alpha=0.75, lw=2, label = str(epoch)+' brks') + plt.xlabel("# breaks") + plt.ylabel("theta") + plt.legend(loc='upper right') + plt.savefig(title+'_test'+str(k)+'.pdf') - plt.plot(x_2, y, 'o', linestyle="dotted", alpha=0.75, lw=2, label = str(epoch)+' brks') - plt.xscale('log') - plt.xlabel("# breaks") - plt.ylabel("theta") - plt.legend(loc='upper right') - plt.savefig(title+'_test'+str(k)+'.pdf') - # + # Plotting (fig 3) which is the same but log scale for x + plt.plot(x_2, y, 'o', linestyle="dotted", alpha=0.75, lw=2, label = str(epoch)+' brks') + plt.xscale('log') + plt.xlabel("# breaks") + plt.ylabel("theta") + plt.legend(loc='upper right') + plt.savefig(title+'_test'+str(k)+'_log.pdf') +def save_multi_image(filename): + pp = PdfPages(filename) + fig_nums = plt.get_fignums() + figs = [plt.figure(n) for n in fig_nums] + for fig in figs: + fig.savefig(pp, format='pdf') + pp.close() + +def combined_plot(folder_path, mu, tgen, breaks, title = "Title", theta_scale = True): + plot_all_epochs_thetafolder(folder_path, mu, tgen, title, theta_scale) + plot_test_theta(folder_path, mu, tgen, title, theta_scale, breaks_max = breaks) + save_multi_image(title+"_combined.pdf") if __name__ == "__main__": if len(sys.argv) != 4: