fixed 2 major PCA bugs

Author: aetherspritee

CCA.py

@@ -2,6 +2,7 @@
 import copy
 import csv
 import numba
+from numpy.core.multiarray import ndarray
 from PCA import PCA_subcluster
 
 def CCA_subcluster(R: np.ndarray, N: int, DF: float, kf: float,iter: int,N_subcl_perc: float ,ext_case: int ,tolerance: float=1e-7) -> tuple[bool,bool]:
@@ -243,28 +244,28 @@ def generate_CCA_pairs(I_t: int, i_orden: np.ndarray, X: np.ndarray,Y: np.ndarra
         ID_agglom[loc,loc] = 1
     return ID_agglom, CCA_ok
 
-def CCA_agg_properties(X: np.ndarray, Y: np.ndarray, Z: np.ndarray, R: np.ndarray, npp: int, Df: float, kf: float):
+def CCA_agg_properties(X: np.ndarray, Y: np.ndarray, Z: np.ndarray, R: np.ndarray, npp: int, Df: float, kf: float) -> tuple[float ,float, float, float, float, float]:
     m_vec = np.zeros((npp))
     R_i = np.zeros((npp))
-    sum_xm = 0
-    sum_ym = 0
-    sum_zm = 0
+    sum_xm: float = 0
+    sum_ym: float = 0
+    sum_zm: float = 0
     for i in range(npp):
         m_vec[i] = 4/3 * np.pi * np.power(R[i], 3)
         sum_xm += X[i]*m_vec[i]
         sum_ym += Y[i]*m_vec[i]
         sum_zm += Z[i]*m_vec[i]
 
-    m = np.sum(m_vec)
+    m = float(np.sum(m_vec))
     X_cm = sum_xm/m
     Y_cm = sum_ym/m
     Z_cm = sum_zm/m
 
-    rg = (np.exp(np.sum(np.log(R)/np.log(R).size)))*np.power(npp/kf, 1/Df)
+    rg = float((np.exp(np.sum(np.log(R)/np.log(R).size)))*np.power(npp/kf, 1/Df))
     for i in range(npp):
         R_i[i] = np.sqrt(np.power(X_cm-X[i],2) + np.power(Y_cm-Y[i],2) + np.power(Z_cm-Z[i],2))
 
-    R_max = np.max(R_i)
+    R_max = float(np.max(R_i))
     return rg, R_max, m, X_cm, Y_cm, Z_cm
 
 def CCA_identify_monomers(i_orden: np.ndarray):
@@ -276,7 +277,7 @@ def CCA_identify_monomers(i_orden: np.ndarray):
     return ID_mon
 
 @numba.njit()
-def CCA_random_select_list(X1: np.ndarray, Y1: np.ndarray, Z1: np.ndarray, R1: np.ndarray, X_cm1: np.ndarray, Y_cm1: np.ndarray, Z_cm1: np.ndarray, X2: np.ndarray, Y2: np.ndarray, Z2: np.ndarray,R2: np.ndarray, X_cm2: np.ndarray, Y_cm2: np.ndarray, Z_cm2: np.ndarray, curr_list: np.ndarray, gamma_pc: float, gamma_real: bool, ext_case: int):
+def CCA_random_select_list(X1: np.ndarray, Y1: np.ndarray, Z1: np.ndarray, R1: np.ndarray, X_cm1: float, Y_cm1: float, Z_cm1: float, X2: np.ndarray, Y2: np.ndarray, Z2: np.ndarray,R2: np.ndarray, X_cm2: float, Y_cm2: float, Z_cm2: float, curr_list: np.ndarray, gamma_pc: float, gamma_real: bool, ext_case: int):
     if gamma_real and ext_case == 1:
         for i in range(curr_list.shape[0]-1):
             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]
@@ -392,6 +393,7 @@ def CCA(X: np.ndarray,Y: np.ndarray,Z: np.ndarray,R: np.ndarray, N: int, ID_mon:
         prev_cand2 = 0
         cov_max = 1
 
+        # TODO: this here should go into its own function. run different tries in parallel, check if it fits, then break out.
         while cov_max > tolerance:
             if np.sum(curr_list) > 1:
                 prev_cand1 = CCA_random_pick(curr_list,prev_cand1)
@@ -434,6 +436,7 @@ def CCA(X: np.ndarray,Y: np.ndarray,Z: np.ndarray,R: np.ndarray, N: int, ID_mon:
                 curr_try += 1
 
                 if int(np.mod(curr_try, 359)) == 0 and np.sum(curr_list[prev_cand1,:]) > 1:
+                    print("OOF")
                     prev_cand2 = CCA_random_pick(curr_list,prev_cand1,prev_cand2)
 
                     COR1[:,0] = X1
@@ -545,7 +548,7 @@ def CCA_random_pick(curr_list: np.ndarray, prev_cand1: int, prev_cand2=None):
         prev_cand2 = selected_real
         return prev_cand2
 
-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):
+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) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray,np.ndarray, np.ndarray]:
     if gamma_real:
         X1 = COR1[:,0]
         Y1 = COR1[:,1]
@@ -605,6 +608,7 @@ def CCA_sticking_process(gamma_real: bool, gamma_pc: float, COR1,COR2,CM1,CM2, p
         sphere1 = np.array([x1_sph1, y1_sph1, z1_sph1, d1_min, d1_max])
         sphere2 = np.array([x2_sph2, y2_sph2, z2_sph2, d2_min, d2_max])
 
+        u_s1_cm1 = np.array([X_cm1, Y_cm1, Z_cm1])
         if ext_case == 1:
             if np.abs(d2_max-d1_max) < gamma_pc:
                 case = 1
@@ -697,7 +701,10 @@ def CCA_sticking_process(gamma_real: bool, gamma_pc: float, COR1,COR2,CM1,CM2, p
 
         return COR1, COR2, CM2, vec0, i_vec, j_vec
     else:
-        pass
+        vec0 = np.zeros(())
+        i_vec = np.zeros(())
+        j_vec = np.zeros(())
+        return COR1, COR2, CM2, vec0, i_vec, j_vec
 
 def random_point_SC(case: int,sphere1: np.ndarray, sphere2: np.ndarray):
     phi_crit_max = 0

PCA.py

@@ -1,40 +1,55 @@
 import numpy as np
 import numba
 import time
+from time import sleep
+import pyvista as pv
+
 
 def PCA(number: int, mass: np.ndarray, r: np.ndarray, Df: float, kf: float, tolerance: float) -> tuple[bool, np.ndarray]:
     PCA_ok = True
     n1,m1,rg1,x_cm,y_cm,z_cm, X,Y,Z = First_two_monomers(r, mass, number, Df, kf)
 
     if number > 2:
         k=3
-        while k < number+1:
+        while k <= number:
             n2 = 1
             m2 = mass[k-1]
 
             rg2 = np.sqrt(0.6)*r[k-1]
 
             n3 = n1+n2
             m3 = m1+m2
-
             rg3 = (np.exp(np.sum(np.log(r))/r.size))*np.power(n3/kf,1/Df)
+
             gamma_ok, gamma_pc = gamma_calc(rg1,rg2,rg3,n1,n2,n3)
             monomer_candidates = np.zeros((number-k+1)) # 0 = not considered
             monomer_candidates[0] = 1 # first one has been considered
+
             candidates, rmax = random_list_selection(gamma_ok, gamma_pc, X,Y,Z,r,n1,x_cm,y_cm,z_cm)
-            list_sum = 0
 
+            list_sum = 0
+            # print("here")
             while list_sum == 0:
-                while np.sum(candidates) == 0 and np.sum(monomer_candidates) < number-k:
+                # print(f"{np.sum(candidates) == 0 = }")
+                # print(f"{np.sum(monomer_candidates) <= number-k = }")
+                # if not np.sum(monomer_candidates) <= number-k-1:
+                #     print(f"{number = }")
+                #     print(f"{k = }")
+                #     print(f"{number-k = }")
+                #     print(f"{list_sum = }")
+                #     print(f"{np.sum(candidates) = }")
+                #     print(f"{monomer_candidates = }")
+                while np.sum(candidates) == 0 and np.sum(monomer_candidates) <= number-k:
                     candidates, rmax = search_monomer_candidates(r,mass,monomer_candidates,number,k,n3,Df,kf,rg1,n1,X,Y,Z,x_cm,y_cm,z_cm)
 
-                previous_candidate = 0 
+                previous_candidate = -1
 
                 if np.sum(candidates) > 0:
                     candidates, selected_real = random_list_selection_one(candidates, previous_candidate)
                     previous_candidate = selected_real
                 elif np.sum(candidates) == number-k+1:
                     PCA_ok = False
+                    exit(-1)
 
                 curr_try = 1
                 X_sel = X[selected_real]
@@ -43,13 +58,28 @@ def PCA(number: int, mass: np.ndarray, r: np.ndarray, Df: float, kf: float, tole
                 R_sel = r[selected_real]
                 r_k = r[k-1]
 
+                # print(f"modified {k-1 = }")
                 x_k, y_k, z_k, r0,x0,y0,z0,i_vec,j_vec = sticking_process(X_sel,Y_sel,Z_sel,R_sel,r_k,x_cm,y_cm,z_cm,gamma_pc)
                 X[k-1] = x_k
                 Y[k-1] = y_k
                 Z[k-1] = z_k
 
                 cov_max = overlap_check(X[0:k], Y[0:k], Z[0:k], r[0:k],k)
 
+                positions = [[X[i], Y[i], Z[i]] for i in range(X.shape[0])]
+                # print(f"{len(positions) = }")
+                # print(f"{positions}")
+                # point_cloud = pv.PolyData(positions)
+                # point_cloud["radius"] = [2 for i in range(X.shape[0])]
+
+                # geom = pv.Sphere(theta_resolution=8, phi_resolution=8)
+                # glyphed = point_cloud.glyph(scale="radius", geom=geom,)
+                #     # p = pv.Plotter(notebook=False, off_screen=True, window_size=(2000,2000))
+                # p = pv.Plotter(notebook=False)
+                # p.add_mesh(glyphed, color='white',smooth_shading=True)
+                # p.show()
+                # if cov_max > tolerance:
+                #     print("Need to readjust")
                 while cov_max > tolerance and curr_try < 360:
                     x_k, y_k, z_k,_ = sticking_process2(x0,y0,z0,r0,i_vec,j_vec)
 
@@ -61,6 +91,8 @@ def PCA(number: int, mass: np.ndarray, r: np.ndarray, Df: float, kf: float, tole
 
                     if np.mod(curr_try,359) == 0 and np.sum(candidates) > 1:
                         candidates, selected_real = random_list_selection_one(candidates, previous_candidate)
+                        # if np.sum(candidates) == 0:
+                        #     print("FOURTH, candidates all ZERO")
                         X_sel = X[selected_real]
                         Y_sel = Y[selected_real]
                         Z_sel = Z[selected_real]
@@ -73,15 +105,27 @@ def PCA(number: int, mass: np.ndarray, r: np.ndarray, Df: float, kf: float, tole
                         previous_candidate = selected_real
                         curr_try += 1
 
+                        # print("huh1")
                         cov_max = overlap_check(X[0:k], Y[0:k], Z[0:k], r[0:k],k)
+                        # print("huh2")
 
+                # FIXME: candidates may be full of zeros at some point which causes the program to get stuck
                 list_sum = np.sum(candidates)
+                if np.sum(candidates) == 0:
+                    print("FIFTH, candidates all ZERO")
+                    sleep(2)
+
 
                 if cov_max > tolerance:
                     list_sum = 0
                     candidates *= 0
+                    # print("SIXTH, set all candidates to ZERO")
+                    if number == k:
+                        print("Failure -- restarting PCA routine")
+                        return False, np.zeros((number,4))
 
 
+            # print("bruh")
             x_cm = (x_cm*m1 + X[k-1]*m2)/(m1+m2)
             y_cm = (y_cm*m1 + Y[k-1]*m2)/(m1+m2)
             z_cm = (z_cm*m1 + Z[k-1]*m2)/(m1+m2)
@@ -91,10 +135,21 @@ def PCA(number: int, mass: np.ndarray, r: np.ndarray, Df: float, kf: float, tole
             rg1 = (np.exp(np.sum(np.log(r))/np.log(r).size))*(np.power(n1/kf, 1/Df))
             k = k + 1
             # FIXME: this gets stuck for some reason (no apparent connection to the arcos issue)
+            # print("hi")
 
     data_new = np.zeros((number,4))
     for i in range(number):
         data_new[i,:] = np.array([X[i], Y[i], Z[i], r[i]])
+
+    # point_cloud = pv.PolyData(data_new[:,:-1])
+    # point_cloud["radius"] = [2 for i in range(X.shape[0])]
+
+    # geom = pv.Sphere(theta_resolution=8, phi_resolution=8)
+    # glyphed = point_cloud.glyph(scale="radius", geom=geom,)
+    #     # p = pv.Plotter(notebook=False, off_screen=True, window_size=(2000,2000))
+    # p = pv.Plotter(notebook=False)
+    # p.add_mesh(glyphed, color='white',smooth_shading=True)
+    # p.show()
     return PCA_ok, data_new
 
 def PCA_subcluster(N: int, N_subcluster: int, R: np.ndarray, DF: float, kf: float, tolerance: float) -> tuple[bool, np.ndarray, int, np.ndarray]:
@@ -120,7 +175,11 @@ def PCA_subcluster(N: int, N_subcluster: int, R: np.ndarray, DF: float, kf: floa
         for j in range(radius.size):
             mass[j] = 4/3 * np.pi * np.power(R[j],3)
 
-        PCA_OK, data_new = PCA(number,mass,radius,DF,kf,tolerance)
+        PCA_OK = False
+        while not PCA_OK:
+            PCA_OK, data_new = PCA(number,mass,radius,DF,kf,tolerance)
+        print(f"PCA LOOP {i} DONE!")
+        # time.sleep(1)
 
         if i == 0:
             acum = number
@@ -145,10 +204,16 @@ def First_two_monomers(R: np.ndarray,M: np.ndarray,N: int,DF: float,kf:float) ->
 
     u = np.random.rand()
     v = np.random.rand()
+    # print("first two monomers")
+    # print(f"{u = }, {v = }")
+
+
+    # u = 0.007316022089059793
+    # v = 0.43406326698233766
     theta = 2*np.pi*u
     phi = np.arccos(2*v-1)
-    theta = 1
-    phi = 1
+    # theta = 0.5
+    # phi = 0.1
 
     X[1] = X[0] + (R[0]+R[1])*np.cos(theta)*np.sin(phi)
     Y[1] = Y[0] + (R[0]+R[1])*np.sin(theta)*np.sin(phi)
@@ -158,11 +223,17 @@ def First_two_monomers(R: np.ndarray,M: np.ndarray,N: int,DF: float,kf:float) ->
     n1 = 2
 
     rg1 = (np.exp(np.sum(np.log(R[:2]))/2))*np.power(n1/kf,1/DF)
+    # print(f"{rg1 = }")
 
     x_cm = (X[0]*M[0]+X[1]*M[1])/(M[0] + M[1])
     y_cm = (Y[0]*M[0]+Y[1]*M[1])/(M[0] + M[1])
     z_cm = (Z[0]*M[0]+Z[1]*M[1])/(M[0] + M[1])
 
+
+    # print(f"{X[0] = }, {Y[0] = }, {Z[0] = }")
+    # print(f"{x_cm = }, {y_cm = }, {z_cm = }")
+    # print(f"{M = }")
+
     return n1,m1,rg1,x_cm,y_cm,z_cm, X,Y,Z
 
 def gamma_calc(rg1: float,rg2: float,rg3: float,n1: int,n2: int,n3: int) -> tuple[bool,float]:
@@ -178,7 +249,7 @@ def gamma_calc(rg1: float,rg2: float,rg3: float,n1: int,n2: int,n3: int) -> tupl
         gamma_pc = np.sqrt((np.power(n3,2)*np.power(rg3_aux,2)-n3*(n1*np.power(rg1,2)+n2*np.power(rg2,2)))/(n1*n2))
     else:
         gamma_ok = False
-    
+
     return gamma_ok, gamma_pc
 
 @numba.njit()
@@ -202,7 +273,8 @@ def search_monomer_candidates(R: np.ndarray, M: np.ndarray, monomer_candidates:
     M_sl = M 
     vector_search = np.zeros((N-k+1))
     for i in range(vector_search.size):
-        vector_search[i] = i+k
+        vector_search[i] = i+k-2
+        # print(f"{i+k-2 = }")
 
     for i in range(monomer_candidates.size):
         if monomer_candidates[i] == 1:
@@ -211,6 +283,8 @@ def search_monomer_candidates(R: np.ndarray, M: np.ndarray, monomer_candidates:
 
     if vector_search2.size > 1:
         u = np.random.rand()
+        # print("search monomer candidates")
+        # print(f"{u = }")
         RS_1 = int(vector_search2[int(u*(vector_search.size-1))])
     else:
         RS_1 = int(vector_search2[0])
@@ -229,15 +303,16 @@ def search_monomer_candidates(R: np.ndarray, M: np.ndarray, monomer_candidates:
 
     candidates, rmax = random_list_selection(gamma_ok, gamma_pc, X,Y,Z,R,n1,x_cm,y_cm,z_cm)
 
-    # FIXME: this fails once kf > 1.5 LMFAO
-    candidates[RS_1-k-1] = 1
     return candidates, rmax
 
 def random_list_selection_one(candidates: np.ndarray, previous_candidate: int):
-    if previous_candidate > 0:
+    if previous_candidate > -1:
         candidates[previous_candidate] = 0
     candidates2 = candidates[candidates > 0]
     n = np.random.rand()
+    # print("random list selection one")
+    # print(f"{n = }")
+    # n = 0.5
     selected = 1+int(n*(candidates2.size-1))
 
     selected_real = 0
@@ -260,6 +335,9 @@ def sticking_process(x: float,y: float,z: float,r: float,r_k: float, x_cm: float
     z2 = z_cm
     r2 = gamma_pc
 
+    # print(f"{x1 = }, {y1 = }, {z1 = }, {r1 = }")
+    # print(f"{x2 = }, {y2 = }, {z2 = }, {r2 = }")
+
     a = 2*(x2-x1)
     b = 2*(y2-y1)
     c = 2*(z2-z1)
@@ -273,9 +351,24 @@ def sticking_process(x: float,y: float,z: float,r: float,r_k: float, x_cm: float
 
     distance = np.sqrt(np.power(x2-x1,2) + np.power(y2-y1,2) + np.power(z2-z1,2))
 
+
+    # print(f"{distance = }")
+
     alpha = np.arccos((np.power(r1,2) + np.power(distance,2) - np.power(r2,2))/(2*r1*distance))
+    # print(f"{alpha = }")
     if abs((np.power(r1,2) + np.power(distance,2) - np.power(r2,2))/(2*r1*distance)) > 1:
         # FIXME: this should never happen!
+
+        pos = [[x1,y1,z1], [x2,y2,z2]]
+        point_cloud = pv.PolyData(pos)
+        point_cloud["radius"] = [2,2]
+
+        geom = pv.Sphere(theta_resolution=8, phi_resolution=8)
+        glyphed = point_cloud.glyph(scale="radius", geom=geom,)
+            # p = pv.Plotter(notebook=False, off_screen=True, window_size=(2000,2000))
+        p = pv.Plotter(notebook=False)
+        p.add_mesh(glyphed, color='white',smooth_shading=True)
+        p.show()
         print(f"{(np.power(r1,2) + np.power(distance,2) - np.power(r2,2))/(2*r1*distance) = }")
         exit(-1)
     r0 = r1*np.sin(alpha)
@@ -289,6 +382,9 @@ def sticking_process(x: float,y: float,z: float,r: float,r_k: float, x_cm: float
 
     u = np.random.rand()
     v = np.random.rand()
+    # print("sticking process")
+    # print(f"{u = }, {v = }")
+
     theta = 2.*np.pi*u
     phi = np.arccos(2.*v-1.)
 
@@ -300,6 +396,8 @@ def sticking_process(x: float,y: float,z: float,r: float,r_k: float, x_cm: float
 
 def sticking_process2(x0, y0, z0, r0,i_vec,j_vec):
     u = np.random.rand()
+    # print("sticking process2")
+    # print(f"{u = }")
     theta = 2 * np.pi * u 
 
     x_k = x0 + r0*np.cos(theta)*i_vec[0]+r0*np.sin(theta)*j_vec[0]