-
Notifications
You must be signed in to change notification settings - Fork 0
/
TestingNumericalSolns.py
149 lines (119 loc) · 4.25 KB
/
TestingNumericalSolns.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
import scipy.optimize as op
import DeterministicPrediction as det
mutdat32_fromscript = np.loadtxt("mutdat32_fromscript.csv", delimiter = ",")
mm32 = np.transpose(mutdat32_fromscript)
#siloc = input("Location of set of selection coefficients to solve with: ")
#si = np.loadtxt(siloc, delimiter = ",")
for ii in range(mutdat32_fromscript.shape[0]):
mm32[ii,ii]=-(np.sum(mm32[ii])-mm32[ii,ii])
def g(xv) :
x=31;
mut=np.zeros((x,x))
for ii in range(x):
for jj in range(x):
if ii!=jj:
mut[ii,jj]=-xv[ii]*xv[jj]
else:
mut[ii,jj]=xv[ii]*(1.0 - xv[ii])
return mut
def losseq(xg, sv):
#x = np.append(xg, 1.0-np.sum(xg))
vec_eq = np.zeros(xg.shape[0])
for i in range(vec_eq.shape[0]):
vec_eq[i] = np.sum(mm32[i, :]*xg)+xg[i]*sv[i]
sum_xg_isv_i = np.sum(xg*sv)
vec_eq[i] = vec_eq[i] - xg[i]*sum_xg_isv_i
return vec_eq
def losseq2(xg, sv):
'''
Original function from Shamreen
'''
sM=np.diag(sv)
DM=np.identity(32)
x = np.append(xg, 1.0-np.sum(xg))
sv31 = sv[0:31] # Leaving out last element for dimensional compatibility
T1 = np.dot(np.dot(mm32,(DM+sM)),x)[0:31]
T2 = np.dot(g(xg),sv31)
return (T1 + T2)
def catch_negs(Xt):
Xt_min = np.min(Xt)
if Xt_min < 0:
return True
else:
return False
def solveSequence(xeq0, si):
xg = xeq0
Xeq = []
for i in tqdm(range(si.shape[0])): #tqdm here
#if i <= 4/9*si.shape[0]:
Xt = op.fsolve(losseq, xg, args = si[i])
#Xt = np.append(Xt, 1-np.sum(Xt))
#diffs = np.abs(Xt - xg)
# if np.argmax(Xt) != np.argmax(xg):
# print(i)
# X_freqs = det.path(xg, si[i], 1000, track = True)
# Xt = op.fsolve(losseq, X_freqs[-1, :], args = si[i])
while catch_negs(Xt) == True: #or np.abs(np.sum(Xt - xg)) > 1e-15:
#print("Using approx")
Xsol = Xt
#Xt = np.append(xg, 1-np.sum(xg)) ##NOT KOSHER
X_freqs = det.path(Xsol, si[i], 10000, track = False)
Xsol = op.fsolve(losseq, X_freqs[-1, :], args = si[i])
xg = X_freqs[-1, :]
Xt = Xsol
Xeq.append(Xt)
xg = Xt
# print(i)
#else:
# Xeq.append(Xt)
Xeq = np.array(Xeq)
return Xeq
def solveSequenceOldLosseq(xeq0, si):
xg = xeq0[:-1]
Xeq = []
for i in tqdm(range(si.shape[0])): #tqdm here
#if i <= 4/9*si.shape[0]:
Xt = op.fsolve(losseq2, xg, args = si[i])
Xt = np.append(Xt, 1-np.sum(Xt))
#diffs = np.abs(Xt - xg)
# if np.argmax(Xt) != np.argmax(xg):
# print(i)
# X_freqs = det.path(xg, si[i], 1000, track = True)
# Xt = op.fsolve(losseq, X_freqs[-1, :], args = si[i])
while catch_negs(Xt) == True: #or np.abs(np.sum(Xt - xg)) > 1e-15:
Xsol = Xt
#print("Using approx")
#Xt = np.append(xg, 1-np.sum(xg)) ##NOT KOSHER
X_freqs = det.path(Xsol, si[i], 10000, track = False)
Xsol = op.fsolve(losseq2, X_freqs[-1, :-1], args = si[i]) #xt has length 31
Xsol = np.append(Xsol, 1-np.sum(Xsol))
xg = X_freqs[-1, :]
Xt = Xsol
Xeq.append(Xt)
xg = Xt[:-1]
# print(i)
#else:
# Xeq.append(Xt)
Xeq = np.array(Xeq)
return Xeq
def solveSequenceNoApprox(xeq0, si):
xg = xeq0
Xeq = []
for i in tqdm(range(si.shape[0])):
if i <= 4/9*si.shape[0]:
Xt = op.fsolve(losseq, xg, args = si[i])
#Xt = np.append(Xt, 1-np.sum(Xt))
# if catch_negs(Xt) == True:
# print("Using approx")
# #Xt = np.append(xg, 1-np.sum(xg)) ##NOT KOSHER
# X_freqs = det.path(Xeq[-1], si[i], 100000, track = False)
# Xt = op.fsolve(losseq, X_freqs[-1, :], args = si[i])
Xeq.append(Xt)
xg = Xt
else:
Xeq.append(Xt)
Xeq = np.array(Xeq)
return Xeq