diff --git a/src/napari_stress/_measurements/geodesics.py b/src/napari_stress/_measurements/geodesics.py index 1f313c91..c76579ed 100644 --- a/src/napari_stress/_measurements/geodesics.py +++ b/src/napari_stress/_measurements/geodesics.py @@ -138,7 +138,6 @@ def correlation_on_surface( # calculate autocorrelation with spatially averaged feature values for dist_i in range(0, maximal_distance + 1): - sum_mean_curv_corr_d_i, n_points = _avg_around_pt( dist_i, dists_lbdv_non0, Corr_non0_pts, maximal_distance ) @@ -235,7 +234,6 @@ def local_extrema_analysis( delta_feature_nearest_min_max = np.zeros_like(feature) for pt in range(quad_fit): - H_pt = feature[pt] # set to True, change if False: diff --git a/src/napari_stress/_measurements/stresses.py b/src/napari_stress/_measurements/stresses.py index d6dff5e1..e10d0aa7 100644 --- a/src/napari_stress/_measurements/stresses.py +++ b/src/napari_stress/_measurements/stresses.py @@ -172,6 +172,7 @@ def calculate_anisotropy( every group in the dataframe """ + # write a function to apply to every group in the dataframe def anisotropy( df: pd.DataFrame, alpha: float = 0.05, column: str = "anisotropic_stress" diff --git a/src/napari_stress/_reconstruction/reconstruct_surface.py b/src/napari_stress/_reconstruction/reconstruct_surface.py index a68c6bc0..3453301b 100644 --- a/src/napari_stress/_reconstruction/reconstruct_surface.py +++ b/src/napari_stress/_reconstruction/reconstruct_surface.py @@ -46,7 +46,6 @@ def reconstruct_surface_from_quadrature_points(points: PointsData) -> SurfaceDat vertex = tetras[tri_i, tetra_vert] if vertex != n_quadrature_points and vert_ind < 3: - delauney_triangles[tri_i, vert_ind] = vertex vert_ind = vert_ind + 1 diff --git a/src/napari_stress/_reconstruction/refine_surfaces.py b/src/napari_stress/_reconstruction/refine_surfaces.py index af3efcc2..f678bedf 100644 --- a/src/napari_stress/_reconstruction/refine_surfaces.py +++ b/src/napari_stress/_reconstruction/refine_surfaces.py @@ -147,7 +147,6 @@ def trace_refinement_of_surface( # Iterate over all provided target points for idx in range(n_points): - array = np.array(intensity_along_vector.loc[idx].to_numpy()) # Simple or fancy fit? if selected_fit_type == fit_types.quick_edge_fit: @@ -329,7 +328,6 @@ def _fancy_edge_fit( array = [x for x in array if not np.isnan(x)] # filter out nans try: if selected_edge_func == _sigmoid: - # trim begin of trace to get rid of rising intensity slope ind_max = np.argmax(array) array = array[ind_max:] diff --git a/src/napari_stress/_spherical_harmonics/spherical_harmonics.py b/src/napari_stress/_spherical_harmonics/spherical_harmonics.py index e439b343..a9098416 100644 --- a/src/napari_stress/_spherical_harmonics/spherical_harmonics.py +++ b/src/napari_stress/_spherical_harmonics/spherical_harmonics.py @@ -55,7 +55,6 @@ def shtools_spherical_harmonics_expansion( longitude = np.rad2deg(spherical_coordinates[2]) if expansion_type == "radial": - # Find spherical harmonics expansion coefficients until specified degree opt_fit_params = pyshtools._SHTOOLS.SHExpandLSQ( radius, latitude, longitude, lmax=max_degree @@ -74,7 +73,6 @@ def shtools_spherical_harmonics_expansion( return points, coefficients elif expansion_type == "cartesian": - optimal_fit_params = [] fitted_coordinates = np.zeros_like(points) for i in range(3): @@ -267,7 +265,6 @@ def lebedev_quadrature( def create_manifold( points: PointsData, lebedev_fit: lebedev_info.lbdv_info, max_degree: int ) -> mnfd.manifold: - # add information to manifold on what type of expansion was used Manny_Dict = {} Manny_Name_Dict = {} # sph point cloud at lbdv diff --git a/src/napari_stress/_stress/charts_SPB.py b/src/napari_stress/_stress/charts_SPB.py index 169ca557..6861e83d 100644 --- a/src/napari_stress/_stress/charts_SPB.py +++ b/src/napari_stress/_stress/charts_SPB.py @@ -14,7 +14,6 @@ def Domain(Theta, Phi): # Which coordinates to use if Chart_Min_Polar <= Phi and Phi <= Chart_Max_Polar: # Use Coordinates A - Bar_Coors = Coor_A_To_B(Theta, Phi) Phi_Bar = Bar_Coors[1] @@ -32,7 +31,6 @@ def Domain(Theta, Phi): # Which coordinates to use def Domain_Unaffected(Theta, Phi): # Where fns in each chart arent affected by eta fns if eta_Min_Polar <= Phi and Phi <= eta_Max_Polar: # Use Coordinates A - Bar_Coors = Coor_A_To_B(Theta, Phi) Phi_Bar = Bar_Coors[1] @@ -97,7 +95,6 @@ def eta_A_const( def eta_B( func, Theta_Bar, Phi_Bar ): # Cutoff fn for Phi_Bar in [eta_Min_Polar, eta_Max_Polar] - # func needs to be COUNTER rotated back to phi, theta, to be consistent def func_counter_rot(theta_bar, phi_bar): Orig_Coors = Coor_B_To_A(theta_bar, phi_bar) @@ -111,7 +108,6 @@ def func_counter_rot(theta_bar, phi_bar): # NOTE: (Theta, Phi) -> # (Theta_Bar, Phi_Bar) by rotating the north pole (0,0) to the east pole (0,pi/2) def Coor_A_To_B(Theta, Phi): # Converts from A to B - X1 = np.cos(Theta) * np.sin(Phi) Y1 = np.sin(Theta) * np.sin(Phi) Z1 = np.cos(Phi) @@ -141,7 +137,6 @@ def Coor_A_To_B(Theta, Phi): # Converts from A to B def Coor_B_To_A(Theta_Bar, Phi_Bar): # Converts from B back to A - # X2 = cos(Theta_Bar)*sin(Phi_Bar) # Y2 = sin(Theta_Bar)*sin(Phi_Bar) # Z2 = cos(Phi_Bar) @@ -173,14 +168,12 @@ def Coor_B_To_A(Theta_Bar, Phi_Bar): # Converts from B back to A # Uses F(theta, phi) to find F_{Rot}(theta_bar, phi_bar) def Rotate_Fn(func, Theta_Bar, Phi_Bar): - Theta, Phi = Coor_B_To_A(Theta_Bar, Phi_Bar) return func(Theta, Phi) # (x,y,z) -> (theta, phi) def Cart_To_Coor_A(x, y, z): - r = np.sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) ##Calculate Theta @@ -200,7 +193,6 @@ def Cart_To_Coor_A(x, y, z): # (x,y,z) -> (theta_bar, phi_bar) def Cart_To_Coor_B(x, y, z): - r = np.sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)) ##Calculate Theta_Bar diff --git a/src/napari_stress/_stress/euclidian_k_form_SPB.py b/src/napari_stress/_stress/euclidian_k_form_SPB.py index 44ce1f91..187bd69a 100644 --- a/src/napari_stress/_stress/euclidian_k_form_SPB.py +++ b/src/napari_stress/_stress/euclidian_k_form_SPB.py @@ -12,14 +12,12 @@ def zero_quad_array(Q_val): # Return list of quad vals at each quad pt, given func def get_quadrature_points_from_function(Func, lbdv): - Extracted_Quad_vals = Func(lbdv.theta_pts, lbdv.phi_pts) return Extracted_Quad_vals # Return list of quad vals at each quad pt, given SPH_Func (CHART A order) def get_quadrature_points_from_sh_function(SPH_Func, lbdv, Chart): - Extracted_Quad_vals = zero_quad_array(lbdv.lbdv_quad_pts) quad_pts = range(lbdv.lbdv_quad_pts) @@ -37,13 +35,11 @@ def get_quadrature_points_from_sh_function(SPH_Func, lbdv, Chart): # Return list of quad vals at each quad pt within Chart (CHART A order) def Extract_dPhi_Quad_Pt_Vals_From_SPH_Fn(sph_func, lbdv, Chart): - # Extracted_dPhi_Quad_Vals = zero_quad_array(lbdv.lbdv_quad_pts) quad_pts = range(lbdv.lbdv_quad_pts) if Chart == "A": - return np.where( lbdv.Chart_of_Quad_Pts > 0, sph_func.Eval_SPH_Der_Phi_Coef(quad_pts, lbdv), @@ -51,7 +47,6 @@ def Extract_dPhi_Quad_Pt_Vals_From_SPH_Fn(sph_func, lbdv, Chart): ) if Chart == "B": - quad_pts_rot = lbdv.Eval_Rot_Lbdv_Quad_vals(quad_pts) return np.where( lbdv.Chart_of_Quad_Pts[quad_pts_rot] > 0, @@ -62,13 +57,11 @@ def Extract_dPhi_Quad_Pt_Vals_From_SPH_Fn(sph_func, lbdv, Chart): # Return list of quad vals at each quad pt within Chart (CHART A order) def Extract_dPhi_Phi_Quad_Pt_Vals_From_SPH_Fn(sph_func, lbdv, Chart): - # Extracted_dPhi_Phi_Quad_Vals = zero_quad_array(lbdv.lbdv_quad_pts) quad_pts = range(lbdv.lbdv_quad_pts) if Chart == "A": - return np.where( lbdv.Chart_of_Quad_Pts > 0, sph_func.Eval_SPH_Der_Phi_Phi_Coef(quad_pts, lbdv), @@ -76,7 +69,6 @@ def Extract_dPhi_Phi_Quad_Pt_Vals_From_SPH_Fn(sph_func, lbdv, Chart): ) if Chart == "B": - quad_pts_rot = lbdv.Eval_Rot_Lbdv_Quad_vals(quad_pts) return np.where( lbdv.Chart_of_Quad_Pts[quad_pts_rot] > 0, @@ -87,7 +79,6 @@ def Extract_dPhi_Phi_Quad_Pt_Vals_From_SPH_Fn(sph_func, lbdv, Chart): # Combines info from each chart to give good vals at every point def Combine_Chart_Quad_Vals(Quad_Vals_A, Quad_Vals_B, lbdv): - quad_pts = range(lbdv.lbdv_quad_pts) return np.where( lbdv.Chart_of_Quad_Pts > 0, @@ -98,7 +89,6 @@ def Combine_Chart_Quad_Vals(Quad_Vals_A, Quad_Vals_B, lbdv): # Combines K_A and K_B into Chart A vector def Combine_Manny_Gauss_Curvatures(Manny, lbdv, verbose=False): - K_cominbed_A_pts = Combine_Chart_Quad_Vals(Manny.K_A_pts, Manny.K_B_pts, lbdv) if verbose is True: @@ -116,7 +106,6 @@ def Combine_Manny_Gauss_Curvatures(Manny, lbdv, verbose=False): # Converts 1-form alpha, to vector field at quad pts: def Sharp(dtheta_A_vals, dphi_A_vals, dtheta_B_vals, dphi_B_vals, p_deg, lbdv, Manny): - # alpha_i = (dtheta_vals_i)dtheta + (d_phi_vals_i)dphi # Store x,y,z vals of sharped vector field at each quad_pt @@ -263,7 +252,6 @@ def Sharp(dtheta_A_vals, dphi_A_vals, dtheta_B_vals, dphi_B_vals, p_deg, lbdv, M # Converts 2-forms from Euclidean to Polar at a Quad Pt def Two_Form_Conv_to_Polar_pt(quad_pt, lbdv, Manny, Chart): - # BJG: DOES NOT NEED TO CHANGE FOR SPB FRAME: theta_pt = lbdv.theta_pts[quad_pt] @@ -274,7 +262,6 @@ def Two_Form_Conv_to_Polar_pt(quad_pt, lbdv, Manny, Chart): dx_dy_comp = [] if np.isscalar(quad_pt) is True: - theta_pt = np.asscalar(theta_pt) phi_pt = np.asscalar(phi_pt) @@ -288,7 +275,6 @@ def Two_Form_Conv_to_Polar_pt(quad_pt, lbdv, Manny, Chart): dx_dz_comp = -1 * Cross_of_Basis[1] dx_dy_comp = Cross_of_Basis[2] else: - # Cross_of_Basis = Manny.Normal_Dir_Quad_Pts(lbdv, Chart) # Manny.Normal_Dir(theta_pt, phi_pt, Chart) @@ -311,7 +297,6 @@ def Two_Form_Conv_to_Polar_pt(quad_pt, lbdv, Manny, Chart): def Gen_Curl_1( theta_C_A_vals, phi_C_A_vals, theta_C_B_vals, phi_C_B_vals, p_deg, lbdv, Manny ): - # Project into SPH within Charts: theta_C_A_SPH_Fn = sph_f.Faster_Double_Proj(theta_C_A_vals, p_deg, lbdv) phi_C_A_SPH_Fn = sph_f.Faster_Double_Proj(phi_C_A_vals, p_deg, lbdv) @@ -367,7 +352,6 @@ def Gen_Curl_1( # takes in quadrature points and weights them by metric: def Integral_on_Manny_Eta_Z(vals_at_quad_pts, Manny, lbdv): # Orriginal Version - met_fac_over_sin_phi_pts_A = Manny.Metric_Factor_A_over_sin_phi_pts met_fac_over_sin_phi_pts_B = Manny.Metric_Factor_B_over_sin_phi_bar_pts @@ -381,7 +365,6 @@ def Integral_on_Manny_Eta_Z(vals_at_quad_pts, Manny, lbdv): # Orriginal Version rotated_vals_at_quad_pts = np.zeros((num_quad_pts, 1)) for quad_pt in range(num_quad_pts): - eta_Z_pt = lbdv.eta_z(quad_pt) eta_Z_Chart_A[quad_pt, 0] = eta_Z_pt @@ -406,7 +389,6 @@ def Integral_on_Manny_Eta_Z(vals_at_quad_pts, Manny, lbdv): # Orriginal Version def Integral_on_Manny(vals_at_quad_pts, Manny, lbdv): # New Version - vals_at_quad_pts = vals_at_quad_pts.reshape( len(vals_at_quad_pts), 1 ) # ASSERT THIS IS THE RIGHT SIZE @@ -453,7 +435,6 @@ def lebedev_quad_adj_Manny(Manny, lbdv): def Lp_Rel_Error_At_Quad_On_Manny( approx_f_vals, f_vals, lbdv, p, Manny ): # Assumes f NOT 0 - # Lp_Err = 0 # ||self - f||_p # Lp_f = 0 # || f ||_p @@ -467,7 +448,6 @@ def Lp_Rel_Error_At_Quad_On_Manny( # Get rid of normal components in a elegant way: def Tangent_Projection(Vectors_at_Quad, lbdv, Manny): - num_quad_pts = lbdv.lbdv_quad_pts Tangent_Vecs_at_Quad = np.zeros((num_quad_pts, 3)) @@ -478,14 +458,12 @@ def Tangent_Projection(Vectors_at_Quad, lbdv, Manny): sigma_phi_B_pts = Manny.Sigma_Phi_B_Pts for quad_pt in range(num_quad_pts): - Vec_pt = Vectors_at_Quad[quad_pt, :].flatten() Normal_Vec_pt = [] # Chart A: if lbdv.Chart_of_Quad_Pts[quad_pt] > 0: - Normal_Vec_pt = np.cross( sigma_theta_A_pts[quad_pt, :].flatten(), sigma_phi_A_pts[quad_pt, :].flatten(), @@ -530,7 +508,6 @@ def Tangent_Projection(Vectors_at_Quad, lbdv, Manny): # Returns Riemannian Inner Product of 1-Forms input, returns Integral of def Riemann_L2_Inner_Product_One_Form(One_Form_pts_1, One_Form_pts_2, lbdv, manny): - One_Form_1 = euc_k_form( 1, lbdv.lbdv_quad_pts, manny.Man_SPH_Deg, manny, One_Form_pts_1 ) @@ -567,7 +544,6 @@ def Riemann_L2_Inner_Product_Zero_Form(Zero_Form_pts_1, Zero_Form_pts_2, lbdv, m class euc_k_form(object): def __init__(self, k_val, Q_val, P_deg, Manny, array_of_quad_vals): - self.k_value = k_val # k = 0, 1, 2 self.Q_value = Q_val # lbdv.lbdv_quad_pts # Number of Quad Pts self.P_degree = P_deg # Degree of Basis @@ -611,7 +587,6 @@ def times_const(self, Const): # Multiplies k_form by a function def times_fn(self, fn_vals_at_quad_pts): - num_cols = self.array_size product_array = np.zeros((self.Q_value, num_cols)) @@ -638,10 +613,8 @@ def Flat(self, G_deriv, A_inv_deriv, lbdv): One_Form_Phi_Comp_B_vals = zero_quad_array(self.Q_value) for quad_pt in range(self.Q_value): - # We use quad pts within each chart if lbdv.Chart_of_Quad_Pts[quad_pt] > 0: - quad_pt_rot = lbdv.Eval_Inv_Rot_Lbdv_Quad_vals(quad_pt) alpha_sharp_X_A_pt = self.quad_pt_array[quad_pt, :].T @@ -710,9 +683,7 @@ def Flat(self, G_deriv, A_inv_deriv, lbdv): # d: k_form -> (k+1)_form def Ext_Der(self, lbdv): - if self.k_value == 0: - # print("Doing Ext_Der on 0-Forms") f_dtheta_A_quad_vals = zero_quad_array(lbdv.lbdv_quad_pts) @@ -767,7 +738,6 @@ def Ext_Der(self, lbdv): ) if self.k_value == 1: - # print("Doing Ext_Der on 1-Forms") # Store values of Ext Det of (1-forms) alpha, at quad pts in each Chart @@ -925,7 +895,6 @@ def Ext_Der(self, lbdv): # Compute Y = -dA_i * A^(-1) def Y_mat(q_pt, Chart, deriv): - Basis_Mat = self.Manifold.Change_Basis_Mat(q_pt, Chart) dBasis_Mat_deriv = [] @@ -1043,7 +1012,6 @@ def Y_mat(q_pt, Chart, deriv): ) if self.k_value == 2: # ext der of of 2-form is 0: - # print("Doing Ext_Der on 2-Forms") return euc_k_form( @@ -1056,9 +1024,7 @@ def Y_mat(q_pt, Chart, deriv): # Compute *: k_form -> (n-k)_form def Hodge_Star(self, lbdv): - if self.k_value == 0: - # print("Doing Hodge Star on 0-Forms") # store values of Hodge star of fn (0-form) f, at quad pts in each Chart @@ -1125,7 +1091,6 @@ def Hodge_Star(self, lbdv): ) if self.k_value == 1: - # print("Doing Hodge Star on 1-Forms") # Store values of Hodge Star of (1-forms) alpha, at quad pts in each Chart @@ -1198,7 +1163,6 @@ def Hodge_Star(self, lbdv): ) if self.k_value == 2: - # print("Doing Hodge Star on 2-Forms") # store values of Hodge star of fn (2-form) beta, at quad pts in each Chart @@ -1270,9 +1234,7 @@ def Hodge_Star(self, lbdv): # Computes -*d directly; # Takes 0-Forms -> 1-Forms (in local chart) def Gen_Curl_0(self, lbdv): - if self.k_value == 0: - # print("Doing Ext_Der on 0-Forms") f_dtheta_A_quad_vals = zero_quad_array(lbdv.lbdv_quad_pts) @@ -1378,9 +1340,7 @@ def Gen_Curl_0(self, lbdv): # Computes -*d directly; Will give output k-forms def Gen_Curl_to_K_Form(self, lbdv): - if self.k_value == 0: - # print("Doing Ext_Der on 0-Forms") f_dtheta_A_quad_vals = zero_quad_array(lbdv.lbdv_quad_pts) @@ -1466,7 +1426,6 @@ def Gen_Curl_to_K_Form(self, lbdv): # div(v) = -\delta v^{\flat} def Divergence_1_Form(self, lbdv, debug_mode=False): if self.k_value == 1: - # get X derivatives at quad pts: Vec_X_SPH_Fn_A, Vec_X_SPH_Fn_B = sph_f.Proj_Into_SPH_Charts_At_Quad_Pts( self.quad_pt_array[:, 0], self.P_degree, lbdv @@ -1894,7 +1853,6 @@ def Divergence_1_Form(self, lbdv, debug_mode=False): # Write Expressions in 1 line (COPIED FROM LB_POINT_VALS_SPB): def Explicit_LB(self, lbdv): - Y_mn_SPH_Fn_A, Y_mn_SPH_Fn_B = sph_f.Proj_Into_SPH_Charts_At_Quad_Pts( self.quad_pt_array, self.P_degree, lbdv ) @@ -1980,7 +1938,6 @@ def Explicit_LB(self, lbdv): # we get better convergence this way # (by ignoring intermediate k-form stuctures) def LB_Zero_Form_From_Curl(self, lbdv): - if self.k_value != 0: print("Error: Zero Forms needed for this inner-product fn to work") @@ -2023,10 +1980,8 @@ def Laplace_Beltrami(self, lbdv): # Returns Riemannian Inner Product of 1-Forms input, # returns vector of scalar inner-products at each point def Riemann_Inner_Prouct_One_Form(self, One_Form_2, lbdv): - # if we have faulty input if self.k_value != 1 or One_Form_2.k_value != 1: - print("Error: One Forms needed for this inner-product fn to work") else: @@ -2042,7 +1997,6 @@ def Riemann_Inner_Prouct_One_Form(self, One_Form_2, lbdv): # Returns Riemannian Inner Product of 2-Forms with ITSELF input, # returns vector of scalar inner-products at each point def Riemann_Self_Inner_Prouct_Two_Form(self, lbdv): - # Corresponding quad pt in chart B (in Chart A Coors) quad_pts = range(self.Q_value) quad_pts_inv_rot = lbdv.Eval_Inv_Rot_Lbdv_Quad_vals(quad_pts) @@ -2103,13 +2057,10 @@ def Riemann_Self_Inner_Prouct_Two_Form(self, lbdv): # Pushes forward a 1-form onto Manny_2, using our chart convention def Push_Forward_One_Form(self, Manny_2, lbdv): - if self.k_value != 1: - print("Error: One Form needed for this Push-Forward fn to work") else: - push_forward_quad_vals = np.zeros((self.Q_value, 3)) sigma_theta_A_vecs = Manny_2.sigma_theta_Quad_Pts(lbdv, "A") @@ -2119,12 +2070,10 @@ def Push_Forward_One_Form(self, Manny_2, lbdv): sigma_phi_B_vecs = Manny_2.sigma_phi_Quad_Pts(lbdv, "B") for quad_pt in range(self.Q_value): - # print("q = "+str(quad_pt)) # If Chart A: if lbdv.Chart_of_Quad_Pts[quad_pt] > 0: - alpha_sharp_X_A_pt = self.quad_pt_array[quad_pt, :].T A_Mat_pt = self.Manifold.Change_Basis_Mat(quad_pt, "A") @@ -2148,7 +2097,6 @@ def Push_Forward_One_Form(self, Manny_2, lbdv): push_forward_quad_vals[quad_pt, :] = push_forward_vec_A_pt.T else: - quad_pt_rot = lbdv.Eval_Rot_Lbdv_Quad_vals(quad_pt) alpha_sharp_X_B_pt = self.quad_pt_array[quad_pt, :].T diff --git a/src/napari_stress/_stress/lebedev_info_SPB.py b/src/napari_stress/_stress/lebedev_info_SPB.py index 8b5e1a5f..d4eec9a1 100644 --- a/src/napari_stress/_stress/lebedev_info_SPB.py +++ b/src/napari_stress/_stress/lebedev_info_SPB.py @@ -86,6 +86,7 @@ 66: 5810, } + # Finds Lowest order basis for degree: def look_up_lbdv_pts(Degree): if Degree > 66 or Degree < 2: @@ -115,7 +116,6 @@ def Eval_SPH_Basis(M_Coef, N_Coef, Theta, Phi): # Evaluates d_phi(Y^m_n) at SINGLE PT def Der_Phi_Basis_Fn(M_Coef, N_Coef, Theta, Phi): # M_Coef < 0 Corresonds to Z^|M|_N - Der_Phi_Val = [] # For Scalar Case, we use usual vectorization: @@ -133,7 +133,6 @@ def Der_Phi_Basis_Fn(M_Coef, N_Coef, Theta, Phi): # M_Coef < 0 Corresonds to Z^ return 0 # d_phi Y^0_0 = 0 elif M_Coef < 0: - m_sph = -1 * M_Coef Der_Phi_Val += (m_sph * mpmath.cot(Phi)) * sph_harm( m_sph, N_Coef, Theta, Phi @@ -151,7 +150,6 @@ def Der_Phi_Basis_Fn(M_Coef, N_Coef, Theta, Phi): # M_Coef < 0 Corresonds to Z^ ) else: # M_Coef >= 0 - m_sph = M_Coef Der_Phi_Val += (m_sph * mpmath.cot(Phi)) * sph_harm( m_sph, N_Coef, Theta, Phi @@ -169,7 +167,6 @@ def Der_Phi_Basis_Fn(M_Coef, N_Coef, Theta, Phi): # M_Coef < 0 Corresonds to Z^ ) else: # array input case: - # COPIED FROM SPH_DER_PHI_FN Der_Phi_Val = np.zeros_like(Theta) @@ -183,7 +180,6 @@ def Der_Phi_Basis_Fn(M_Coef, N_Coef, Theta, Phi): # M_Coef < 0 Corresonds to Z^ return 0 # d_phi Y^0_0 = 0 elif M_Coef < 0: - m_sph = -1 * M_Coef Der_Phi_Val += (m_sph * (1.0 / np.tan(Phi))) * sph_harm( m_sph, N_Coef, Theta, Phi @@ -201,7 +197,6 @@ def Der_Phi_Basis_Fn(M_Coef, N_Coef, Theta, Phi): # M_Coef < 0 Corresonds to Z^ ) else: # M_Coef >= 0 - m_sph = M_Coef Der_Phi_Val += (m_sph * (1.0 / np.tan(Phi))) * sph_harm( m_sph, N_Coef, Theta, Phi @@ -225,11 +220,9 @@ def Der_Phi_Basis_Fn(M_Coef, N_Coef, Theta, Phi): # M_Coef < 0 Corresonds to Z^ def Der_Phi_Phi_Basis_Fn( M_Coef, N_Coef, Theta, Phi ): # M_Coef < 0 Corresonds to Z^|M|_N - Der_Phi_Phi_Val = 0 if M_Coef < 0: - m_sph = -1 * M_Coef Der_Phi_Phi_Val += ( m_sph @@ -261,7 +254,6 @@ def Der_Phi_Phi_Basis_Fn( ) else: # M_Coef >= 0 - m_sph = M_Coef Der_Phi_Phi_Val += ( m_sph @@ -297,13 +289,11 @@ def Der_Phi_Phi_Basis_Fn( def Lbdv_Cart_To_Sph( Cart_Pts_Wts, ): # takes matrix with rows [x,y,z,w] -> [theta, phi, w] (r=1) - num_lbdv_pts = np.shape(Cart_Pts_Wts)[0] Sph_Pts_Wts = np.zeros((num_lbdv_pts, 3)) for pt in range(num_lbdv_pts): - x = Cart_Pts_Wts[pt][0] y = Cart_Pts_Wts[pt][1] z = Cart_Pts_Wts[pt][2] @@ -321,7 +311,6 @@ def Lbdv_Cart_To_Sph( # Get max degree quad pts for plotting/interpolation: def get_5810_quad_pts(): - euc_quad_pts = Lebedev(5810) sph_quad_pts = Lbdv_Cart_To_Sph(euc_quad_pts) @@ -338,7 +327,6 @@ class lbdv_info(object): # Generates (ONCE) and stores Lebedev Info def __init__( self, Max_SPH_Deg, Num_Quad_Pts ): # This generates lebedev and quad pt info upon instantiation - MY_DIR = os.path.realpath(os.path.dirname(__file__)) # current folder PICKLE_DIR = os.path.join( MY_DIR, "Pickled_LBDV_Files" @@ -378,7 +366,6 @@ def __init__( if os.path.isfile( LBDV_Basis_at_Quad_Pts_Mats_filepath ): # If already pickled, we load it, and split it into the needed arrays: - # print("\n"+"Loading Pickled LBDV data Mats"+"\n") # # BJG: only notify if NEW mats are needed (time-consuming) @@ -400,7 +387,6 @@ def __init__( ] else: # If not pickled, we generate and pickle these files: - print("\n" + "Pickling This LBDV data:" + "\n") ### Generate EVAL_SPH Vals At Quad Pts ONCE ######## @@ -418,7 +404,6 @@ def __init__( for N_Coef in range(Max_SPH_Deg + 1): # 0,...,Max_SPH_Deg for M_Coef in range(-1 * N_Coef, N_Coef + 1): # -N,...,N for quad_pt in range(self.lbdv_quad_pts): - Theta_Quad_Pt = self.Lbdv_Sph_Pts_Quad[quad_pt][0] Phi_Quad_Pt = self.Lbdv_Sph_Pts_Quad[quad_pt][1] Weight_Quad_Pt = self.Lbdv_Sph_Pts_Quad[quad_pt][2] @@ -472,13 +457,11 @@ def __init__( for N_Coef in range(Max_SPH_Deg + 1): # 0,...,Max_SPH_Deg for M_Coef in range(-1 * N_Coef, N_Coef + 1): # -N,...,N for quad_pt in range(self.lbdv_quad_pts): - Theta_Quad_Pt = self.Lbdv_Sph_Pts_Quad[quad_pt][0] Phi_Quad_Pt = self.Lbdv_Sph_Pts_Quad[quad_pt][1] # Dont need Weight, since Quadrature covers that if M_Coef == 0: # We Dont need eta_A for first der in this case - self.SPH_Phi_Der_At_Quad_Pts[N_Coef - M_Coef][N_Coef][ quad_pt ] = Der_Phi_Basis_Fn( @@ -496,7 +479,6 @@ def __init__( ) elif M_Coef >= 0: - self.SPH_Phi_Der_At_Quad_Pts[N_Coef - M_Coef][N_Coef][ quad_pt ] = eta_A( @@ -571,7 +553,6 @@ def __init__( if os.path.isfile( LBDV_Chart_of_Quad_Pts_Mats_filepath ): # If already pickled, we load it, and split it into the needed arrays: - # print("\n"+"Loading Pickled LBDV Chart Mats"+"\n") # # BJG: only notify if NEW mats are needed (time-consuming) @@ -588,7 +569,6 @@ def __init__( ) = np.hsplit(Pickled_LBDV_Charts_Quad_Pts_Mats, 3) else: # If not pickled, we generate and pickle these files: - print("\n" + "Pickling This LBDV Chart data:" + "\n") ### Generate Chart LBDV Data ONCE ######## @@ -608,13 +588,11 @@ def __init__( # Choose which values we use and where: for quad_pt in range(self.lbdv_quad_pts): - theta_pt = self.theta_pts[quad_pt] phi_pt = self.phi_pts[quad_pt] # Determine which Chart Each Quad Pt is in: if Domain(theta_pt, phi_pt) >= 0: - self.Chart_of_Quad_Pts[quad_pt] = 1 else: @@ -629,7 +607,6 @@ def __init__( # Find Rotated Quad Pt at same location: for quad_pt_rot in range(self.lbdv_quad_pts): - theta_bar_pt_rot = self.theta_pts[quad_pt_rot] phi_bar_pt_rot = self.phi_pts[quad_pt_rot] @@ -641,7 +618,6 @@ def __init__( if abs(y_pt - y_pt_rot) < 1e-7: if abs(z_pt - z_pt_rot) < 1e-7: if rot_pt_found is False: - rot_pt_found = True self.Rot_Lbdv_Quad_vals[quad_pt] = quad_pt_rot @@ -697,7 +673,6 @@ def Eval_SPH_Der_Phi_Phi_At_Quad_Pts(self, M, N, Quad_Pt): # Fn that will retrieve all values to be used in quadrature, # in nice format (of All quad points, for Y^M_N) def Eval_SPH_Basis_Wt_M_N(self, M, N): - if M >= 0: return self.SPH_Basis_Wt_At_Quad_Pts[N - M, N, :] else: # If M<0 @@ -705,17 +680,14 @@ def Eval_SPH_Basis_Wt_M_N(self, M, N): # Fn that will retrieve all basis vals at quad pt for product projections def Eval_SPH_At_Quad_Pt_Mat(self, Quad_Pt): - return self.SPH_Basis_At_Quad_Pts[:, :, Quad_Pt] # Fn that will retrieve ALL Der Phi values in nice format def Eval_SPH_Der_Phi_At_Quad_Pt_Mat(self, Quad_Pt): - return self.SPH_Phi_Der_At_Quad_Pts[:, :, Quad_Pt] # Fn that will retrieve ALL 2nd Der Phi values in nice format def Eval_SPH_Der_Phi_Phi_At_Quad_Pt_Mat(self, Quad_Pt): - return self.SPH_Phi_Phi_Der_At_Quad_Pts[:, :, Quad_Pt] # Fn that will retrive vals of Rot_Lbdv_Quad_vals, USE astpye(int)!! @@ -732,7 +704,6 @@ def Eval_Chart_of_Quad_Pts(self, Quad_Pt): # For Splitting Integration between charts: def eta_z(self, quad_pt): - z_pt = self.Z[quad_pt, 0] if abs(z_pt) <= 0.25: diff --git a/src/napari_stress/_stress/lebedev_write_SPB.py b/src/napari_stress/_stress/lebedev_write_SPB.py index 2dae4c66..00aab0cf 100644 --- a/src/napari_stress/_stress/lebedev_write_SPB.py +++ b/src/napari_stress/_stress/lebedev_write_SPB.py @@ -323,7 +323,6 @@ def leb194(): def leb230(): - v = -0.5522639919727325e-1 out = genOh_a00(v) v = 0.4450274607445226e-2 @@ -436,7 +435,6 @@ def leb302(): def leb350(): - v = 0.3006796749453936e-2 out = genOh_a00(v) v = 0.3050627745650771e-2 @@ -481,7 +479,6 @@ def leb350(): def leb434(): - v = 0.5265897968224436e-3 out = genOh_a00(v) v = 0.2548219972002607e-2 @@ -535,7 +532,6 @@ def leb434(): def leb590(): - v = 0.3095121295306187e-3 out = genOh_a00(v) v = 0.1852379698597489e-2 @@ -2915,7 +2911,6 @@ def leb3890(): def leb4334(): - v = 0.1449063022537883e-4 out = genOh_a00(v) v = 0.2546377329828424e-3 @@ -3320,7 +3315,6 @@ def leb4334(): def leb4802(): - v = 0.9687521879420705e-4 out = genOh_a00(v) v = 0.2307897895367918e-3 @@ -3766,7 +3760,6 @@ def leb4802(): def leb5294(): - v = 0.9080510764308163e-4 out = genOh_a00(v) v = 0.2084824361987793e-3 @@ -4255,7 +4248,6 @@ def leb5294(): def leb5810(): - v = 0.9735347946175486e-5 out = genOh_a00(v) v = 0.1907581241803167e-3 diff --git a/src/napari_stress/_stress/manifold_SPB.py b/src/napari_stress/_stress/manifold_SPB.py index 3616a871..0841da12 100644 --- a/src/napari_stress/_stress/manifold_SPB.py +++ b/src/napari_stress/_stress/manifold_SPB.py @@ -19,7 +19,6 @@ # Use PREVIOUS definitions below to define x,y,z maps for SPB coors: def Manny_Fn_Def(theta, phi, r_0, Manny_Name, IsRadial): - x_val = Manny_Fn_Def_X(theta, phi, r_0, Manny_Name, IsRadial) y_val = Manny_Fn_Def_Y(theta, phi, r_0, Manny_Name, IsRadial) z_val = Manny_Fn_Def_Z(theta, phi, r_0, Manny_Name, IsRadial) @@ -29,7 +28,6 @@ def Manny_Fn_Def(theta, phi, r_0, Manny_Name, IsRadial): # We may want seperate defs for each coordinate: def Manny_Fn_Def_X(theta, phi, r_0, Manny_Name, IsRadial): - if IsRadial is True: return ( Radial_Manifold_R_Def(theta, phi, r_0, Manny_Name) @@ -41,7 +39,6 @@ def Manny_Fn_Def_X(theta, phi, r_0, Manny_Name, IsRadial): def Manny_Fn_Def_Y(theta, phi, r_0, Manny_Name, IsRadial): - if IsRadial is True: return ( Radial_Manifold_R_Def(theta, phi, r_0, Manny_Name) @@ -53,7 +50,6 @@ def Manny_Fn_Def_Y(theta, phi, r_0, Manny_Name, IsRadial): def Manny_Fn_Def_Z(theta, phi, r_0, Manny_Name, IsRadial): - if IsRadial is True: return Radial_Manifold_R_Def(theta, phi, r_0, Manny_Name) * np.cos(phi) else: @@ -80,7 +76,6 @@ def Non_Radial_Manifold_X_Def(theta, phi, r_0, Manny_Name): elif Manny_Name == "Cone_Head_Logis": return np.cos(theta) * np.sin(phi) elif Manny_Name == "Fission_Yeast_R2" or Manny_Name == "Fission_Yeast_R1pt2": - R = 0.0 if Manny_Name == "Fission_Yeast_R2": R = 2.0 # fix for this geo @@ -143,7 +138,6 @@ def Non_Radial_Manifold_Y_Def(theta, phi, r_0, Manny_Name): elif Manny_Name == "Cone_Head_Logis": return np.sin(theta) * np.sin(phi) elif Manny_Name == "Fission_Yeast_R2" or Manny_Name == "Fission_Yeast_R1pt2": - R = 0.0 if Manny_Name == "Fission_Yeast_R2": R = 2.0 # fix for this geo @@ -204,7 +198,6 @@ def Non_Radial_Manifold_Z_Def(theta, phi, r_0, Manny_Name): elif Manny_Name == "Cone_Head_Logis": return np.cos(phi) * (1.0 + r_0 / (1.0 + np.exp(-20.0 * (np.pi / 16 - phi)))) elif Manny_Name == "Fission_Yeast_R2" or Manny_Name == "Fission_Yeast_R1pt2": - R = 0.0 if Manny_Name == "Fission_Yeast_R2": R = 2.0 # fix for this geo @@ -269,7 +262,6 @@ def Man_Name_Radial(Man_Name): # Define RADIAL Manifolds ONCE here: def Radial_Manifold_R_Def(theta, phi, r_0, Manny_Name): - if Manny_Name == "S2": # Unit Sphere return 1.0 @@ -475,7 +467,6 @@ def Radial_Manifold_R_Def(theta, phi, r_0, Manny_Name): elif ( Manny_Name == "Fission_Yeast_Rad_R2" or Manny_Name == "Fission_Yeast_Rad_R1pt2" ): - R = 0.0 if Manny_Name == "Fission_Yeast_Rad_R2": R = 2.0 @@ -511,7 +502,6 @@ def Radial_Manifold_R_Def(theta, phi, r_0, Manny_Name): ) elif Manny_Name == "Divet": - if np.isscalar(phi): if np.cos(phi) > 0.9: return 1.0 - r_0 * np.exp( @@ -530,7 +520,6 @@ def Radial_Manifold_R_Def(theta, phi, r_0, Manny_Name): ) elif Manny_Name == "Double_Divet": - if np.isscalar(phi): if np.cos(phi) > 0.9: return 1.0 - r_0 * np.exp( @@ -564,7 +553,6 @@ def Radial_Manifold_R_Def(theta, phi, r_0, Manny_Name): ) elif Manny_Name == "Divets_Around_Pimple": - if np.isscalar(phi): if abs(np.cos(phi)) > 0.9: return 1.0 - r_0 * np.exp( @@ -620,11 +608,9 @@ def Radial_Manifold_R_Def(theta, phi, r_0, Manny_Name): # Takes an array of R0 Values and returns string format # for Rayleigh Dissipation Studies (!!!!VTU NOT PICKLE!!!!) def Max_Decimal_R0_Array(R0_Values_Array): - max_dec_length = 0 for R0_i in range(len(R0_Values_Array)): - R0_Value_i = R0_Values_Array[R0_i] str_index = str(R0_Value_i)[::-1].find(".") @@ -634,7 +620,6 @@ def Max_Decimal_R0_Array(R0_Values_Array): new_R0_array_labels = [] for R0_j in range(len(R0_Values_Array)): - a, b = str(R0_Values_Array[R0_j]).split(".", 1) b_adjusted = b.ljust(max_dec_length, "0") @@ -696,7 +681,6 @@ def __init__( # BJG: given name and r0 value, we can use function to get points we need: if self.Use_Man_Name is True: - self.IsManRad = self.Man_Shape_Dict["Is_Manifold_Radial"] self.Man_R0_Val = self.Man_Shape_Dict["Maniold_R0_Value"] self.Man_Shape_Name = self.Man_Shape_Dict[ @@ -1513,7 +1497,6 @@ def zero_vector_of_basis_mats(): or self.Man_Official_Name == [] or self.pickling is False ): # If we need to (re)generate these: - self.rho_A_Mats = zero_vector_of_basis_mats() # rho = G*(A^-1) self.rho_B_Mats = zero_vector_of_basis_mats() @@ -1571,7 +1554,6 @@ def zero_vector_of_basis_mats(): for quad_pt in range(self.num_quad_pts): if lbdv.Chart_of_Quad_Pts[quad_pt] > 0: - # Change of Basis, metric tensor mats, and their first derivatives A_Mat_pt = self.Change_Basis_Mat(quad_pt, "A") B_Mat_pt = self.Change_Basis_Mat(quad_pt, "B") @@ -1643,7 +1625,6 @@ def zero_vector_of_basis_mats(): and self.Man_Official_Name != [] and self.pickling is True ): - print("pickling Manny Inv Mats for re-use" + "\n") # We save matricies as a list: @@ -1678,7 +1659,6 @@ def zero_vector_of_basis_mats(): # If we have already pickled the above matricies, we load them: else: - print("\n" + "loading pickled Manny Inv Mats" + "\n") # print("Pickled_Inverse_Mats.shape = "+str(Pickled_Inverse_Mats.shape)) @@ -1876,7 +1856,6 @@ def W_22(self, quad_pt, Chart): # First Fundamental Form def I_Mat(self, quad_pt, Chart): - return np.array( [ [self.E(quad_pt, Chart), self.F(quad_pt, Chart)], @@ -1886,7 +1865,6 @@ def I_Mat(self, quad_pt, Chart): # First Fundamental Form def I_theta_Mat(self, quad_pt, Chart): - return np.array( [ [self.E_theta(quad_pt, Chart), self.F_theta(quad_pt, Chart)], @@ -1896,7 +1874,6 @@ def I_theta_Mat(self, quad_pt, Chart): # First Fundamental Form def I_phi_Mat(self, quad_pt, Chart): - return np.array( [ [self.E_phi(quad_pt, Chart), self.F_phi(quad_pt, Chart)], @@ -1908,7 +1885,6 @@ def I_phi_Mat(self, quad_pt, Chart): # Weingarten Map, derivative of (inward) normal vector def Wein_Map(self, quad_pt, Chart): - return np.array( [ [self.W_11(quad_pt, Chart), self.W_12(quad_pt, Chart)], @@ -1918,7 +1894,6 @@ def Wein_Map(self, quad_pt, Chart): # Derivative of inward normal, wrt e_theta def Normal_vec_dtheta(self, quad_pt, Chart): - W_map = self.Wein_Map(quad_pt, Chart) sigma_theta_C_pt = self.sigma_theta(quad_pt, Chart) sigma_phi_C_pt = self.sigma_phi(quad_pt, Chart) @@ -1939,7 +1914,6 @@ def Normal_vec_dtheta(self, quad_pt, Chart): # Derivative of inward normal, wrt e_phi def Normal_vec_dphi(self, quad_pt, Chart): - W_map = self.Wein_Map(quad_pt, Chart) sigma_theta_C_pt = self.sigma_theta(quad_pt, Chart) sigma_phi_C_pt = self.sigma_phi(quad_pt, Chart) @@ -1960,7 +1934,6 @@ def Normal_vec_dphi(self, quad_pt, Chart): # Matrix Used to Convert to Polar def Change_Basis_Mat(self, quad_pt, Chart): - Basis_mat = np.zeros((3, 3)) Basis_mat[:, 0] = self.sigma_theta(quad_pt, Chart).T @@ -1971,7 +1944,6 @@ def Change_Basis_Mat(self, quad_pt, Chart): # Entry-wise deriv of Basis mat, dtheta def dChange_Basis_Mat_theta(self, quad_pt, Chart): - dBasis_mat_theta = np.zeros((3, 3)) dBasis_mat_theta[:, 0] = self.sigma_theta_theta(quad_pt, Chart).T @@ -1982,7 +1954,6 @@ def dChange_Basis_Mat_theta(self, quad_pt, Chart): # Entry-wise deriv of Basis mat, dphi def dChange_Basis_Mat_phi(self, quad_pt, Chart): - dBasis_mat_phi = np.zeros((3, 3)) dBasis_mat_phi[:, 0] = self.sigma_theta_phi(quad_pt, Chart).T @@ -1993,21 +1964,18 @@ def dChange_Basis_Mat_phi(self, quad_pt, Chart): # I mat is embedded: def G_Mat(self, quad_pt, Chart): - G_mat = np.zeros((3, 3)) G_mat[0:2, 0:2] = self.I_Mat(quad_pt, Chart) return G_mat def dG_Mat_theta(self, quad_pt, Chart): - G_theta_mat = np.zeros((3, 3)) G_theta_mat[0:2, 0:2] = self.I_theta_Mat(quad_pt, Chart) return G_theta_mat def dG_Mat_phi(self, quad_pt, Chart): - G_phi_mat = np.zeros((3, 3)) G_phi_mat[0:2, 0:2] = self.I_phi_Mat(quad_pt, Chart) diff --git a/src/napari_stress/_stress/sph_func_SPB.py b/src/napari_stress/_stress/sph_func_SPB.py index dcdd3001..8ec32776 100644 --- a/src/napari_stress/_stress/sph_func_SPB.py +++ b/src/napari_stress/_stress/sph_func_SPB.py @@ -20,7 +20,6 @@ def Mass_Mat_Exact(m_coef, n_coef): # Inverse Mass Matrix on S2 pullback (NOT INCLUDING constant mode def Inv_Mass_Mat_on_Pullback(SPH_Deg): - mat_dim = (SPH_Deg + 1) ** 2 - 1 inv_mass_mat = np.zeros((mat_dim, mat_dim)) @@ -28,7 +27,6 @@ def Inv_Mass_Mat_on_Pullback(SPH_Deg): for n in range(1, SPH_Deg + 1): for m in range(-1 * n, n + 1): - inv_mass_mat[diag_entry, diag_entry] = 1.0 / Mass_Mat_Exact(m, n) return inv_mass_mat @@ -48,7 +46,6 @@ def Create_Basis_Fn(m, n, basis_degree): # Approximate func(theta, phi) in FE space def Proj_Func(func, SPH_Deg, lbdv): - Proj_Coef = np.zeros([SPH_Deg + 1, SPH_Deg + 1]) theta_quad_pts = lbdv.Lbdv_Sph_Pts_Quad[:, 0] @@ -64,7 +61,6 @@ def Proj_Func(func, SPH_Deg, lbdv): # Compute inner product of theta der with each basis elt for n in range(SPH_Deg + 1): for m in range(-1 * n, n + 1): - Proj_Coef_mn = sum( np.multiply(lbdv.Eval_SPH_Basis_Wt_M_N(m, n), func_vals_at_quad_pts) ) / Mass_Mat_Exact(m, n) @@ -80,14 +76,12 @@ def Proj_Func(func, SPH_Deg, lbdv): # Projects function into both charts def Proj_Into_SPH_Charts(func, Coef_Deg, lbdv): - # In chart A # Coef_A = Proj_Func(lambda theta, phi: eta_A(func, theta, phi), Coef_Deg) sph_func_A = Proj_Func(func, Coef_Deg, lbdv) # using class structure # rotate func to Chart B Coors def func_rot(theta_bar, phi_bar): - A_Coors = Coor_B_To_A(theta_bar, phi_bar) Theta_Vals = A_Coors[0] @@ -104,7 +98,6 @@ def func_rot(theta_bar, phi_bar): # Projects function into both charts using quad_pts def Proj_Into_SPH_Charts_At_Quad_Pts(func_quad_vals, Proj_Deg, lbdv): - # In chart A sph_func_A = Faster_Double_Proj( func_quad_vals, Proj_Deg, lbdv @@ -116,7 +109,6 @@ def Proj_Into_SPH_Charts_At_Quad_Pts(func_quad_vals, Proj_Deg, lbdv): # rotate func to Chart B Coors for quad_pt in range(lbdv.lbdv_quad_pts): - quad_pt_vals_B[quad_pt] = func_quad_vals[ lbdv.Eval_Inv_Rot_Lbdv_Quad_vals(quad_pt) ] @@ -130,7 +122,6 @@ def Proj_Into_SPH_Charts_At_Quad_Pts(func_quad_vals, Proj_Deg, lbdv): # TAKES VALS AT QUAD PTS, to Approximate func(theta, phi) in SPH Basis def Faster_Double_Proj(func_quad_vals, Proj_Deg, lbdv): - Proj_Coef = np.zeros( [Proj_Deg + 1, Proj_Deg + 1] ) # Size of basis used to represent derivative @@ -138,7 +129,6 @@ def Faster_Double_Proj(func_quad_vals, Proj_Deg, lbdv): # Compute inner product of theta der with each basis elt for n in range(Proj_Deg + 1): for m in range(-1 * n, n + 1): - I_mn = 0 for quad_pt in range(lbdv.lbdv_quad_pts): @@ -165,7 +155,6 @@ def Faster_Double_Proj(func_quad_vals, Proj_Deg, lbdv): # TAKES VALS AT QUAD PTS, to Approximate # func(theta, phi)*Coef_Mat(theta, phi) in SPH Basis def Faster_Double_Proj_Product(func1_quad_vals, func2_quad_vals, Proj_Deg, lbdv): - Proj_Product_Coef = np.zeros( [Proj_Deg + 1, Proj_Deg + 1] ) # Size of basis used to represent derivative @@ -173,7 +162,6 @@ def Faster_Double_Proj_Product(func1_quad_vals, func2_quad_vals, Proj_Deg, lbdv) # Compute inner product of theta der with each basis elt for n in range(Proj_Deg + 1): for m in range(-1 * n, n + 1): - I_mn = 0 for quad_pt in range(lbdv.lbdv_quad_pts): @@ -201,12 +189,10 @@ def Faster_Double_Proj_Product(func1_quad_vals, func2_quad_vals, Proj_Deg, lbdv) # Inputs quad vals for f, f_approx, integrates on SPHERE: def Lp_Rel_Error_At_Quad(approx_f_vals, f_vals, lbdv, p): # Assumes f NOT 0 - Lp_Err = 0 # ||self - f||_p Lp_f = 0 # || f ||_p for quad_pt in range(lbdv.lbdv_quad_pts): - weight_pt = lbdv.Lbdv_Sph_Pts_Quad[quad_pt][2] Lp_Err_Pt = abs((approx_f_vals[quad_pt] - f_vals[quad_pt]) ** p) * weight_pt @@ -223,12 +209,10 @@ def Lp_Rel_Error_At_Quad(approx_f_vals, f_vals, lbdv, p): # Assumes f NOT 0 # Inputs quad vals for f, f_approx, integrates on CHART: def Lp_Rel_Error_At_Quad_In_Chart(approx_f_vals, f_vals, lbdv, p): # Assumes f NOT 0 - Lp_Err = 0 # ||self - f||_p Lp_f = 0 # || f ||_p for quad_pt in range(lbdv.lbdv_quad_pts): - if lbdv.Chart_of_Quad_Pts[quad_pt] > 0: weight_pt = lbdv.Lbdv_Sph_Pts_Quad[quad_pt][2] @@ -243,7 +227,6 @@ def Lp_Rel_Error_At_Quad_In_Chart(approx_f_vals, f_vals, lbdv, p): # Assumes f # Computes Coef of constant mode of Fn def Const_SPH_Mode_Of_Func(func, lbdv): - theta_quad_pts = lbdv.Lbdv_Sph_Pts_Quad[:, 0] phi_quad_pts = lbdv.Lbdv_Sph_Pts_Quad[:, 1] @@ -258,7 +241,6 @@ def Const_SPH_Mode_Of_Func(func, lbdv): # Computes Average of sph proj of quad pts def Avg_of_SPH_Proj_of_Func(func_vals_at_quad_pts, lbdv): - Proj_Coef_Const_Mode = sum( np.multiply(lbdv.Eval_SPH_Basis_Wt_M_N(0, 0), func_vals_at_quad_pts.flatten()) ) / Mass_Mat_Exact(0, 0) @@ -280,7 +262,6 @@ def Un_Flatten_Coef_Vec(coefficient_vector, basis_degree): row = 0 for n in range(basis_degree + 1): for m in range(-1 * n, n + 1): - if m > 0: coefficient_matrix[n - m][n] = coefficient_vector[row] @@ -350,11 +331,9 @@ def convert_coeffcients_stress_to_pyshtools( # Gives L_1 Integral on SPHERE pullback: def L1_Integral(f_quad_vals, lbdv): - L1_Int = 0 for quad_pt in range(lbdv.lbdv_quad_pts): - weight_pt = lbdv.Lbdv_Sph_Pts_Quad[quad_pt][2] L1_Int_Pt = abs(f_quad_vals[quad_pt]) * weight_pt @@ -365,11 +344,9 @@ def L1_Integral(f_quad_vals, lbdv): # Gives Integral on SPHERE pullback: def S2_Integral(f_quad_vals, lbdv): - S2_Int = 0 for quad_pt in range(lbdv.lbdv_quad_pts): - weight_pt = lbdv.Lbdv_Sph_Pts_Quad[quad_pt][2] S2_Int_Pt = f_quad_vals[quad_pt, 0] * weight_pt @@ -395,12 +372,10 @@ def __init__(self, SPH_Coef, SPH_Deg): # Evaluates the Coef Matrix Representing SPH Basis of a Fn def Eval_SPH(self, Theta, Phi): - SPH_Val = 0 for L_Coef in range(0, self.sph_deg + 1): for M_Coef in range(0, L_Coef + 1): - SPH_Val += self.sph_coef[L_Coef - M_Coef][L_Coef] * Eval_SPH_Basis( M_Coef, L_Coef, Theta, Phi ) @@ -414,12 +389,10 @@ def Eval_SPH(self, Theta, Phi): # Evaluates the Phi Derivatibve of Coef Matrix Representing SPH Basis of a Fn def Eval_SPH_Der_Phi(self, Theta, Phi): - Der_SPH_Val = 0 for L_Coef in range(0, self.sph_deg + 1): for M_Coef in range(0, L_Coef + 1): - Der_SPH_Val += self.sph_coef[L_Coef - M_Coef][ L_Coef ] * Der_Phi_Basis_Fn(M_Coef, L_Coef, Theta, Phi) @@ -433,12 +406,10 @@ def Eval_SPH_Der_Phi(self, Theta, Phi): # Evaluates the Second Phi Derivatibve of Coef Matrix Representing SPH Basis of a Fn def Eval_SPH_Der_Phi_Phi(self, Theta, Phi): - Sec_Der_SPH_Val = 0 for L_Coef in range(0, self.sph_deg + 1): for M_Coef in range(0, L_Coef + 1): - Sec_Der_SPH_Val += self.sph_coef[L_Coef - M_Coef][ L_Coef ] * Der_Phi_Phi_Basis_Fn(M_Coef, L_Coef, Theta, Phi) @@ -453,7 +424,6 @@ def Eval_SPH_Der_Phi_Phi(self, Theta, Phi): # Should Evalaute Phi Der of Matrix at Quad_Pt (times weight), # to be used in quadrature, def Eval_SPH_Der_Phi_Coef(self, Quad_Pt, lbdv): - # For Scalar Case, we use usual vectorization: if np.isscalar(Quad_Pt): return sum( @@ -474,7 +444,6 @@ def Eval_SPH_Der_Phi_Coef(self, Quad_Pt, lbdv): # Should Evalaute 2nd Phi Der of Matrix at Quad_Pt (times weight), # to be used in quadrature, def Eval_SPH_Der_Phi_Phi_Coef(self, Quad_Pt, lbdv): - # For Scalar Case, we use usual vectorization: if np.isscalar(Quad_Pt): return sum( @@ -494,10 +463,8 @@ def Eval_SPH_Der_Phi_Phi_Coef(self, Quad_Pt, lbdv): # Evaluates A SPH fn at a quad pt(s) def Eval_SPH_Coef_Mat(self, Quad_Pt, lbdv): - # For Scalar Case, we use usual vectorization: if np.isscalar(Quad_Pt): - return sum( np.multiply(self.sph_coef, lbdv.Eval_SPH_At_Quad_Pt_Mat(Quad_Pt)) ) @@ -505,19 +472,16 @@ def Eval_SPH_Coef_Mat(self, Quad_Pt, lbdv): # For multiple quad pts, we use einstein sumation to # output vector of solutions at each point: else: - return np.einsum( "ij, ijk -> k", self.sph_coef, lbdv.Eval_SPH_At_Quad_Pt_Mat(Quad_Pt) ).reshape((len(Quad_Pt), 1)) # Use the fact that Theta Der Formula is EXACT for our basis def Quick_Theta_Der(self): - Der_Coef = np.zeros(np.shape(self.sph_coef)) for n_coef in range(self.sph_deg + 1): for m_coef in range(1, n_coef + 1): # Theta Der of Y^0_n = 0 - # D_theta X^m_n = -m*Z^m_n Der_Coef[n_coef - m_coef][n_coef] = ( m_coef * self.sph_coef[n_coef][n_coef - m_coef] @@ -535,7 +499,6 @@ def Quick_Theta_Bar_Der(self): # For rotated coordinate frame # Approximate func(theta, phi)*Coef_Mat(theta, phi) in SPH Basis def Fast_Proj_Product(self, func, Proj_Deg, lbdv): - Proj_Product_Coef = np.zeros( [Proj_Deg + 1, Proj_Deg + 1] ) # Size of basis used to represent derivative @@ -543,7 +506,6 @@ def Fast_Proj_Product(self, func, Proj_Deg, lbdv): # Compute inner product of theta der with each basis elt for n in range(Proj_Deg + 1): for m in range(-1 * n, n + 1): - I_mn = 0 for quad_pt in range(lbdv.lbdv_quad_pts): @@ -571,7 +533,6 @@ def Fast_Proj_Product(self, func, Proj_Deg, lbdv): # TAKES VALS AT QUAD PTS, to Approximate # func(theta, phi)*Coef_Mat(theta, phi) in SPH Basis def Faster_Proj_Product(self, func_quad_vals, Proj_Deg, lbdv): - Proj_Product_Coef = np.zeros( [Proj_Deg + 1, Proj_Deg + 1] ) # Size of basis used to represent derivative @@ -579,7 +540,6 @@ def Faster_Proj_Product(self, func_quad_vals, Proj_Deg, lbdv): # Compute inner product of theta der with each basis elt for n in range(Proj_Deg + 1): for m in range(-1 * n, n + 1): - I_mn = 0 for quad_pt in range(lbdv.lbdv_quad_pts): @@ -606,7 +566,6 @@ def Faster_Proj_Product(self, func_quad_vals, Proj_Deg, lbdv): # Approximate Coef_Mat2(theta, phi)*Coef_Mat(theta, phi) in SPH Basis def Fast_Proj_Product_SPH(self, SPH_Func2, Proj_Deg, lbdv): - Proj_Product_Coef = np.zeros( [Proj_Deg + 1, Proj_Deg + 1] ) # Size of basis used to represent derivative @@ -614,7 +573,6 @@ def Fast_Proj_Product_SPH(self, SPH_Func2, Proj_Deg, lbdv): # Compute inner product of theta der with each basis elt for n in range(Proj_Deg + 1): for m in range(-1 * n, n + 1): - I_mn = 0 for quad_pt in range(lbdv.lbdv_quad_pts): @@ -640,7 +598,6 @@ def Fast_Proj_Product_SPH(self, SPH_Func2, Proj_Deg, lbdv): return spherical_harmonics_function(Proj_Product_Coef, Proj_Deg) def Inner_Product_SPH(self, Other_SPH): - Vec1 = self.sph_coef.flatten() Vec2 = Other_SPH.sph_coef.flatten() @@ -660,7 +617,6 @@ def L2_Norm_SPH(self): return np.sqrt(self.Inner_Product_SPH(self)) def Lp_Rel_Error_in_Chart(self, f, lbdv, p): # Assumes f NOT 0 - Lp_Err = 0 # ||self - f||_p Lp_f = 0 # || f ||_p @@ -683,7 +639,6 @@ def Lp_Rel_Error_in_Chart(self, f, lbdv, p): # Assumes f NOT 0 # Looks at Rel Error in ALL of S^2 def Lp_Rel_Error_in_S2(self, f, lbdv, p): # Assumes f NOT 0 - Lp_Err = 0 # ||self - f||_p Lp_f = 0 # || f ||_p @@ -705,7 +660,6 @@ def Lp_Rel_Error_in_S2(self, f, lbdv, p): # Assumes f NOT 0 # Returns properly ordered vector of SPH Coef def Flatten_Coef_Mat(self): - coef_mat = self.sph_coef sph_degree = self.sph_deg @@ -714,7 +668,6 @@ def Flatten_Coef_Mat(self): for n in range(sph_degree + 1): for m in range(-1 * n, n + 1): - if m > 0: coef_vec[row, 0] = coef_mat[n - m][n] diff --git a/src/napari_stress/_surface.py b/src/napari_stress/_surface.py index 97783b3e..b4b4f53f 100644 --- a/src/napari_stress/_surface.py +++ b/src/napari_stress/_surface.py @@ -143,7 +143,6 @@ def smooth_sinc( feature_angle: float = 60, boundary: bool = False, ) -> SurfaceData: - mesh = vedo.mesh.Mesh((surface[0], surface[1])) mesh.smooth( niter=niter, @@ -160,7 +159,6 @@ def smooth_sinc( def smoothMLS2D( points: PointsData, factor: float = 0.5, radius: float = None ) -> PointsData: - pointcloud = vedo.pointcloud.Points(points) pointcloud.smoothMLS2D(f=factor, radius=radius) @@ -173,7 +171,6 @@ def smoothMLS2D( @register_function(menu="Surfaces > Simplify (decimate, vedo, n-STRESS)") @frame_by_frame def decimate(surface: SurfaceData, fraction: float = 0.1) -> SurfaceData: - mesh = vedo.mesh.Mesh((surface[0], surface[1])) n_vertices = mesh.N() diff --git a/src/napari_stress/_tests/test_spherical_harmonic_fit.py b/src/napari_stress/_tests/test_spherical_harmonic_fit.py index badce9dd..bb3d300c 100644 --- a/src/napari_stress/_tests/test_spherical_harmonic_fit.py +++ b/src/napari_stress/_tests/test_spherical_harmonic_fit.py @@ -90,7 +90,6 @@ def test_spherical_harmonics(): def test_interoperatibility(): - from napari_stress._spherical_harmonics import spherical_harmonics as sh from napari_stress._stress.sph_func_SPB import ( convert_coeffcients_stress_to_pyshtools, @@ -120,7 +119,6 @@ def test_interoperatibility(): def test_lebedev_points(): - from napari_stress._stress.lebedev_write_SPB import LebFunc for i in [ diff --git a/src/napari_stress/_tests/test_types.py b/src/napari_stress/_tests/test_types.py index 18d0d292..b88e096f 100644 --- a/src/napari_stress/_tests/test_types.py +++ b/src/napari_stress/_tests/test_types.py @@ -15,7 +15,6 @@ def test_custom_types(make_napari_viewer): ) def test_function(argument: manifold) -> manifold: - layer = None if isinstance(argument, layers.Layer): layer = argument diff --git a/src/napari_stress/_tests/test_utils.py b/src/napari_stress/_tests/test_utils.py index f53b278d..8b97bdf2 100644 --- a/src/napari_stress/_tests/test_utils.py +++ b/src/napari_stress/_tests/test_utils.py @@ -94,7 +94,6 @@ def test_decorator_points_layerdatatuple(): def test_decorator_points_layers(): - from napari_stress import TimelapseConverter, get_droplet_point_cloud_4d Converter = TimelapseConverter() @@ -148,7 +147,6 @@ def test_decorator_surfaces(): def test_decorator_images(): - from napari_stress import TimelapseConverter Converter = TimelapseConverter() @@ -168,7 +166,6 @@ def test_decorator_images(): def test_frame_by_frame_vectors(): - from napari_stress import TimelapseConverter Converter = TimelapseConverter() diff --git a/src/napari_stress/_utils/coordinate_conversion.py b/src/napari_stress/_utils/coordinate_conversion.py index 9ea5310e..38a9d9ff 100644 --- a/src/napari_stress/_utils/coordinate_conversion.py +++ b/src/napari_stress/_utils/coordinate_conversion.py @@ -85,7 +85,6 @@ def polynomial_to_parameters3D(coefficients: np.ndarray): # Convert input R^3 points into Ellipsoidal coors: def cartesian_to_elliptical(ellipsoid: VectorsData, points: PointsData) -> np.ndarray: - # first, get lengths of ellipsoid axes: lengths = np.linalg.norm(ellipsoid[:, 1], axis=1) @@ -135,7 +134,6 @@ def cartesian_to_elliptical(ellipsoid: VectorsData, points: PointsData) -> np.nd def elliptical_to_cartesian( U_pts_cloud: np.ndarray, V_pts_cloud: np.ndarray, ellipsoid: VectorsData ) -> PointsData: - R_inverse = _orientation_from_ellipsoid(ellipsoid).T lengths = _axes_lengths_from_ellipsoid(ellipsoid) center = _center_from_ellipsoid(ellipsoid) diff --git a/src/napari_stress/_utils/fit_utils.py b/src/napari_stress/_utils/fit_utils.py index c1533c1a..8c4e3742 100644 --- a/src/napari_stress/_utils/fit_utils.py +++ b/src/napari_stress/_utils/fit_utils.py @@ -48,7 +48,6 @@ def _detect_max_gradient(array: np.ndarray, center: float = None): def _function_args_to_list(function: callable) -> list: - sig = inspect.signature(function) return list(sig.parameters.keys()) diff --git a/src/napari_stress/_utils/frame_by_frame.py b/src/napari_stress/_utils/frame_by_frame.py index fc4fd0e2..e1d19f98 100644 --- a/src/napari_stress/_utils/frame_by_frame.py +++ b/src/napari_stress/_utils/frame_by_frame.py @@ -23,7 +23,6 @@ def frame_by_frame(function: callable, progress_bar: bool = False): @wraps(function) def wrapper(*args, **kwargs): - sig = inspect.signature(function) annotations = [sig.parameters[key].annotation for key in sig.parameters.keys()] @@ -98,7 +97,6 @@ class TimelapseConverter: """ def __init__(self): - # Supported LayerData types self.data_to_list_conversion_functions = { PointsData: self._points_to_list_of_points, @@ -479,7 +477,6 @@ def _surface_to_list_of_surfaces(self, surface: SurfaceData) -> list: # as previously determined surfaces = [None] * n_frames for t in range(n_frames): - # Find points with correct frame index _points = points[points[:, 0] == t, 1:] @@ -519,7 +516,6 @@ def _list_of_surfaces_to_surface(self, surfaces: list) -> tuple: n_vertices = 0 for idx, surface in enumerate(surfaces): - # Offset indices in faces list by previous amount of points faces[idx] = n_vertices + np.array(faces[idx]) @@ -567,7 +563,6 @@ def _points_to_list_of_points(self, points: PointsData) -> list: # Fill the respective entries in the list with coordinates in the # original data where time-coordinate matches the current frame for t in range(n_frames): - # Find points with correct frame index points_out[t] = points[points[:, 0] == t, 1:]