first round of review comments

Author: Gregory Roberts

tests/test_components/test_autograd_rf_box.py

@@ -25,7 +25,7 @@
 SAVE_ADJ_LOC = 1
 LOCAL_GRADIENT = True
 VERBOSE = False
-SHOW_PRINT_STATEMENTS = False
+SHOW_PRINT_STATEMENTS = True
 USE_POLYSLAB_FOR_BOX = False
 NUMERICAL_RESULTS_DATA_DIR = (
     "./numerical_rf_box_polyslab_test/" if USE_POLYSLAB_FOR_BOX else "./numerical_rf_box_box_test/"
@@ -366,7 +366,7 @@ def run_and_process_fd(all_box_parameters, fd_step, objective):
         fd_grad = np.squeeze(all_fd_grad_parameters[argmin_convergence_test])
     else:
         fd_grad = np.squeeze(all_fd_grad_parameters)
-        valid_mask = np.ones(fd_grad, dtype=bool)
+        valid_mask = np.ones(fd_grad.shape, dtype=bool)
 
     return fd_grad, valid_mask
 
@@ -473,7 +473,6 @@ def test_finite_difference_2d_box_pec(rf_2d_test_parameters, rng, tmp_path, crea
         "test_number",
     )(rf_2d_test_parameters)
 
-    dim_um = 3 * mesh_wvl_um
     dim_um = 3 * mesh_wvl_um
 
     thickness_box_placement_um = 0.6 * mesh_wvl_um  # 1.2 * mesh_wvl_um
@@ -497,8 +496,6 @@ def test_finite_difference_2d_box_pec(rf_2d_test_parameters, rng, tmp_path, crea
     else:
         fd_step = np.array([mesh_cell_override_size])
 
-    all_perm = []
-
     z_planes = np.array([0.0])
 
     all_box_parameters, box_z_values = generate_boxes(
@@ -632,7 +629,7 @@ def test_finite_difference_2d_box_pec(rf_2d_test_parameters, rng, tmp_path, crea
     ),
     indirect=["dir_name"],
 )
-def test_finite_difference_3d_pec(rf_3d_test_parameters, rng, tmp_path, create_directory):
+def test_finite_difference_3d_box_pec(rf_3d_test_parameters, rng, tmp_path, create_directory):
     """Test a variety of autograd permittivity gradients for 3D PEC boxes by"""
     """comparing them to numerical finite difference."""
 
@@ -656,7 +653,6 @@ def test_finite_difference_3d_pec(rf_3d_test_parameters, rng, tmp_path, create_d
         "test_number",
     )(rf_3d_test_parameters)
 
-    dim_um = 3 * mesh_wvl_um
     dim_um = 3 * mesh_wvl_um
 
     thickness_box_placement_um = 1.2 * mesh_wvl_um
@@ -680,8 +676,6 @@ def test_finite_difference_3d_pec(rf_3d_test_parameters, rng, tmp_path, create_d
     else:
         fd_step = np.array([mesh_cell_override_size])
 
-    all_perm = []
-
     all_box_parameters = generate_boxes(
         lateral_dim_bounds=[-0.5 * dim_um, 0.5 * dim_um],
         mesh_cell_override_size=mesh_cell_override_size,

tests/test_components/test_autograd_rf_polyslab.py

@@ -327,7 +327,7 @@ def run_and_process_fd(polyslab_parameters, fd_step, objective):
         fd_grad = np.squeeze(all_fd_grad_parameters[argmin_convergence_test])
     else:
         fd_grad = np.squeeze(all_fd_grad_parameters)
-        valid_mask = np.ones(fd_grad, dtype=bool)
+        valid_mask = np.ones(fd_grad.shape, dtype=bool)
 
     return fd_grad, valid_mask
 
@@ -398,6 +398,10 @@ def run_and_process_fd(polyslab_parameters, fd_step, objective):
 
             test_number += 1
 
+import sys
+
+sys.stdout = sys.stderr
+
 
 @pytest.mark.numerical
 @pytest.mark.parametrize(
@@ -408,7 +412,7 @@ def run_and_process_fd(polyslab_parameters, fd_step, objective):
     ),
     indirect=["dir_name"],
 )
-def test_finite_difference_2d_pec(rf_2d_test_parameters, rng, tmp_path, create_directory):
+def test_finite_difference_2d_polyslab_pec(rf_2d_test_parameters, rng, tmp_path, create_directory):
     """Test a variety of autograd permittivity gradients for 2D `PolySlab` PEC by"""
     """comparing them to numerical finite difference."""
 
@@ -432,7 +436,6 @@ def test_finite_difference_2d_pec(rf_2d_test_parameters, rng, tmp_path, create_d
         "test_number",
     )(rf_2d_test_parameters)
 
-    dim_um = 1.5 * mesh_wvl_um
     dim_um = 1.5 * mesh_wvl_um
 
     thickness_box_placement_um = 1.2 * mesh_wvl_um
@@ -456,8 +459,6 @@ def test_finite_difference_2d_pec(rf_2d_test_parameters, rng, tmp_path, create_d
     else:
         fd_step = np.array([mesh_cell_override_size])
 
-    all_perm = []
-
     z_planes = np.linspace(-0.5 * thickness_box_placement_um, 0.5 * thickness_box_placement_um, 3)
 
     polyslab = generate_polyslab(0.5 * dim_um, dim_um, rng)
@@ -485,16 +486,16 @@ def test_finite_difference_2d_pec(rf_2d_test_parameters, rng, tmp_path, create_d
 
     obj, adj_grad = obj_val_and_grad([polyslab])
 
-    if test_number > 3:
-        old_data = np.load(
-            f"{NUMERICAL_RESULTS_DATA_DIR}/rf_polyslab_2d_vertex_fd_adj_data_test_{test_number}.npy"
-        )
-        fd_grad = old_data[SAVE_FD_LOC, :]
-        valid_mask = np.ones(len(fd_grad)) > 0
-    else:
-        fd_grad, valid_mask = run_and_process_fd(
-            polyslab_parameters=polyslab, fd_step=fd_step, objective=objective
-        )
+    # if test_number > 3:
+    #     old_data = np.load(
+    #         f"{NUMERICAL_RESULTS_DATA_DIR}/rf_polyslab_2d_vertex_fd_adj_data_test_{test_number}.npy"
+    #     )
+    #     fd_grad = old_data[SAVE_FD_LOC, :]
+    #     valid_mask = np.ones(len(fd_grad)) > 0
+    # else:
+    fd_grad, valid_mask = run_and_process_fd(
+        polyslab_parameters=polyslab, fd_step=fd_step, objective=objective
+    )
 
     fd_grad = np.array(fd_grad)
     adj_grad = np.squeeze(np.array(adj_grad))
@@ -575,7 +576,7 @@ def test_finite_difference_2d_pec(rf_2d_test_parameters, rng, tmp_path, create_d
     ),
     indirect=["dir_name"],
 )
-def test_finite_difference_3d_pec(rf_3d_test_parameters, rng, tmp_path, create_directory):
+def test_finite_difference_3d_polyslab_pec(rf_3d_test_parameters, rng, tmp_path, create_directory):
     """Test a variety of autograd permittivity gradients for 3D PEC `PolySlab` by"""
     """comparing them to numerical finite difference."""
 
@@ -599,7 +600,6 @@ def test_finite_difference_3d_pec(rf_3d_test_parameters, rng, tmp_path, create_d
         "test_number",
     )(rf_3d_test_parameters)
 
-    dim_um = 1.5 * mesh_wvl_um
     dim_um = 1.5 * mesh_wvl_um
 
     thickness_box_placement_um = 1.0 * mesh_wvl_um
@@ -623,8 +623,6 @@ def test_finite_difference_3d_pec(rf_3d_test_parameters, rng, tmp_path, create_d
     else:
         fd_step = np.array([mesh_cell_override_size])
 
-    all_perm = []
-
     polyslab = generate_polyslab(0.5 * dim_um, dim_um, rng)
     polyslab_z_value = 0
     polyslab_z_thickness = polyslab_z_thickness_wvl * adj_wvl_um
@@ -651,11 +649,9 @@ def test_finite_difference_3d_pec(rf_3d_test_parameters, rng, tmp_path, create_d
 
     obj, adj_grad = obj_val_and_grad([polyslab])
 
-    old_data = np.load(
-        f"{NUMERICAL_RESULTS_DATA_DIR}/rf_polyslab_3d_vertex_fd_adj_data_test_{test_number}.npy"
+    fd_grad, valid_mask = run_and_process_fd(
+        polyslab_parameters=polyslab, fd_step=fd_step, objective=objective
     )
-    fd_grad = old_data[SAVE_FD_LOC, :]
-    valid_mask = np.ones(len(fd_grad)) > 0
 
     fd_grad = np.array(fd_grad)
     adj_grad = np.squeeze(np.array(adj_grad))
@@ -693,7 +689,7 @@ def test_finite_difference_3d_pec(rf_3d_test_parameters, rng, tmp_path, create_d
 
     if SAVE_FD_ADJ_DATA:
         np.save(
-            f"{NUMERICAL_RESULTS_DATA_DIR}/rf_polyslab_3d_vertex_fd_adj_data_test_{test_number}_.npy",
+            f"{NUMERICAL_RESULTS_DATA_DIR}/rf_polyslab_3d_vertex_fd_adj_data_test_{test_number}.npy",
             vertex_data,
         )
 
@@ -723,5 +719,5 @@ def test_finite_difference_3d_pec(rf_3d_test_parameters, rng, tmp_path, create_d
             )
 
         plt.savefig(
-            f"{NUMERICAL_RESULTS_DATA_DIR}/rf_3d_polyslab_summary_plot_test_{test_number}_.png"
+            f"{NUMERICAL_RESULTS_DATA_DIR}/rf_3d_polyslab_summary_plot_test_{test_number}.png"
         )

tests/test_components/test_geometry.py

@@ -1298,3 +1298,136 @@ def test_cleanup_shapely_object():
     orig_polygon = shapely.Polygon(exterior_coords)
     new_polygon = cleanup_shapely_object(orig_polygon, tolerance_ratio=1e-12)
     assert len(new_polygon.exterior.coords) == 0  # empty / collinear polygons should get deleted
+
+
+def test_snap_coords_outside():
+    """Test coordinate snapping for box geometries."""
+    num_coords = 101
+    coords_to_snap = np.linspace(-1, 1, num_coords)
+
+    coord_face_min = np.mean(coords_to_snap[0:2])
+    coord_face_max = np.mean(coords_to_snap[-2:])
+
+    snap_pt_min = td.Box._snap_coords_outside(
+        min_max_index=0, snap_coords_values=coords_to_snap, coord_normal_face=coord_face_min
+    )
+    snap_pt_max = td.Box._snap_coords_outside(
+        min_max_index=1, snap_coords_values=coords_to_snap, coord_normal_face=coord_face_max
+    )
+
+    assert snap_pt_min == coords_to_snap[0], "Unexpected snapping point minimum"
+    assert snap_pt_max == coords_to_snap[-1], "Unexpected snapping point maximum"
+
+    with AssertLogLevel("WARNING", contains_str="Unable to snap coordinates outside"):
+        snap_pt_min = td.Box._snap_coords_outside(
+            min_max_index=0,
+            snap_coords_values=coords_to_snap,
+            coord_normal_face=coords_to_snap[0] - 0.1,
+        )
+
+        assert snap_pt_min == coords_to_snap[0], "Unexpected snapping point minimum"
+
+    with AssertLogLevel("WARNING", contains_str="Unable to snap coordinates outside"):
+        snap_pt_min = td.Box._snap_coords_outside(
+            min_max_index=1,
+            snap_coords_values=coords_to_snap,
+            coord_normal_face=coords_to_snap[-1] + 0.1,
+        )
+
+        assert snap_pt_min == coords_to_snap[-1], "Unexpected snapping point maxium"
+
+
+def test_singularity_correction_pec():
+    """Tests singularity correction detection used for PEC gradients of `Box` geometries."""
+
+    check_3d_size = (1.0, 2.0, 3.0)
+
+    for axis in range(3):
+        is_2d, zero_dimension, do_singularity_correction = td.Box._check_singularity_correction_pec(
+            size=check_3d_size, axis_normal=axis
+        )
+        assert not is_2d, "Unexpectedly found 2D geometry"
+        assert not zero_dimension, "Unexpectedly found zero_dimension"
+        assert not do_singularity_correction, (
+            "Unexpectedly identifying need for singularity correction"
+        )
+
+    check_1d_error_size = (1.0, 0.0, 0.0)
+
+    for axis in range(3):
+        with AssertLogLevel(
+            "ERROR", contains_str="Derivative of PEC material with less than 2 dimensions"
+        ):
+            is_2d, zero_dimension, do_singularity_correction = (
+                td.Box._check_singularity_correction_pec(size=check_1d_error_size, axis_normal=axis)
+            )
+
+    check_2d_size = 1.0
+
+    for index_2d in range(3):
+        check_2d_size = tuple(0.0 if idx == index_2d else 1.0 for idx in range(3))
+
+        for axis in range(3):
+            print(f"check 2d size = {check_2d_size}")
+
+            is_2d, zero_dimension, do_singularity_correction = (
+                td.Box._check_singularity_correction_pec(size=check_2d_size, axis_normal=axis)
+            )
+
+            print(f"is 2d = {is_2d}")
+
+            assert is_2d, "Expected 2D geometry detection"
+            assert zero_dimension == "xyz"[index_2d], "Wrong zero dimension identified"
+
+            if axis == index_2d:
+                assert not do_singularity_correction, (
+                    "Unexpectedly detecting need for singularity correction"
+                )
+            else:
+                assert do_singularity_correction, (
+                    "Unexpectedly not detecting need for singularity correction"
+                )
+
+
+def test_trim_dims_and_bounds_edge():
+    """Tests the dimension and bound trimming for 2D PEC gradient integrations in `Box`."""
+
+    for zero_dimension_idx in range(3):
+        box = td.Box(
+            center=(0.0, 0.0, 0.0),
+            size=tuple(0.0 if zero_dimension_idx == idx else 1.0 * (idx + 1) for idx in range(3)),
+        )
+
+        for axis_normal in range(3):
+            if axis_normal == zero_dimension_idx:
+                continue
+
+            dim_normal, dims_perp = box.pop_axis("xyz", axis=axis_normal)
+            bounds_normal, bounds_perp = box.pop_axis(np.array(box.bounds).T, axis=axis_normal)
+            bounds_perp = np.array(bounds_perp).T
+
+            trimmed_dims, trimmed_bounds = td.Box._trim_dims_and_bounds_edge(
+                dims_perp=dims_perp,
+                bounds_perp=bounds_perp,
+                zero_dimension="xyz"[zero_dimension_idx],
+            )
+
+            expected_trimmed_idx = 0
+            for idx in range(3):
+                if (idx == axis_normal) or (idx == zero_dimension_idx):
+                    continue
+
+                expected_trimmed_idx = idx
+
+            expected_trimmed_dims = ("xyz"[expected_trimmed_idx],)
+            expected_trimmed_bounds = (
+                (-0.5 * (expected_trimmed_idx + 1),),
+                (0.5 * (expected_trimmed_idx + 1),),
+            )
+
+            assert np.all(np.array(expected_trimmed_dims) == np.array(trimmed_dims)), (
+                "Unexpected trimmed dims"
+            )
+            assert np.all(np.array(expected_trimmed_bounds) == np.array(trimmed_bounds)), (
+                "Unexpected trimmed bounds"
+            )