From 6a6d4bf6f94a970e15df700c959a2afdaf094197 Mon Sep 17 00:00:00 2001 From: tforest Date: Fri, 8 Dec 2023 22:09:53 +0100 Subject: [PATCH] Remove unused functions for swp2 --- swp2.py | 71 +++++++-------------------------------------------------- 1 file changed, 8 insertions(+), 63 deletions(-) diff --git a/swp2.py b/swp2.py index 6c4ac18..499c481 100644 --- a/swp2.py +++ b/swp2.py @@ -3,13 +3,11 @@ import os import numpy as np import math import json -import io -from scipy.special import gammaln -from matplotlib.backends.backend_pdf import PdfPages -from matplotlib.ticker import MaxNLocator -from mpl_toolkits.axes_grid1.inset_locator import inset_axes -from matplotlib.ticker import MultipleLocator + def log_facto(k): + """ + Using the Stirling's approximation + """ k = int(k) if k > 1e6: return k * np.log(k) - k + np.log(2*math.pi*k)/2 @@ -18,14 +16,6 @@ def log_facto(k): 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 @@ -128,51 +118,6 @@ def return_x_y_from_stwp_theta_file(stwp_theta_file, breaks, mu, tgen, relative_ # # x[i] = x[i]/N0 return x,y,likelihood,thetas,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 - if line.startswith("dim") : - dim = int(line.split()[1]) - if dim == breaks+1: - likelihood = line.split()[5] - groups = line.split()[6:6+dim] - theta_site = line.split()[6+dim:6+dim+1+dim] - elif dim < breaks+1: - line = swp_file.readline() - continue - elif dim > breaks+1: - break - #return 0,0,0 - # Read the next line - line = swp_file.readline() - #### END of parsing - # quit this file if the number of dimensions is incorrect - if dim < breaks+1: - 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. - groups = groups[::-1] - theta_site = theta_site[::-1] - - thetas = {} - - for i in range(len(groups)): - groups[i] = groups[i].split(',') - # print(groups[i], len(groups[i])) - thetas[i] = [float(theta_site[i]), groups[i], likelihood] - return thetas, sfs def plot_k_epochs_thetafolder(folder_path, mu, tgen, breaks = 2, title = "Title", theta_scale = True): scenari = {} @@ -180,7 +125,7 @@ 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,sfs,L = return_x_y_from_stwp_theta_file(folder_path+file_name, breaks = breaks, + x, y, likelihood, theta, sfs, L = return_x_y_from_stwp_theta_file(folder_path+file_name, breaks = breaks, tgen = tgen, mu = mu, relative_theta_scale = theta_scale) if x == 0 or y == 0: @@ -396,7 +341,7 @@ def save_k_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, cpt +=1 if os.path.isfile(os.path.join(folder_path, file_name)): for k in range(breaks_max): - thetas,sfs = return_x_y_from_stwp_theta_file_as_is(folder_path+file_name, breaks = k, + x,y,likelihood,thetas,sfs,L = return_x_y_from_stwp_theta_file(folder_path+file_name, breaks = k, tgen = tgen, mu = mu, relative_theta_scale = theta_scale) if thetas == 0: @@ -591,7 +536,7 @@ def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, cpt +=1 if os.path.isfile(os.path.join(folder_path, file_name)): for k in range(breaks_max): - thetas,sfs = return_x_y_from_stwp_theta_file_as_is(folder_path+file_name, breaks = k, + x, y, likelihood, theta, sfs, L = return_x_y_from_stwp_theta_file(folder_path+file_name, breaks = k, tgen = tgen, mu = mu, relative_theta_scale = theta_scale) if thetas == 0: @@ -751,7 +696,7 @@ def combined_plot(folder_path, mu, tgen, breaks, title = "Title", theta_scale = # # 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() - # save_k_theta(folder_path, mu, tgen, title, theta_scale, breaks_max = breaks, output = title+"_plotdata.json") + save_k_theta(folder_path, mu, tgen, title, theta_scale, breaks_max = breaks, output = title+"_plotdata.json") with open(title+"_plotdata.json", 'r') as json_file: loaded_data = json.load(json_file) # plot page 1 of summary