fixed bug where -2 direction from points on DB edge returned None; ma…

…de ways to evaluate how different point placement methods look in terms of variation of direction and angle to neighbors

scripts/DistortionGeometry.py

@@ -1,21 +1,11 @@
 # figuring out the correspondence between point code / peel coordinates, xyz, and location projected onto the face plane
+# compare different methods of placing points onto face planes based on their raw peel coordinates (trying to get nice distribution of lattice on the sphere)
 
 import numpy as np
-import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
-from functools import reduce
 
-import icosalattice.StartingPoints as sp
-import icosalattice.IcosahedronMath as icm
-import icosalattice.CoordinatesByAncestry as anc
-import icosalattice.PeelCoordinates as pe
-import icosalattice.PointCodeArithmetic as pca
-import icosalattice.MathUtil as mu
-import icosalattice.GeneratePointCodes as gpc
-import icosalattice.PlotPointLocations as ppl
-from icosalattice.PointPaths import get_point_path, get_stepwise_path_distances_and_angles_2d
-from icosalattice.PlotPaths import plot_distances_and_angles_2d
 from icosalattice.CoordinatesOfPointCode import METHOD_NAME_TO_FUNCTION_POINT_CODE_TO_XYZ
+import icosalattice.EvaluatePointPlacementMethods as epp
+import icosalattice.Adjacency as adj
 
 
 
@@ -37,73 +27,21 @@ def vectors_between_pairs_of_point_codes_are_parallel(pair0, pair1, func_pc_to_x
 
 
 if __name__ == "__main__":
-    # plot a half-peel and some basic points
-    ls = [0, 0, 0, 0.5, 0.5, 0.5, 1, 1, 1]
-    ds = [0, 0.5, 1] * 3
-    labels = ["C", "C3", "L", "C1", "C2", "L1", "A", "K1", "K"]
-    ppl.plot_points_by_peel_coordinates_on_one_half_peel(ls, ds, labels)
-
-    # # test if line segments are parallel on the face plane
-    # reference_pair_l = ["C", "A"]
-    # test_pairs_l = [
-    #     ["C2", "K1"], ["C02", "K11"], ["C22", "K01"], ["C002", "K111"], ["C022", "K101"], ["C202", "K011"], ["C222", "K001"],  # both points on edges
-    #     ["C13", "K11"], ["C2", "C21"],  # one point on an edge and the other inside face
-    #     ["C13", "C12"], ["C012", "C112"],  # both points inside face
-    # ]
-    # reference_pair_dl = ["C", "K"]
-    # test_pairs_dl = [
-    #     ["C1", "K1"], ["C11", "K11"], ["C111", "K111"], ["C01", "K01"], ["C101", "K101"], ["C011", "K011"], ["C001", "K001"],  # both points on edges
-    #     ["C1", "C12"], ["C01", "C13"], ["C12", "K1"],  # one point on an edge and the other inside face
-    #     ["C21", "C212"], ["C13", "C21"],  # both points inside face
-    # ]
-    # reference_pair_d = ["A", "K"]  # this is only in the D direction from the perspective of the CAKL half-peel; from K's perspective it is the R direction!
-    # test_pairs_d = [
-    #     ["C1", "C2"], ["C01", "C02"], ["C11", "C22"], ["C001", "C002"], ["C011", "C022"], ["C101", "C202"], ["C111", "C222"],  # both points on edges
-    #     ["C1", "C13"], ["C12", "C22"],  # one point on an edge and the other inside face
-    #     ["C12", "C21"], ["C0103", "C0133"],  # both points inside face
-    # ]
-    # for reference_pair, test_pairs in zip(
-    #     [reference_pair_l, reference_pair_dl, reference_pair_d],
-    #     [test_pairs_l, test_pairs_dl, test_pairs_d],
-    # ):
-    #     pair0 = reference_pair
-    #     for pair1 in test_pairs:
-    #         if vectors_between_pairs_of_point_codes_are_parallel(pair0, pair1):
-    #             print(f"{pair0} vs {pair1}: IS parallel")
-    #         else:
-    #             print(f"{pair0} vs {pair1}: is NOT parallel")
+    # # plot a half-peel and some basic points
+    # ls = [0, 0, 0, 0.5, 0.5, 0.5, 1, 1, 1]
+    # ds = [0, 0.5, 1] * 3
+    # labels = ["C", "C3", "L", "C1", "C2", "L1", "A", "K1", "K"]
+    # ppl.plot_points_by_peel_coordinates_on_one_half_peel(ls, ds, labels)
 
     # CHOSEN_METHOD = "ebs1"
     # CHOSEN_METHOD = "upg1"
-    # CHOSEN_METHOD = "cpg1"
-    CHOSEN_METHOD = "rta1"
+    CHOSEN_METHOD = "cpg1"
+    # CHOSEN_METHOD = "rta1"
     func_pc_to_xyz = METHOD_NAME_TO_FUNCTION_POINT_CODE_TO_XYZ[CHOSEN_METHOD]
+
+    # epp.plot_point_distribution_on_one_face(func_pc_to_xyz, iterations=6)
 
-    # observe how the points lie along lines or curves on the face
-    # pcs = gpc.get_all_point_codes_from_ancestor_at_iteration(ancestor_pc="C", iterations=6)
-    pcs = gpc.get_all_point_codes_on_face_at_iteration(face_name="CAKX", iterations=4, with_edges=True, with_trailing_zeros=True)
-    ppl.plot_point_codes_on_half_peel_face_planes(pcs, face_name="CAKX", func_pc_to_xyz=func_pc_to_xyz, with_labels=False)
-
-    pcs2 = reduce((lambda x,y: x+y), [gpc.get_all_point_codes_from_ancestor_at_iteration(ancestor_pc=apc, iterations=6) for apc in "CEGIK"]) + ["A"]
-    ppl.plot_point_codes_on_sphere_3d(pcs2, func_pc_to_xyz=func_pc_to_xyz, with_labels=False)
-
-    # TODO write a function that maps from an idealized face plane triangle (ACK centered on the origin)
-    # - to an idealized sphere surface section (ACK where the vertices are centered on the origin and it bows upward in the z direction)
-    # - and plot some points on both of those views
-    # - there should only be one such mapping between the sphere surface and the face plane: the one defined by projection to/from the sphere center
-
-    # TODO write a function that maps from the triangle to itself subject to the constraints that:
-    # - it is one-to-one
-    # - it is symmetric for any of the symmetries of the triangle (D3 group): 120 deg turns about the centroid, reflections
-    # - things that map to themselves: the vertices, the centroid, the midpoints of the three edges
-    # - the edges also map to themselves (but not necessarily pointwise, and in fact they shouldn't: there needs to be distortion between vertex and midpoint)
-    # - there can be many such functions, play with them and see how the points look on the sphere surface when plotted there
-
-    # # plot the points used in the line segment tests, so I can look at where they actually are after distortion induced by projection onto the face plane
-    # pcs = sorted(set(reduce(lambda x,y: x+y, [reference_pair_l, reference_pair_dl, reference_pair_d] + test_pairs_l + test_pairs_dl + test_pairs_d)))
-    # plot_point_codes_on_half_peel_face_planes(pcs, face_name="CAKX", func_pc_to_xyz=func_pc_to_xyz)
-    # plot_point_codes_on_sphere_3d(pcs, with_labels=True)
-
+    # epp.plot_point_distribution_on_northern_half_peels(func_pc_to_xyz, iterations=6)
 
     # try a "line" of points only going in one icosa direction
     pc_init, pc_final, direction = "C010100000", "K010100000", 2
@@ -113,8 +51,6 @@ def vectors_between_pairs_of_point_codes_are_parallel(pair0, pair1, func_pc_to_x
     # pc_init, pc_final, direction = "C02020202", "K10101011", 1  # alternating staircase fractal
     # pc_init, pc_final, direction = "C0101010101", "K0101010101", 2  # alternating staircase fractal
     # pc_init, pc_final, direction = "C0100011001", "K0100011001", 2  # mostly straight with some abrupt jumps and looking a little similar to the staircase
-    pcs = get_point_path(pc_init, pc_final, direction)[:-1]
-    ppl.plot_point_codes_on_half_peel_face_planes(pcs, face_name="CAKX", func_pc_to_xyz=func_pc_to_xyz, with_labels=False)
-    ls, ds = pe.get_adjusted_peel_coordinates_of_point_codes_on_face(pcs, face_name="CAKX", func_pc_to_xyz=func_pc_to_xyz)
-    distances, angles = get_stepwise_path_distances_and_angles_2d(ls, ds)
-    plot_distances_and_angles_2d(distances, angles)
+    # epp.plot_trajectory_between_two_points_on_face_plane(pc_init, pc_final, direction, func_pc_to_xyz)
+
+    epp.report_neighbor_angle_and_distance_statistics(func_pc_to_xyz, iterations=6)

src/icosalattice/Adjacency.py

@@ -1,4 +1,5 @@
 import icosalattice.PointCodeArithmetic as pca
+import icosalattice.StartingPoints as sp
 
 
 def get_adjacency_from_point_code(pc, fix_edge_polarity=True):
@@ -12,3 +13,30 @@ def get_adjacency_from_point_code(pc, fix_edge_polarity=True):
     # print(f"{pc=} has {adj=}")
     return adj
 
+
+def get_neighbors_of_point_code(pc):
+    # go in counterclockwise direction, which is how the [1, 2, 3, -1, -2, -3] vectors go for a point with 6 neighbors
+    # around the north pole, this is the order CEGIK
+    # around the south pole, this is the order LJHFD
+    if pc[0] in sp.POLES:
+        pca.validate_point_code(pc)
+        iterations = len(pc) - 1
+        if pc[0] == sp.NORTH_POLE:
+            return [x + "1"*iterations for x in sp.NORTHERN_RING]
+        elif pc[0] == sp.SOUTH_POLE:
+            return [x + "3"*iterations for x in sp.SOUTHERN_RING[::-1]]
+        else:
+            raise Exception(f"invalid pole; shouldn't happen")
+    else:
+        adj = get_adjacency_from_point_code(pc)
+        res = []
+        for direction in [1, 2, 3, -1, -2, -3]:
+            val = adj[direction]
+            if val is not None:
+                res.append(val)
+        if all(x == "0" for x in pc[1:]):
+            n_neighbors_expected = 5
+        else:
+            n_neighbors_expected = 6
+        assert len(res) == n_neighbors_expected, f"{pc = } has wrong number of neighbors: {adj}"
+        return res

src/icosalattice/AnglesOnSphere.py

@@ -0,0 +1,38 @@
+# functions relating to measuring angles on the sphere surface
+
+import math
+import numpy as np
+
+from icosalattice.DistancesOnSphere import distance_great_circle
+
+
+# TODO will probably be good to implement generic methods for solving spherical triangles
+# TODO implement measurement of area occupied by a set of points that approximate a 2D region, e.g. a country; use sum of a bunch of triangles for this?
+
+# https://en.wikipedia.org/wiki/Spherical_trigonometry
+
+
+def measure_angles_on_sphere(xyz_A, xyz_B, xyz_C):
+    # angle ABC with vertex at B; what angle do the paths AB and BC make when they meet at B?
+    # assumes radius = 1
+    a = distance_great_circle(xyz_B, xyz_C)
+    b = distance_great_circle(xyz_A, xyz_C)
+    c = distance_great_circle(xyz_A, xyz_B)
+    cos_a = np.cos(a)
+    cos_b = np.cos(b)
+    cos_c = np.cos(c)
+    sin_a = np.sin(a)
+    sin_b = np.sin(b)
+    sin_c = np.sin(c)
+    cos_A = (cos_a - cos_b*cos_c)/(sin_b*sin_c)
+    cos_B = (cos_b - cos_c*cos_a)/(sin_c*sin_a)
+    cos_C = (cos_c - cos_a*cos_b)/(sin_a*sin_b)
+    return [np.arccos(cos_A), np.arccos(cos_B), np.arccos(cos_C)]
+
+
+def get_area_of_triangle_on_sphere(xyz1, xyz2, xyz3):
+    # Girard's theorem
+    # assumes radius = 1
+    angles = measure_angles_on_sphere(xyz1, xyz2, xyz3)
+    return sum(angles) - np.pi
+

src/icosalattice/CoordinatesByPlaneGridding.py

@@ -9,3 +9,9 @@ def get_xyz_from_point_code_using_corrected_plane_gridding(pc, as_array=True):
     l_modified, d_modified = adjust_ld_using_lp_transformation_in_triangle_coordinates(l_raw, d_raw)
     sld = spc, l_modified, d_modified
     return pe.get_xyz_from_adjusted_peel_coordinates(sld, as_array=as_array)
+
+
+def get_xyz_from_point_code_using_uncorrected_plane_gridding(pc, as_array=True):
+    # don't move the l and d coordinates, take them as-is from the point code
+    sld = pe.get_raw_peel_coordinates_from_point_code(pc)
+    return pe.get_xyz_from_adjusted_peel_coordinates(sld, as_array=as_array)

src/icosalattice/CoordinatesOfPointCode.py

@@ -25,7 +25,7 @@
 import matplotlib.pyplot as plt
 
 from icosalattice.CoordinatesByAncestry import get_xyz_from_point_code_using_ancestry
-from icosalattice.CoordinatesByPlaneGridding import get_xyz_from_point_code_using_corrected_plane_gridding
+from icosalattice.CoordinatesByPlaneGridding import get_xyz_from_point_code_using_corrected_plane_gridding, get_xyz_from_point_code_using_uncorrected_plane_gridding
 from icosalattice.CoordinatesByRThetaAdjustment import get_xyz_from_point_code_using_r_theta_adjustment
 from icosalattice.GeneratePointCodes import get_all_point_codes_from_ancestor_at_iteration, get_all_point_codes_at_iteration
 import icosalattice.PeelCoordinates as pe
@@ -40,11 +40,11 @@
 
 METHOD_NAME_TO_FUNCTION_POINT_CODE_TO_XYZ = {
     "ebs1": get_xyz_from_point_code_using_ancestry,  # "edge bisection"
-    "upg1": pe.get_xyz_from_point_code_using_peel_coordinates,  # "uncorrected plane gridding"
+    "upg1": get_xyz_from_point_code_using_uncorrected_plane_gridding,  # "uncorrected plane gridding"
     "cpg1": get_xyz_from_point_code_using_corrected_plane_gridding,  # "corrected plane gridding"
     "rta1": get_xyz_from_point_code_using_r_theta_adjustment,  # "r-theta adjustment"
 }
-CHOSEN_METHOD = "rta1"
+CHOSEN_METHOD = "cpg1"
 
 
 def get_xyz_from_point_code(pc, as_array=True):

src/icosalattice/DistancesOnSphere.py

@@ -0,0 +1,33 @@
+# functions relating to measuring distance on sphere surface
+
+import math
+import numpy as np
+
+
+
+def convert_distance_3d_to_great_circle(d0, r=1):
+    theta = 2 * np.arcsin(d0 / (2*r))
+    d_gc = r * theta
+    # assert (0 <= d_gc).all(), f"bad great circle distance {d_gc} from d0={d0}, r={r}"
+    # assert (d_gc <= np.pi * r).all(), f"bad great circle distance {d_gc} from d0={d0}, r={r}"
+    # assert ((d_gc > d0) | (abs(d_gc - d0) < 1e-9)).all(), f"shortest distance should be a straight line, but got great-circle {d_gc} from Euclidean {d0}"
+    # print(f"d0 = {d0}, r = {r} -> great circle distance {d_gc}")
+    return d_gc
+    # return (np.vectorize(lambda d: UnitSpherePoint.convert_distance_3d_to_great_circle_single_value(d, radius=radius)))(d0)
+
+
+def convert_distance_great_circle_to_3d(d_gc, r=1):
+    theta = d_gc / r
+    d0 = 2 * r * np.sin(theta)
+    assert 0 <= d0 <= 2*r, f"bad 3d distance {d0} from d_gc={d_gc}, r={r}"
+    assert d0 <= d_gc, "shortest distance should be a straight line"
+    return d0
+
+
+def distance_3d(xyz1, xyz2):
+    return np.linalg.norm(xyz1 - xyz2)
+
+
+def distance_great_circle(xyz1, xyz2, r=1):
+    d0 = distance_3d(xyz1, xyz2)
+    return convert_distance_3d_to_great_circle(d0, r)

src/icosalattice/EvaluatePointPlacementMethods.py

@@ -0,0 +1,106 @@
+from functools import reduce
+import numpy as np
+import matplotlib.pyplot as plt
+
+import icosalattice.GeneratePointCodes as gpc
+import icosalattice.PlotPointLocations as ppl
+from icosalattice.PointPaths import get_point_path, get_stepwise_path_distances_and_angles_2d
+from icosalattice.PlotPaths import plot_distances_and_angles_2d
+import icosalattice.PeelCoordinates as pe
+from icosalattice.Adjacency import get_neighbors_of_point_code
+from icosalattice.AnglesOnSphere import measure_angles_on_sphere
+from icosalattice.DistancesOnSphere import distance_great_circle
+
+
+
+def plot_point_distribution_on_one_face(func_pc_to_xyz, iterations=4):
+    # observe how the points lie along lines or curves on the face
+    # pcs = gpc.get_all_point_codes_from_ancestor_at_iteration(ancestor_pc="C", iterations=6)
+    pcs = gpc.get_all_point_codes_on_face_at_iteration(face_name="CAKX", iterations=iterations, with_edges=True, with_trailing_zeros=True)
+    ppl.plot_point_codes_on_half_peel_face_planes(pcs, face_name="CAKX", func_pc_to_xyz=func_pc_to_xyz, with_labels=False)
+    # for future self searching for code that made plots:
+    # CAKX_i4_{method}.png
+
+
+def plot_point_distribution_on_northern_half_peels(func_pc_to_xyz, iterations=6):
+    pcs = reduce((lambda x,y: x+y), [gpc.get_all_point_codes_from_ancestor_at_iteration(ancestor_pc=apc, iterations=iterations) for apc in "CEGIK"]) + ["A"]
+    ppl.plot_point_codes_on_sphere_3d(pcs, func_pc_to_xyz=func_pc_to_xyz, with_labels=False)
+
+
+def plot_trajectory_between_two_points_on_face_plane(pc_init, pc_final, direction, func_pc_to_xyz):
+    pcs = get_point_path(pc_init, pc_final, direction)[:-1]
+    ppl.plot_point_codes_on_half_peel_face_planes(pcs, face_name="CAKX", func_pc_to_xyz=func_pc_to_xyz, with_labels=False)
+    ls, ds = pe.get_adjusted_peel_coordinates_of_point_codes_on_face(pcs, face_name="CAKX", func_pc_to_xyz=func_pc_to_xyz)
+    distances, angles = get_stepwise_path_distances_and_angles_2d(ls, ds)
+    plot_distances_and_angles_2d(distances, angles)
+    # for future self searching for code that made plots:
+    # C010100000_to_K010100000_steps_{method}.png
+    # {pc_init}_to_{pc_final}_steps_{method}.png
+
+
+def report_neighbor_angle_and_distance_statistics(func_pc_to_xyz, iterations=6):
+    pcs = gpc.get_all_point_codes_at_iteration(iterations=iterations, with_trailing_zeros=True)
+    pc_to_xyz = {pc: func_pc_to_xyz(pc) for pc in pcs}
+    angles = []
+
+    # debug: finding aberrantly large angles
+    angles_larger_than_72 = {}
+
+    distances = []
+    for i, pc in enumerate(pcs):
+        if i % 1000 == 0:
+            print(f"getting distances and angles: {i = } / {len(pcs)}")
+        xyz0 = pc_to_xyz[pc]
+        neighbors = get_neighbors_of_point_code(pc)
+        # print(pc, neighbors)
+        these_angles_deg = []
+        these_distances = []
+
+        # there is double counting and redundant computation going on here
+        # (since the distances are used in computing the angles, we compute the same angles more than once, etc.) 
+        # but I don't care right now unless it takes very long to finish
+
+        for i in range(len(neighbors)):
+            n1 = neighbors[i]
+            n2 = neighbors[(i+1) % len(neighbors)]
+            xyz1 = pc_to_xyz[n1]
+            xyz2 = pc_to_xyz[n2]
+            angle1, angle0, angle2 = measure_angles_on_sphere(xyz1, xyz0, xyz2)
+            angle0_deg = angle0 * 180/np.pi
+
+            # debug: finding aberrantly large angles
+            if angle0_deg > 72 + 1e-9:
+                a = round(angle0_deg, 6)
+                if a not in angles_larger_than_72:
+                    angles_larger_than_72[a] = []
+                angles_larger_than_72[a].append((n1, pc, n2))
+
+            distance = distance_great_circle(xyz0, xyz1)
+            these_angles_deg.append(angle0_deg)
+            these_distances.append(distance)
+        angles += these_angles_deg
+        distances += these_distances
+
+    print("\nlarge angles:\n" + ("\n".join(f"  {a:f} deg: {tup}" for a, tups in sorted(angles_larger_than_72.items(), reverse=True) for tup in tups) if len(angles_larger_than_72) > 0 else "(none)") + "\n")
+
+    print(f"{min(angles) = :f}")
+    print(f"mean(angles) = {np.mean(angles):f}")
+    print(f"median(angles) = {np.median(angles):f}")
+    print(f"{max(angles) = :f}\n")
+
+    print(f"{min(distances) = :f}")
+    print(f"mean(distances) = {np.mean(distances):f}")
+    print(f"median(distances) = {np.median(distances):f}")
+    print(f"{max(distances) = :f}\n")
+
+    plt.hist(angles, bins=100)
+    plt.title("angles")
+    plt.show()
+    # for future self searching for code that made plots:
+    # neighbor_angle_distribution_{method}.png
+
+    plt.hist(distances, bins=100)
+    plt.title("distances")
+    plt.show()
+    # for future self searching for code that made plots:
+    # neighbor_distance_distribution_{method}.png

src/icosalattice/PeelCoordinates.py

@@ -296,10 +296,11 @@ def get_adjusted_peel_coordinates_from_xyz(xyz):
 
 
 def get_xyz_from_point_code_using_peel_coordinates(pc, as_array=True):
-    sld_raw = get_raw_peel_coordinates_from_point_code(pc)
-    sld_adjusted = get_adjusted_peel_coordinates_from_raw_peel_coordinates(sld_raw)
-    xyz = get_xyz_from_adjusted_peel_coordinates(sld_adjusted, as_array=as_array)
-    return xyz
+    raise DeprecationWarning("use functions in icosalattice.CoordinatesByPlaneGridding instead")
+    # sld_raw = get_raw_peel_coordinates_from_point_code(pc)
+    # sld_adjusted = get_adjusted_peel_coordinates_from_raw_peel_coordinates(sld_raw)
+    # xyz = get_xyz_from_adjusted_peel_coordinates(sld_adjusted, as_array=as_array)
+    # return xyz
 
 
 def get_point_code_from_xyz_using_peel_coordinates(xyz):