Skip to content

Commit eb855c1

Browse files
aetherspriteeaetherspritee
committed
last fixes + clean up
1 parent bf265ee commit eb855c1

File tree

1 file changed

+4
-79
lines changed

1 file changed

+4
-79
lines changed

‎CCA.py‎

Lines changed: 4 additions & 79 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,6 @@ def CCA_subcluster(R: np.ndarray, N: int, DF: float, kf: float,iter: int,N_subcl
2121

2222
PCA_OK, data, n_clusters, i_orden = PCA_subcluster(N, N_subcluster, R, DF, kf, tolerance)
2323

24-
# print(f"{data}")
2524
if not PCA_OK:
2625
return CCA_OK, PCA_OK
2726

@@ -58,7 +57,7 @@ def CCA_subcluster(R: np.ndarray, N: int, DF: float, kf: float,iter: int,N_subcl
5857

5958
other = 0
6059
i_orden = np.zeros((int(number_pairs),3))
61-
while k < I_total:
60+
while k <= I_total:
6261
for i in range(ID_agglom[k-1,:].size):
6362
if ID_agglom[k-1,i] == 1:
6463
other = i+1
@@ -71,9 +70,6 @@ def CCA_subcluster(R: np.ndarray, N: int, DF: float, kf: float,iter: int,N_subcl
7170
break
7271

7372
if k != other and IS_EMPTY:
74-
print("==================================")
75-
print(f"{I_total = }")
76-
print("==================================")
7773
Xn, Yn, Zn, Rn, CCA_OK = CCA(X,Y,Z,R, N, ID_mon, k,other,DF,kf,ext_case, tolerance)
7874

7975
considered[k-1] = 1
@@ -149,7 +145,6 @@ def CCA_subcluster(R: np.ndarray, N: int, DF: float, kf: float,iter: int,N_subcl
149145
R = R_next
150146

151147
iteration += 1
152-
print(f"{iteration = }")
153148

154149
for k in range(N-1):
155150
if np.isnan(X[k]) or np.isnan(Y[k]):
@@ -211,22 +206,15 @@ def generate_CCA_pairs(I_t: int, i_orden: np.ndarray, X: np.ndarray,Y: np.ndarra
211206
R2[jjj] = R[jj-1]
212207
jjj += 1
213208

214-
# print("WAZZUP")
215-
# print(f"{X2 = }, {Y2 = }, {Z2 = }, {R2 = }")
216209
rg2, r2_max, m2, _, _, _ = CCA_agg_properties(X2,Y2, Z2, R2,jjj, Df, kf)
217210

218211
m3 = m1+m2
219212
r_com = np.hstack((R1,R2))
220213
rg3 = (np.exp(np.sum(np.log(r_com))/((np.log(r_com).size))))*(((R1.size+R2.size))/kf)**(1./Df)
221214
if np.power(m3,2)*np.power(rg3,2) > m3 * (m1*np.power(rg1,2) + m2 * np.power(rg2,2)):
222-
print("GAMMA REAL")
223-
# print(f"{m1 = }, {m2 = }, {m3 = }")
224-
# print(f"{rg1 = }, {rg2 = }, {rg3 = }")
225215
gamma_pc = np.sqrt((np.power(m3,2)*np.power(rg3,2) - m3*(m1 * np.power(rg1,2) + m2 * np.power(rg2,2)))/(m1*m2))
226-
print(f"{gamma_pc = }")
227216
gamma_real = True
228217
else:
229-
print("GAMMA FALSE")
230218
gamma_pc = np.inf
231219
gamma_real = False
232220

@@ -281,7 +269,6 @@ def CCA_identify_monomers(i_orden: np.ndarray):
281269

282270
def CCA_random_select_list(X1, Y1, Z1, R1, X_cm1, Y_cm1, Z_cm1, X2, Y2, Z2,R2, X_cm2, Y_cm2, Z_cm2, curr_list: np.ndarray, gamma_pc: float, gamma_real: bool, ext_case):
283271
if gamma_real and ext_case == 1:
284-
print("EXT CASE 1 YOOOOOOOOOO")
285272
for i in range(curr_list.shape[0]-1):
286273
d_i_min = np.sqrt(np.power(X1[i]-X_cm1,2) + np.power(Y1[i]-Y_cm1,2) + np.power(Z1[i]-Z_cm1,2)) - R1[i]
287274
d_i_max = np.sqrt(np.power(X1[i]-X_cm1,2) + np.power(Y1[i]-Y_cm1,2) + np.power(Z1[i]-Z_cm1,2)) + R1[i]
@@ -338,11 +325,13 @@ def CCA(X: np.ndarray,Y: np.ndarray,Z: np.ndarray,R: np.ndarray, N: int, ID_mon:
338325
if ID_mon[i]+1 == other:
339326
monomers_2 += 1
340327

328+
341329
X2 = np.zeros((monomers_2))
342330
Y2 = np.zeros((monomers_2))
343331
Z2 = np.zeros((monomers_2))
344332
R2 = np.zeros((monomers_2))
345333

334+
346335
monomers_2 = 0
347336

348337
for i in range(N):
@@ -381,7 +370,6 @@ def CCA(X: np.ndarray,Y: np.ndarray,Z: np.ndarray,R: np.ndarray, N: int, ID_mon:
381370
CM2[1] = Y_cm2
382371
CM2[2] = Z_cm2
383372

384-
print(f"{CM2 = }")
385373
curr_list = np.zeros((n1,n2))
386374
curr_list = CCA_random_select_list(X1,Y1,Z1,R1,X_cm1,Y_cm1,Z_cm1,X2,Y2,Z2,R2,X_cm2,Y_cm2,Z_cm2,curr_list,gamma_pc,gamma_real,ext_case)
387375

@@ -417,10 +405,6 @@ def CCA(X: np.ndarray,Y: np.ndarray,Z: np.ndarray,R: np.ndarray, N: int, ID_mon:
417405

418406
COR1, COR2, CM2, vec0, i_vec, j_vec = CCA_sticking_process(gamma_real, gamma_pc, COR1, COR2, CM1, CM2, prev_cand1, prev_cand2, ext_case,n1,n2)
419407

420-
# FIXME: the following point coordinates differ !!!!!!!
421-
# this leads to cov_max being off
422-
423-
424408
X1 = COR1[:,0]
425409
Y1 = COR1[:,1]
426410
Z1 = COR1[:,2]
@@ -434,23 +418,15 @@ def CCA(X: np.ndarray,Y: np.ndarray,Z: np.ndarray,R: np.ndarray, N: int, ID_mon:
434418
Z_cm2 = CM2[2]
435419

436420
cov_max = CCA_overlap_check(n1,n2,X1,X2,Y1,Y2,Z1,Z2,R1,R2)
437-
print("PRE WHILE YOOOOOOO")
438-
print(f"{cov_max = }")
439-
print(f"{tolerance = }")
440-
print(f"{curr_try = }")
441421

442422
while cov_max > tolerance and curr_try < 360:
443-
# print("MMMMMMMOOOOOOOOOOOIN")
444423
X2,Y2, Z2 = CCA_sticking_process_v2(CM2, vec0, X2,Y2,Z2,i_vec,j_vec, prev_cand2)
445424
cov_max = CCA_overlap_check(n1,n2,X1,X2,Y1,Y2,Z1,Z2,R1,R2)
446425
curr_try += 1
447426

448427
if int(np.mod(curr_try, 359)) == 0 and np.sum(curr_list[prev_cand1,:]) > 1:
449-
print("YOOOOOOOOOOOO WHATS GOOOOOOOOOOOOOOO")
450428
prev_cand2 = CCA_random_pick(curr_list,prev_cand1,prev_cand2)
451429

452-
print(f"{X1.shape = }")
453-
print(f"{COR1[0].shape = }")
454430
COR1[:,0] = X1
455431
COR2[:,0] = X2
456432

@@ -487,9 +463,6 @@ def CCA(X: np.ndarray,Y: np.ndarray,Z: np.ndarray,R: np.ndarray, N: int, ID_mon:
487463
monomers_1 += 1
488464

489465
monomers_2 = 0
490-
print(f"{N-1 = }")
491-
print(f"{ID_mon = }")
492-
print(f"{k = }")
493466
for i in range(N-1):
494467
if ID_mon[i]+1 == other:
495468
X[i] = X2[monomers_2]
@@ -526,12 +499,9 @@ def CCA_random_pick(curr_list: np.ndarray, prev_cand1: int, prev_cand2=None):
526499
if prev_cand1 > 0:
527500
curr_list[prev_cand1,:] = curr_list[prev_cand1,:]*0
528501
list_sum = np.array([np.sum(curr_list[i,:]) for i in range(curr_list.shape[0])])
529-
# list_sum = curr_list[prev_cand1,:]
530502
curr_list2 = list_sum[list_sum > 0]
531503

532-
print("RANDOM!!!666666")
533504
uu = np.random.rand()
534-
# uu = 0.5
535505
selected = int(uu * (curr_list2.size-1))+1
536506
sel = 0
537507
jj = 0
@@ -544,18 +514,14 @@ def CCA_random_pick(curr_list: np.ndarray, prev_cand1: int, prev_cand2=None):
544514
selected_real = i
545515
break
546516
prev_cand1 = selected_real
547-
print(f"{prev_cand1 = }")
548517
return prev_cand1
549518
else:
550519
if prev_cand2 > 0:
551520
curr_list[prev_cand1,prev_cand2] = curr_list[prev_cand1,prev_cand2]*0
552-
# list_sum = np.array([np.sum(curr_list[:,i]) for i in range(curr_list.shape[1])])
553521
list_sum = curr_list[prev_cand1,:]
554522
curr_list2 = list_sum[list_sum > 0]
555523

556-
print("RANDOM!!111111111")
557524
uu = np.random.rand()
558-
# uu = 0.5
559525
selected = 1+int(uu * (curr_list2.size-1))
560526
sel = 0
561527
jj = 0
@@ -568,7 +534,6 @@ def CCA_random_pick(curr_list: np.ndarray, prev_cand1: int, prev_cand2=None):
568534
selected_real = i
569535
break
570536
prev_cand2 = selected_real
571-
print(f"{prev_cand2 = }")
572537
return prev_cand2
573538

574539
def CCA_sticking_process(gamma_real: bool, gamma_pc: float, COR1,COR2,CM1,CM2, prev_cand1: int, prev_cand2: int, ext_case: int,n1: int, n2: int):
@@ -600,10 +565,6 @@ def CCA_sticking_process(gamma_real: bool, gamma_pc: float, COR1,COR2,CM1,CM2, p
600565
vect_y /= vect_mag
601566
vect_z /= vect_mag
602567

603-
# center of mass of aggregate 2
604-
# print(f"{X_cm1 = }, {Y_cm1 = }, {Z_cm1 = }")
605-
# print(f"{vect_x = }, {vect_y = }, {vect_z = }")
606-
# print(f"{gamma_pc = }")
607568
x_cm22 = X_cm1 + gamma_pc*vect_x
608569
y_cm22 = Y_cm1 + gamma_pc*vect_y
609570
z_cm22 = Z_cm1 + gamma_pc*vect_z
@@ -622,28 +583,15 @@ def CCA_sticking_process(gamma_real: bool, gamma_pc: float, COR1,COR2,CM1,CM2, p
622583
y1_sph1 = Y_cm1
623584
z1_sph1 = Z_cm1
624585

625-
# NOTE: this matches with the original program
626586
d1_min = np.sqrt(np.power(X1[prev_cand1]-X_cm1,2) + np.power(Y1[prev_cand1]-Y_cm1,2) + np.power(Z1[prev_cand1]-Z_cm1,2)) - R1[prev_cand1]
627587
d1_max = np.sqrt(np.power(X1[prev_cand1]-X_cm1,2) + np.power(Y1[prev_cand1]-Y_cm1,2) + np.power(Z1[prev_cand1]-Z_cm1,2)) + R1[prev_cand1]
628588

629589
# sphere 2
630590
x2_sph2 = x_cm22
631591
y2_sph2 = y_cm22
632592
z2_sph2 = z_cm22
633-
# print(f"{R2 = }")
634-
# print(f"{X2_new = }")
635-
# print(f"{Y2_new = }")
636-
# print(f"{Z2_new = }")
637-
# print(f"{prev_cand2 = }")
638593
d2_min = np.sqrt(np.power(X2_new[prev_cand2]-x_cm22,2) + np.power(Y2_new[prev_cand2]-y_cm22,2) + np.power(Z2_new[prev_cand2]-z_cm22,2)) - R2[prev_cand2]
639-
# print(f"{d2_min = }")
640594
d2_max = np.sqrt(np.power(X2_new[prev_cand2]-x_cm22,2) + np.power(Y2_new[prev_cand2]-y_cm22,2) + np.power(Z2_new[prev_cand2]-z_cm22,2)) + R2[prev_cand2]
641-
# print(f"{np.power(X2_new[prev_cand2]-x_cm22,2) = }")
642-
# print(f"{np.power(Y2_new[prev_cand2]-y_cm22,2) = }")
643-
# print(f"{np.power(Z2_new[prev_cand2]-z_cm22,2) = }")
644-
# print(f"{np.sqrt(np.power(X2_new[prev_cand2]-x_cm22,2) + np.power(Y2_new[prev_cand2]-y_cm22,2) + np.power(Z2_new[prev_cand2]-z_cm22,2)) = }")
645-
# print(f"{d2_max = }")
646-
# print(f"{x_cm22 = }, {y_cm22 = }, {z_cm22 = }")
647595

648596
sphere1 = np.array([x1_sph1, y1_sph1, z1_sph1, d1_min, d1_max])
649597
sphere2 = np.array([x2_sph2, y2_sph2, z2_sph2, d2_min, d2_max])
@@ -667,7 +615,6 @@ def CCA_sticking_process(gamma_real: bool, gamma_pc: float, COR1,COR2,CM1,CM2, p
667615

668616
sphere1 = np.array([x1_sph1, y1_sph1, z1_sph1, sph1_r])
669617
sphere2 = np.array([x2_sph2, y2_sph2, z2_sph2, sph2_r])
670-
# print(f"{sphere2 = }")
671618
x,y,z,_,_,_ = CCA_2_sphere_intersec(sphere1,sphere2)
672619
u_s1_cm1 = np.array([X_cm1-x, Y_cm1-y, Z_cm1-z])
673620

@@ -701,8 +648,6 @@ def CCA_sticking_process(gamma_real: bool, gamma_pc: float, COR1,COR2,CM1,CM2, p
701648
Y1_new[i] = y_cm11 + new_c[1]
702649
Z1_new[i] = z_cm11 + new_c[2]
703650

704-
# print(f"{x_cm11 = }")
705-
# print(f"{new_c[0] = }")
706651
sph2_r = np.sqrt(np.power(X2_new[prev_cand2]-x_cm22,2) + np.power(Y2_new[prev_cand2]-y_cm22,2) + np.power(Z2_new[prev_cand2]-z_cm22,2))
707652
sph2_x = x_cm22
708653
sph2_y = y_cm22
@@ -733,8 +678,6 @@ def CCA_sticking_process(gamma_real: bool, gamma_pc: float, COR1,COR2,CM1,CM2, p
733678
Z2_new[i] = z_cm22 + new_c[2]
734679

735680
CM2 = np.array([x_cm22, y_cm22, z_cm22])
736-
# COR1 = np.zeros((X1_new.shape[0],3))
737-
# COR2 = np.zeros((X1_new.shape[0],3))
738681
COR1[:,0] = X1_new
739682
COR1[:,1] = Y1_new
740683
COR1[:,2] = Z1_new
@@ -746,7 +689,6 @@ def CCA_sticking_process(gamma_real: bool, gamma_pc: float, COR1,COR2,CM1,CM2, p
746689
return COR1, COR2, CM2, vec0, i_vec, j_vec
747690
else:
748691
pass
749-
# return COR1, COR2, CM2, vec0, i_vec, j_vec
750692

751693
def random_point_SC(case: int,sphere1: np.ndarray, sphere2: np.ndarray):
752694
phi_crit_max = 0
@@ -781,13 +723,9 @@ def random_point_SC(case: int,sphere1: np.ndarray, sphere2: np.ndarray):
781723
phi_crit_max = spherical_cap_angle(sphere1,sphere2)
782724
phi_crit_min = 0
783725
sph1_r = sph1_r_min
784-
print("RANDOM!!99999")
785726
uu = np.random.rand()
786-
# uu = 0.5
787727
theta_r = 2*np.pi*uu
788-
print("RANDOM!!8888")
789728
uu = np.random.rand()
790-
# uu = 0.5
791729
phi_r = phi_crit_min + (phi_crit_max-phi_crit_min)*uu
792730

793731
sph1_x = sphere1[0]
@@ -856,14 +794,7 @@ def CCA_2_sphere_intersec(sphere1: np.ndarray, sphere2: np.ndarray):
856794

857795
distance = np.sqrt(np.power(sphere2[0]-sphere1[0],2) + np.power(sphere2[1]-sphere1[1],2) + np.power(sphere2[2]-sphere1[2],2))
858796

859-
if np.abs((np.power(sphere1[3],2) + np.power(distance,2) - np.power(sphere2[3],2))/(2*sphere1[3]*distance)) > 1:
860-
# print(f"{sphere1 = }")
861-
# print(f"{sphere2 = }")
862-
# print(f"{distance = }")
863-
print(f"ACOS ARG TOO LARGE: {np.abs((np.power(sphere1[3],2) + np.power(distance,2) - np.power(sphere2[3],2))/(2*sphere1[3]*distance))}! EXITING!")
864-
exit(-1)
865797
alpha_0 = np.arccos((np.power(sphere1[3],2) + np.power(distance,2) - np.power(sphere2[3],2))/(2*sphere1[3]*distance))
866-
# print(f"{(np.power(sphere1[3],2) + np.power(distance,2) - np.power(sphere2[3],2))/(2*sphere1[3]*distance) = }")
867798
r0 = sphere1[3]*np.sin(alpha_0)
868799

869800
AmBdC = -(A+B)/C
@@ -872,9 +803,7 @@ def CCA_2_sphere_intersec(sphere1: np.ndarray, sphere2: np.ndarray):
872803
i_vec = np.array([1,1,AmBdC])/np.sqrt(1+1+np.power(AmBdC,2))
873804
j_vec = np.cross(k_vec,i_vec)
874805

875-
print("RANDOM!00000")
876806
uu = np.random.rand()
877-
# uu = 0.5
878807
theta = np.pi * 2 * uu
879808

880809
x = x0 + r0*np.cos(theta)*i_vec[0] + r0 * np.sin(theta)*j_vec[0]
@@ -897,9 +826,7 @@ def CCA_overlap_check(n1: int, n2: int, X1,X2,Y1,Y2,Z1,Z2,R1,R2):
897826
return cov_max
898827

899828
def CCA_sticking_process_v2(CM2: np.ndarray, vec0: np.ndarray, X2_new,Y2_new,Z2_new, i_vec, j_vec, prev_cand):
900-
# print("RANDOM!!77777777")
901829
uu = np.random.rand()
902-
# uu = 0.5
903830
theta_a = 2*np.pi * uu
904831

905832
x = vec0[0] + vec0[3] * np.cos(theta_a) * i_vec[0] + vec0[3] * np.sin(theta_a) * j_vec[0]
@@ -910,7 +837,6 @@ def CCA_sticking_process_v2(CM2: np.ndarray, vec0: np.ndarray, X2_new,Y2_new,Z2_
910837
v2 = np.array([x-CM2[0], y-CM2[1], z-CM2[2]])
911838
s_vec = np.cross(v1,v2)/np.linalg.norm(np.cross(v1,v2))
912839

913-
# FIXME: check for numerical issues here!
914840
if np.dot(v1,v2)/np.linalg.norm(np.dot(v1,v2)) > 1 or np.dot(v1,v2)/np.linalg.norm(np.dot(v1,v2)) < -1:
915841
angle = np.arccos(1)
916842
else:
@@ -931,7 +857,7 @@ def CCA_sticking_process_v2(CM2: np.ndarray, vec0: np.ndarray, X2_new,Y2_new,Z2_
931857
def save_results(X: np.ndarray, Y: np.ndarray, Z: np.ndarray, R: np.ndarray, iteration: int, res_name: str="test"):
932858
with open(res_name+str(iteration)+".csv", "w") as f:
933859
writer = csv.writer(f)
934-
for i in range(X.size-1):
860+
for i in range(X.size):
935861
writer.writerow([X[i], Y[i], Z[i], R[i]])
936862

937863
def sort_rows(i_orden: np.ndarray):
@@ -943,5 +869,4 @@ def sort_rows(i_orden: np.ndarray):
943869
i_orden[irow,:] = i_orden[krow,:]
944870
i_orden[krow,:] = temp
945871

946-
print(f"sortrows = {i_orden}")
947872
return i_orden

0 commit comments

Comments
 (0)