richard
/
geinf
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Remove unused functions for swp2

master
tforest 2023-12-08 22:09:53 +01:00
parent 89813468b5
commit 6a6d4bf6f9
1 changed files with 8 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

71
swp2.py
View File

@ -3,13 +3,11 @@ import os
import numpy as np import numpy as np
import math import math
import json 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): def log_facto(k):
"""
Using the Stirling's approximation
"""
k = int(k) k = int(k)
if k > 1e6: if k > 1e6:
return k * np.log(k) - k + np.log(2*math.pi*k)/2 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) val += np.log(i)
return val 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): 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: with open(stwp_theta_file, "r") as swp_file:
# Read the first line # 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 # # x[i] = x[i]/N0
return x,y,likelihood,thetas,sfs,L 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): def plot_k_epochs_thetafolder(folder_path, mu, tgen, breaks = 2, title = "Title", theta_scale = True):
scenari = {} 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): for file_name in os.listdir(folder_path):
if os.path.isfile(os.path.join(folder_path, file_name)): if os.path.isfile(os.path.join(folder_path, file_name)):
# Perform actions on each file # 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, tgen = tgen,
mu = mu, relative_theta_scale = theta_scale) mu = mu, relative_theta_scale = theta_scale)
if x == 0 or y == 0: 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 cpt +=1
if os.path.isfile(os.path.join(folder_path, file_name)): if os.path.isfile(os.path.join(folder_path, file_name)):
for k in range(breaks_max): 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, tgen = tgen,
mu = mu, relative_theta_scale = theta_scale) mu = mu, relative_theta_scale = theta_scale)
if thetas == 0: if thetas == 0:
@ -591,7 +536,7 @@ def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True,
cpt +=1 cpt +=1
if os.path.isfile(os.path.join(folder_path, file_name)): if os.path.isfile(os.path.join(folder_path, file_name)):
for k in range(breaks_max): 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, tgen = tgen,
mu = mu, relative_theta_scale = theta_scale) mu = mu, relative_theta_scale = theta_scale)
if thetas == 0: 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_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) # # plot_test_theta(folder_path, mu, tgen, title, theta_scale, breaks_max = breaks, ax = None)
# # plt.clf() # # 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: with open(title+"_plotdata.json", 'r') as json_file:
loaded_data = json.load(json_file) loaded_data = json.load(json_file)
# plot page 1 of summary # plot page 1 of summary