-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path1d_daisy_scale.py
93 lines (79 loc) · 3.35 KB
/
1d_daisy_scale.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
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
from tqdm import tqdm
from manifold_utils.iga import iga, arccos_catch_nan
from ilc_data.ilc_loader import get_sct_var_scale, get_index_to_gene, get_radii_scale
trail_indices = np.load("data/trail_indices.npy")
scale_var = get_sct_var_scale()
indices = list(range(17))
num_inds = len(indices)
true_indices = [trail_indices[j] for j in indices]
center_list = []
radii_list = []
eigval_list = []
eigvec_list = []
for j, j_true in zip(indices, true_indices):
center_list.append(scale_var[j_true, :])
radii_list.append(np.load("scale_intermediates/radii_"+str(j)+".npy"))
eigval_list.append(np.load("scale_intermediates/eigvals_"+str(j)+".npy"))
#eigvecs = np.load("scale_intermediates/eigvecs+"+str(j)+".npy")
#eigvec_list.append(np.swapaxes(eigvecs, 1, 2))
eigvec_list.append(np.load("scale_intermediates/eigvecs+"+str(j)+".npy"))
#preRmins = [35.5, 32, 37.5, 36, 35, 37.5, 37.5, 33.5, 38, 38, 37, 36, 42, 37, 40, 38, 38]
#preRmaxs = [37.5, 33.5, 39, 37.5, 37.5, 39, 38, 34, 39.5, 39, 38.5, 39, 44, 39, 42, 40, 40]
preRmins, preRmaxs = get_radii_scale()
Rmins = [preRmins[j] for j in indices]
Rmaxs = [preRmaxs[j] for j in indices]
Rmin_inds = []
Rmax_inds = []
for radii, Rmin, Rmax in zip(radii_list, Rmins, Rmaxs):
Rmin_dists = np.square(Rmin - radii)
Rmax_dists = np.square(Rmax - radii)
Rmin_dist_min = np.min(Rmin_dists)
Rmin_ind = np.where(Rmin_dists == Rmin_dist_min)
Rmin_inds.append(Rmin_ind[0][0])
Rmax_dist_min = np.min(Rmax_dists)
Rmax_ind = np.where(Rmax_dists == Rmax_dist_min)
Rmax_inds.append(Rmax_ind[0][0])
azalias = []
for eigvecs, Rmin_ind, Rmax_ind in zip(eigvec_list, Rmin_inds, Rmax_inds):
print("azalia planted")
top_one = eigvecs[Rmin_ind:Rmax_ind, :, -1:]
hyperplanes = [top_one[j, :, :] for j in range(top_one.shape[0])]
azalias.append(iga(hyperplanes))
num_genes = 3000
ind_to_gene = get_index_to_gene()
def plot_loadings(ax, azalia, num_loadings=10):
gene_indices = list(range(num_genes))
mags = np.sqrt(np.sum(np.square(azalia), axis=1))
sorted_mags = np.sort(mags)[::-1]
gene_indices.sort(key=lambda x: mags[x], reverse=True)
inds_for_loading = gene_indices[:num_loadings]
genes_for_loading = [ind_to_gene[ind] for ind in inds_for_loading]
for ind, gene, height in zip(inds_for_loading, genes_for_loading, range(num_loadings)):
ax.plot(azalia[ind, 0], 0, "bo", ms=5)
ax.plot([azalia[ind, 0]]*2, [0, 0.01*height], c="gray", ls="--")
ax.annotate(gene, (azalia[ind, 0], 0.01*height))
ax.axvline(0, color="k")
ax.axhline(0, color="k")
figs = []
axx = []
for macro_ind, azalia in zip(indices, azalias):
fig = plt.figure()
figs.append(fig)
ax = fig.add_subplot(111)
axx.append(ax)
plot_loadings(ax, azalia)
ax.set_title("Cheerio "+str(macro_ind))
for j in range(len(azalias) - 1):
trans_vec = center_list[j+1] - center_list[j]
dir_vec = trans_vec / np.linalg.norm(trans_vec)
azalia_pre = azalias[j]
azalia_post = azalias[j+1]
pre_vec = azalia_pre.T @ dir_vec.T
post_vec = azalia_post.T @ dir_vec.T
axx[j].arrow(0, -0.01, pre_vec[0, 0], 0, color="r", head_width=0.01, head_starts_at_zero=True, length_includes_head=True)
axx[j+1].arrow(post_vec[0, 0], -0.02, -post_vec[0, 0], 0, color="b", head_width=0.01, head_starts_at_zero=True, length_includes_head=True)
for j, fig in zip(indices, figs):
fig.savefig("chain_1d_scale/"+str(j)+".jpg")