In [2]:
from pyCodeLib import *
import warnings
import glob
import re
import numpy as np
import pandas as pd
from Bio.PDB.Polypeptide import one_to_three
warnings.filterwarnings('ignore')
# sys.path.insert(0, MYHOME)
%load_ext autoreload
%autoreload 2
In [3]:
# plt.rcParams['figure.figsize'] = [16.18033, 10]
plt.rcParams['figure.figsize'] = 0.5*np.array([16.18033, 10])
# def gamma_format_convertion_iteration_to_simulation(iteration_gamma, gamma_for_simulation):
# # gamma_location = "/Users/weilu/Research/server_backup/jan_2019/optimization/gammas_dec30/cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_gamma"
# # gamma_for_simulation = "/Users/weilu/Research/server_backup/jan_2019/optimization/iteration_gamma.dat"
# gamma = iteration_gamma
# gamma = -gamma # caused by tradition.
# # convert gamma to gamma used by simulation
# with open(gamma_for_simulation, "w") as out:
# c = 0
# for i in range(20):
# for j in range(i, 20):
# out.write(f"{gamma[c]:<.5f} {gamma[c]:10.5f}\n")
# c += 1
# out.write("\n")
# for i in range(20):
# for j in range(i, 20):
# # protein, water
# out.write(f"{gamma[c]:<.5f} {gamma[c+210]:10.5f}\n")
# c += 1
res_type_map = {
'A': 0,
'C': 4,
'D': 3,
'E': 6,
'F': 13,
'G': 7,
'H': 8,
'I': 9,
'K': 11,
'L': 10,
'M': 12,
'N': 2,
'P': 14,
'Q': 5,
'R': 1,
'S': 15,
'T': 16,
'V': 19,
'W': 17,
'Y': 18
}
# res_type_map = gamma_se_map_1_letter = { 'A': 0, 'R': 1, 'N': 2, 'D': 3, 'C': 4,
# 'Q': 5, 'E': 6, 'G': 7, 'H': 8, 'I': 9,
# 'L': 10, 'K': 11, 'M': 12, 'F': 13, 'P': 14,
# 'S': 15, 'T': 16, 'W': 17, 'Y': 18, 'V': 19}
res_type_map_letters = ['A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G',
'H', 'I', 'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V']
inverse_res_type_map = dict(list(zip(list(range(20)), res_type_map_letters)))
def gamma_format_convertion_iteration_to_simulation(iteration_gamma, gamma_for_simulation, burial_gamma_for_simulation=None):
# gamma_location = "/Users/weilu/Research/server_backup/jan_2019/optimization/gammas_dec30/cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_gamma"
# gamma_for_simulation = "/Users/weilu/Research/server_backup/jan_2019/optimization/iteration_gamma.dat"
gamma = iteration_gamma
gamma = -gamma # caused by tradition.
# convert gamma to gamma used by simulation
with open(gamma_for_simulation, "w") as out:
c = 0
for i in range(20):
for j in range(i, 20):
out.write(f"{gamma[c]:<.5f} {gamma[c]:10.5f}\n")
c += 1
out.write("\n")
for i in range(20):
for j in range(i, 20):
# protein, water
out.write(f"{gamma[c]:<.5f} {gamma[c+210]:10.5f}\n")
c += 1
if burial_gamma_for_simulation:
rhoGamma = pd.DataFrame(gamma[630:].reshape(3,20).T, columns=["rho1", "rho2", "rho3"]).reset_index()
rhoGamma["oneLetter"] = rhoGamma["index"].apply(lambda x: inverse_res_type_map[x])
rhoGamma["Residue"] = rhoGamma["index"].apply(lambda x: one_to_three(inverse_res_type_map[x]))
rhoGamma = rhoGamma[["Residue", "rho1", "rho2", "rho3", "index", "oneLetter"]]
g = rhoGamma[["rho1", "rho2", "rho3"]].values
np.savetxt(burial_gamma_for_simulation, g, fmt='%7.4f')
In [73]:
# pre = "/Users/weilu/Research/server_backup/feb_2019/jan_optimization/gammas/"
pre = "/Users/weilu/Research/server/april_2019/optimization_test/gammas/"
# pp = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0"
pp = "proteins_name_list_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0phi_burial_well4.0"
A_name = pp + "_A"
B_name = pp + "_B"
B_filtered_name = pp + "_B_filtered"
P_name = pp + "_P"
Gamma_name = pp + "_gamma"
Gamma_filtered_name = pp + "_gamma_filtered"
Lamb_name = pp + "_lamb"
Lamb_filtered_name = pp + "_lamb_filtered"
A = np.loadtxt(pre+A_name)
B = np.loadtxt(pre+B_name)
B_filtered = np.loadtxt(pre+B_filtered_name, dtype=complex, converters={
0: lambda s: complex(s.decode().replace('+-', '-'))})
Gamma = np.loadtxt(pre+Gamma_name)
Gamma_filtered = np.loadtxt(pre+Gamma_filtered_name, dtype=complex, converters={
0: lambda s: complex(s.decode().replace('+-', '-'))})
Lamb = np.loadtxt(pre+Lamb_name, dtype=complex, converters={
0: lambda s: complex(s.decode().replace('+-', '-'))})
Lamb_filtered = np.loadtxt(pre+Lamb_filtered_name, dtype=complex, converters={
0: lambda s: complex(s.decode().replace('+-', '-'))})
half_B_name = pp + "_half_B"
half_B = np.loadtxt(pre+half_B_name)
other_half_B_name = pp + "_other_half_B"
other_half_B = np.loadtxt(pre+other_half_B_name)
std_half_B_name = pp + "_std_half_B"
std_half_B = np.loadtxt(pre+std_half_B_name)
pre = "/Users/weilu/Research/server/april_2019/"
location = pre + "phis/proteins_name_list_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0phi_burial_well4.0_phi_decoy_summary.txt"
A_prime = np.loadtxt(location)
In [72]:
In [34]:
A_prime.shape
Out[34]:
In [74]:
B_inverse = np.linalg.pinv(B)
In [75]:
up = A.dot(B_inverse).dot(A) - A_prime.dot(B_inverse).dot(A)
down = A.dot(B_inverse).dot(A_prime) - A_prime.dot(B_inverse).dot(A_prime)
lambda_1 = up / down
In [76]:
up
Out[76]:
In [77]:
down
Out[77]:
In [78]:
lambda_1
Out[78]:
In [79]:
g = B_inverse.dot(A - A_prime*lambda_1)
In [96]:
np.mean(g)
Out[96]:
In [100]:
np.std(g)
Out[100]:
In [97]:
g1 = B_inverse.dot(A)
In [98]:
np.mean(g1)
Out[98]:
In [99]:
np.std(g1)
Out[99]:
In [80]:
g.dot(A)
Out[80]:
In [88]:
g.dot(A_prime)
Out[88]:
In [101]:
g1.dot(A)
Out[101]:
In [102]:
g1.dot(A_prime)
Out[102]:
In [41]:
plt.plot(A)
Out[41]:
In [94]:
np.sum(abs(lambda_1*A_prime))
Out[94]:
In [95]:
np.sum(abs(A))
Out[95]:
In [111]:
plt.rcParams['figure.figsize'] = np.array([16.18033, 10])
plt.plot(g1)
plt.plot(g, alpha=0.5)
# plt.yscale("log")
Out[111]:
In [105]:
plt.plot(g1)
Out[105]:
In [103]:
plt.plot(A_prime)
plt.plot(lambda_1*A_prime)
plt.plot(A)
# plt.yscale("log")
Out[103]:
In [ ]:
In [16]:
c = A0
A = [[-3, 1], [1, 2]]
b = [6, 4]
x0_bounds = (None, None)
x1_bounds = (-3, None)
from scipy.optimize import linprog
res = linprog(c, A_ub=A, b_ub=b, bounds=(x0_bounds, x1_bounds),
options={"disp": True})
print(res)
In [ ]:
In [13]:
test = np.dot(np.linalg.pinv(B), A)
In [14]:
test[:10]
Out[14]:
In [15]:
Gamma[:10]
Out[15]:
In [10]:
total_phis = 690
num_decoys = 15000
filtered_gamma, filtered_B, filtered_lamb, P, lamb = get_filtered_gamma_B_lamb_P_and_lamb(A, B, half_B, other_half_B, std_half_B, total_phis, num_decoys, mode=2)
In [11]:
Gamma_filtered[:10]
Out[11]:
In [12]:
filtered_gamma[:10]
Out[12]:
In [80]:
Gamma[:10]
Out[80]:
In [79]:
Gamma[630:]
Out[79]:
In [101]:
burialOnly = np.loadtxt("/Users/weilu/Research/server/feb_2019/optimization_debug_only_burial/gammas/cath-dataset-nonredundant-S20Clean_phi_burial_well4.0_gamma")
In [112]:
# now, positive means favored.
rhoGamma = pd.DataFrame(-burialOnly.reshape(3,20).T, columns=["rho1", "rho2", "rho3"]).reset_index()
rhoGamma["oneLetter"] = rhoGamma["index"].apply(lambda x: inverse_res_type_map[x])
rhoGamma["Residue"] = rhoGamma["index"].apply(lambda x: one_to_three(inverse_res_type_map[x]))
rhoGamma = rhoGamma[["Residue", "rho1", "rho2", "rho3", "index", "oneLetter"]]
g = rhoGamma[["rho1", "rho2", "rho3"]].values
# np.savetxt("/Users/weilu/Research/server/feb_2019/burial_only_gamma.dat", g, fmt='%7.4f')
# rhoGamma
rhoGamma["hydrophobicityOrder"] = rhoGamma["oneLetter"].apply(lambda x: hydrophobicity_map[x])
rhoGamma.sort_values("hydrophobicityOrder")
Out[112]:
In [106]:
hydrophobicity_letters = ['R', 'K', 'N', 'Q', 'D', 'E', 'H', 'Y',
'W', 'S', 'T', 'G', 'P', 'A', 'M', 'C', 'F', 'L', 'V', 'I']
hydrophobicity_map = dict(list(zip(hydrophobicity_letters, list(range(20)))))
In [108]:
In [113]:
In [121]:
location = "/Users/weilu/Research/database/queriedPDB.dat"
with open(location) as f:
n = int(next(f))
print(n)
all_lines = f.read().splitlines()
first_20 = all_lines[:20]
print(first_20)
In [102]:
burialOnly.shape
Out[102]:
In [ ]:
Gamma = np.loadtxt(pre+Gamma_name)
In [94]:
rho_table = [[0.0, 3.0], [3.0, 6.0], [6.0, 9.0]]
for i in range(3):
print(rho_table[i][0], rho_table[i][1])
In [97]:
Gamma.shape
Out[97]:
In [89]:
rhoGamma = pd.DataFrame(Gamma[630:].reshape(3,20).T, columns=["rho1", "rho2", "rho3"]).reset_index()
rhoGamma["oneLetter"] = rhoGamma["index"].apply(lambda x: inverse_res_type_map[x])
rhoGamma["Residue"] = rhoGamma["index"].apply(lambda x: one_to_three(inverse_res_type_map[x]))
rhoGamma = rhoGamma[["Residue", "rho1", "rho2", "rho3", "index", "oneLetter"]]
g = rhoGamma[["rho1", "rho2", "rho3"]].values
# np.savetxt("/Users/weilu/Research/server/feb_2019/iter_burial_gamma.dat", g, fmt='%7.4f')
In [90]:
rhoGamma
Out[90]:
In [91]:
plot_contact_well(Gamma[:210], inferBound=True, invert_sign=False)
In [99]:
gamma_for_simulation = "/Users/weilu/Research/server/feb_2019/optimization/iteration_gamma.dat"
burial_gamma_for_simulation = "/Users/weilu/Research/server/feb_2019/optimization/iteration_burial_gamma.dat"
gamma_format_convertion_iteration_to_simulation(Gamma, gamma_for_simulation, burial_gamma_for_simulation=burial_gamma_for_simulation)
In [12]:
pre = "/Users/weilu/Research/server/feb_2019/jan_optimization/gammas/"
pp = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0"
Gamma_name = pp + "_gamma"
withoutBurialGamma = np.loadtxt(pre+Gamma_name)
In [84]:
plot_contact_well((-Gamma[:210]) - (-withoutBurialGamma[:210]), inferBound=True, invert_sign=False)
In [14]:
plot_contact_well(withoutBurialGamma[:210], inferBound=True)
In [85]:
plot_contact_well(Gamma[:210], inferBound=True, invert_sign=False)
In [67]:
plot_contact_well(Gamma_filtered[:210], inferBound=True)
In [8]:
plot_contact_well(Gamma[210:420], inferBound=True)
In [21]:
plot_contact_well(Gamma[420:630], inferBound=True)
In [19]:
plot_contact_well(withoutBurialGamma[210:420] - Gamma[210:420], inferBound=True)
In [20]:
plot_contact_well(withoutBurialGamma[420:630] - Gamma[420:630], inferBound=True)
In [35]:
# pre = "/Users/weilu/Research/server_backup/feb_2019/jan_optimization/gammas/"
pre = "/Users/weilu/Research/server_backup/feb_2019/jan_optimization/back_up_gammas/"
A_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_A"
B_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_B"
B_filtered_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_B_filtered"
P_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_P"
Gamma_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_gamma"
Gamma_filtered_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_gamma_filtered"
Lamb_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_lamb"
Lamb_filtered_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_lamb_filtered"
A = np.loadtxt(pre+A_name)
B = np.loadtxt(pre+B_name)
B_filtered = np.loadtxt(pre+B_filtered_name, dtype=complex, converters={
0: lambda s: complex(s.decode().replace('+-', '-'))})
Gamma = np.loadtxt(pre+Gamma_name)
Gamma_filtered = np.loadtxt(pre+Gamma_filtered_name, dtype=complex, converters={
0: lambda s: complex(s.decode().replace('+-', '-'))})
Lamb = np.loadtxt(pre+Lamb_name, dtype=complex, converters={
0: lambda s: complex(s.decode().replace('+-', '-'))})
Lamb_filtered = np.loadtxt(pre+Lamb_filtered_name, dtype=complex, converters={
0: lambda s: complex(s.decode().replace('+-', '-'))})
half_B_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_half_B"
half_B = np.loadtxt(pre+half_B_name)
other_half_B_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_other_half_B"
other_half_B = np.loadtxt(pre+other_half_B_name)
std_half_B_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_std_half_B"
std_half_B = np.loadtxt(pre+std_half_B_name)
In [57]:
pre = "/Users/weilu/Research/server_backup/feb_2019/jan_optimization/gammas/"
A_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_A"
B_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_B"
B_filtered_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_B_filtered"
P_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_P"
Gamma_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_gamma"
Gamma_filtered_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_gamma_filtered"
Lamb_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_lamb"
Lamb_filtered_name = "cath-dataset-nonredundant-S20Clean_phi_pairwise_contact_well4.5_6.5_5.0_10phi_density_mediated_contact_well6.5_9.5_5.0_10_2.6_7.0_lamb_filtered"
A = np.loadtxt(pre+A_name)
B = np.loadtxt(pre+B_name)
B_filtered = np.loadtxt(pre+B_filtered_name, dtype=complex, converters={
0: lambda s: complex(s.decode().replace('+-', '-'))})
Gamma = np.loadtxt(pre+Gamma_name)
Gamma_filtered = np.loadtxt(pre+Gamma_filtered_name, dtype=complex, converters={
0: lambda s: complex(s.decode().replace('+-', '-'))})
Lamb = np.loadtxt(pre+Lamb_name, dtype=complex, converters={
0: lambda s: complex(s.decode().replace('+-', '-'))})
Lamb_filtered = np.loadtxt(pre+Lamb_filtered_name, dtype=complex, converters={
0: lambda s: complex(s.decode().replace('+-', '-'))})
In [84]:
plot_contact_well(Gamma[210:420], inferBound=True)
In [63]:
gamma_for_simulation = "/Users/weilu/Research/server_backup/feb_2019/jan_optimization/iteration_gamma.dat"
burial_gamma_for_simulation = "/Users/weilu/Research/server_backup/feb_2019/jan_optimization/iteration_burial_gamma.dat"
gamma_format_convertion_iteration_to_simulation(Gamma, gamma_for_simulation, burial_gamma_for_simulation=burial_gamma_for_simulation)
In [59]:
plt.plot(Lamb)
plt.plot(Lamb_filtered)
Out[59]:
In [27]:
def get_filtered_gamma_B_lamb_P_and_lamb(A, B, half_B, other_half_B, std_half_B, total_phis, num_decoys, noise_iterations=10, relative_error_threshold=0.5):
lamb, P = np.linalg.eig(B)
lamb, P = sort_eigenvalues_and_eigenvectors(lamb, P)
cutoff_modes = []
for i_noise in range(noise_iterations):
noisy_B = np.zeros((total_phis, total_phis))
for i in range(total_phis):
for j in range(i, total_phis):
random_B_ij = np.random.normal(
loc=half_B[i][j], scale=std_half_B[i][j] / float(num_decoys))
noisy_B[i][j] = noisy_B[j][i] = random_B_ij - \
other_half_B[i][j]
noisy_lamb, noisy_P = np.linalg.eig(noisy_B)
noisy_lamb, noisy_P = sort_eigenvalues_and_eigenvectors(
noisy_lamb, noisy_P)
try:
cutoff_mode = np.where(
np.abs(lamb - noisy_lamb) / lamb > relative_error_threshold)[0][0]
except IndexError:
cutoff_mode = len(lamb)
cutoff_modes.append(cutoff_mode)
cutoff_mode = min(cutoff_modes)
print(cutoff_mode)
filtered_lamb = np.copy(lamb)
filtered_B_inv, filtered_lamb, P = get_filtered_B_inv_lambda_and_P(
filtered_lamb, cutoff_mode, P)
filtered_gamma = np.dot(filtered_B_inv, A)
filtered_B = np.linalg.inv(filtered_B_inv)
return filtered_gamma, filtered_B, filtered_lamb, P, lamb
def get_filtered_B_inv_lambda_and_P(filtered_lamb, cutoff_mode, P, method='extend_all_after_first_noisy_mode'):
if method == 'zero_all_after_first_noisy_mode':
filtered_lamb_inv = 1 / filtered_lamb
# for "zeroing unreliable eigenvalues"
filtered_lamb_inv[cutoff_mode:] = 0.0
filtered_B_inv = np.dot(
P, np.dot(np.diag(filtered_lamb_inv), np.linalg.inv(P)))
filtered_lamb = 1 / filtered_lamb_inv
if method == 'extend_all_after_first_noisy_mode':
# for "extending lowest reliable eigenvalue"
filtered_lamb[cutoff_mode:] = filtered_lamb[cutoff_mode - 1]
filtered_B_inv = np.dot(
P, np.dot(np.diag(1 / filtered_lamb), np.linalg.inv(P)))
return filtered_B_inv, filtered_lamb, P
def sort_eigenvalues_and_eigenvectors(eigenvalues, eigenvectors):
idx = eigenvalues.argsort()[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
return eigenvalues, eigenvectors
In [61]:
total_phis = 630
num_decoys = 1000
noise_iterations=10
relative_error_threshold=0.5
lamb, P = np.linalg.eig(B)
lamb, P = sort_eigenvalues_and_eigenvectors(lamb, P)
cutoff_modes = []
for i_noise in range(noise_iterations):
noisy_B = np.zeros((total_phis, total_phis))
for i in range(total_phis):
for j in range(i, total_phis):
# random_B_ij = np.random.normal(
# loc=half_B[i][j], scale=std_half_B[i][j] / float(num_decoys))
random_B_ij = np.random.normal(
loc=half_B[i][j], scale=std_half_B[i][j])
noisy_B[i][j] = noisy_B[j][i] = random_B_ij - other_half_B[i][j]
noisy_lamb, noisy_P = np.linalg.eig(noisy_B)
noisy_lamb, noisy_P = sort_eigenvalues_and_eigenvectors(
noisy_lamb, noisy_P)
try:
cutoff_mode = np.where(
np.abs(lamb - noisy_lamb) / lamb > relative_error_threshold)[0][0]
except IndexError:
cutoff_mode = len(lamb)
cutoff_modes.append(cutoff_mode)
cutoff_mode = min(cutoff_modes)
print(cutoff_mode)
filtered_lamb = np.copy(lamb)
filtered_B_inv, filtered_lamb, P = get_filtered_B_inv_lambda_and_P(
filtered_lamb, cutoff_mode, P)
filtered_gamma = np.dot(filtered_B_inv, A)
filtered_B = np.linalg.inv(filtered_B_inv)
# return filtered_gamma, filtered_B, filtered_lamb, P, lamb
In [78]:
x = np.zeros(10000)
for i in range(10):
x1 = np.random.normal(10, scale=10, size=10000)
x += x1
x = x/10
In [82]:
plt.hist(x1, bins=100)
plt.hist(x, bins=100)
Out[82]:
In [70]:
plt.imshow(noisy_B)
plt.colorbar()
Out[70]:
In [71]:
plt.imshow(B)
plt.colorbar()
Out[71]:
In [67]:
plt.plot(noisy_lamb)
Out[67]:
In [68]:
plt.plot(lamb)
Out[68]:
In [49]:
np.std([1,2,3])
Out[49]:
In [52]:
x = np.array([1,2,3])
np.sqrt(np.mean(abs(x - x.mean())**2))
Out[52]:
In [53]:
np.mean(abs(x - x.mean())**2)
Out[53]:
In [ ]:
np.std()
In [48]:
np.abs(lamb - noisy_lamb) / lamb
Out[48]:
In [47]:
np.where(np.abs(lamb - noisy_lamb) / lamb > relative_error_threshold)
Out[47]:
In [37]:
np.sum(filtered_B-B_filtered)
Out[37]:
In [45]:
plt.plot(noisy_lamb)
Out[45]:
In [43]:
plt.imshow(half_B)
plt.colorbar()
Out[43]:
In [42]:
plt.imshow(std_half_B)
plt.colorbar()
Out[42]:
In [22]:
lamb, P = np.linalg.eig(B)
lamb, P = sort_eigenvalues_and_eigenvectors(lamb, P)
In [28]:
plt.plot(lamb)
Out[28]:
In [ ]:
In [54]:
plot_contact_well(Gamma[:210], inferBound=True)
In [55]:
plot_contact_well(Gamma[210:420], inferBound=True)
In [56]:
plot_contact_well(Gamma[420:], inferBound=True)
In [5]:
plot_contact_well(Gamma_filtered[:210], inferBound=True)
In [6]:
plot_contact_well(Gamma[210:420], inferBound=True)
In [7]:
plot_contact_well(Gamma_filtered[210:420], inferBound=True)
In [8]:
plot_contact_well(Gamma[420:], inferBound=True)
In [9]:
plot_contact_well(Gamma_filtered[420:], inferBound=True)
In [10]:
np.sum(Gamma - Gamma_filtered)
Out[10]:
In [17]:
plot_contact_well(Gamma[:210], inferBound=True)
In [18]:
plot_contact_well(Gamma_filtered[:210], inferBound=True)
In [19]:
plot_contact_well(Gamma[210:420], inferBound=True)
In [20]:
plot_contact_well(Gamma_filtered[210:420], inferBound=True)
In [16]:
plot_contact_well(Gamma[420:], inferBound=True)
In [15]:
plot_contact_well(Gamma_filtered[420:], inferBound=True)
In [ ]: