Add test for microstructure macro sensitivities

splinepy/microstructure/microstructure.py

@@ -427,8 +427,8 @@ def _compute_tiling_prerequisites(
     def create(
         self,
         closing_face=None,
-        knot_span_wise=None,
-        macro_sensitivities=None,
+        knot_span_wise=True,
+        macro_sensitivities=False,
         **kwargs,
     ):
         """Create a Microstructure.
@@ -455,11 +455,12 @@ def create(
         if not self._sanity_check():
             raise ValueError("Not enough information provided, abort")
 
-        # Set default values
-        if knot_span_wise is None:
-            knot_span_wise = True
-        if macro_sensitivities is None:
-            macro_sensitivities = False
+        assert isinstance(
+            knot_span_wise, bool
+        ), "knot_span_wise must be a bool"
+        assert isinstance(
+            macro_sensitivities, bool
+        ), "macro_senstivities must be a bool"
 
         # check if user wants closed structure
         if closing_face is not None:
@@ -533,7 +534,7 @@ def create(
         if macro_sensitivities or parameter_sensitivities:
             spline_list_derivs = [
                 []
-                for i in range(
+                for _ in range(
                     n_parameter_sensitivities + n_macro_sensitivities
                 )
             ]

tests/test_microstructure.py

@@ -1,3 +1,5 @@
+import numpy as np
+
 import splinepy.microstructure as ms
 from splinepy.helpme.create import box
 from splinepy.helpme.integrate import volume
@@ -13,6 +15,8 @@
 
 
 def test_closing_face():
+    """Check if closing face is working"""
+
     def check_if_closed(multipatch, closure_direction):
         """Helper function to see if multipatch has a closing surface
 
@@ -65,7 +69,7 @@ def max_identifier(points):
         generator = ms.microstructure.Microstructure(
             deformation_function=box(*BOX_DIMENSIONS[: tile_creator._dim]),
             microtile=tile_creator,
-            tiling=TILING[: tile_creator._para_dim],
+            tiling=TILING[: tile_creator._dim],
         )
         # Go through all implemented closure direction
         closure_directions = {
@@ -77,5 +81,71 @@ def max_identifier(points):
             check_if_closed(multipatch, closure_direction)
 
 
-def test_macro_sensitivities():
-    pass
+def test_macro_sensitivities(np_rng, heps=1e-7, n_test_points=10):
+    """Testing the correctness of the derivatives of the whole microstructure w.r.t.
+    the deformation function's control points. It is tested by evaluating the derivative
+    obtained via finite differences. The values are evaluated at random points.
+
+    Parameters
+    ----------
+    heps: float
+        Perturbation size for finite difference evaluation
+    n_test_points: int
+        Number of testing points int the parametric domain"""
+
+    for tile_class in all_tile_classes:
+        tile_creator = tile_class()
+        deformation_function_orig = box(*BOX_DIMENSIONS[: tile_creator._dim])
+        generator = ms.microstructure.Microstructure(
+            deformation_function=deformation_function_orig,
+            microtile=tile_creator,
+            tiling=TILING[: tile_creator._dim],
+        )
+        multipatch = generator.create(macro_sensitivities=True)
+        dim = multipatch.dim
+        n_cps = deformation_function_orig.cps.shape[0]
+
+        # Set evaluation points as random spots in the parametric space
+        eval_points = np_rng.random((n_test_points, tile_creator._para_dim))
+        microstructure_orig_evaluations = [
+            patch.evaluate(eval_points) for patch in multipatch.patches
+        ]
+        n_patches = len(multipatch.patches)
+
+        # Go through derivatives of every deformation function's control point
+        for ii_ctps in range(n_cps):
+            # Gradient through finite differences
+            deformation_function_perturbed = deformation_function_orig.copy()
+            deformation_function_perturbed.cps[ii_ctps, :] += heps
+            generator.deformation_function = deformation_function_perturbed
+            multipatch_perturbed = generator.create()
+            microstructure_perturbed_evaluations = [
+                patch.evaluate(eval_points)
+                for patch in multipatch_perturbed.patches
+            ]
+            # Evaluate finite difference gradient
+            fd_sensitivity = [
+                (patch_perturbed - patch_orig) / heps
+                for patch_perturbed, patch_orig in zip(
+                    microstructure_orig_evaluations,
+                    microstructure_perturbed_evaluations,
+                )
+            ]
+
+            # Go through each direction
+            for jj_dim in range(dim):
+                deriv_orig = multipatch.fields[ii_ctps * dim + jj_dim]
+                deriv_evaluations = [
+                    patch.evaluate(eval_points) for patch in deriv_orig.patches
+                ]
+                for k_patch, patch_deriv_implemented, patch_deriv_fd in zip(
+                    range(n_patches), deriv_evaluations, fd_sensitivity
+                ):
+                    assert np.allclose(
+                        -patch_deriv_implemented[:, jj_dim],
+                        patch_deriv_fd[:, jj_dim],
+                    ), (
+                        "Implemented derivative calculation for tile class"
+                        + f"{tile_class}, at patch {k_patch+1}/{n_patches} does not "
+                        + "match the derivative obtained using Finite Differences"
+                    )

tests/test_microtiles.py

@@ -76,7 +76,7 @@ def test_tile_class():
                 required_param in create_parameters
             ), f"create_tile() must have '{required_param}' as an input parameter"
 
-        # Ensure closure can be correctly handled
+        # Ensure closure can be handled correctly
         if "closure" in create_parameters:
             assert "_closure_directions" in members, (
                 "Tile class has closure ability. The available closure directions "
@@ -130,7 +130,7 @@ def test_tile_bounds():
 
 
 def test_tile_closure():
-    """Check if closing tiles also lie in unit cube. TODO: also check derivatives."""
+    """Check if closing tiles also lie in unit cube."""
 
     for tile_class in all_tile_classes:
         # Skip tile if if does not support closure
@@ -218,7 +218,7 @@ def test_tile_derivatives(np_rng, heps=1e-7, n_test_points=10):
             splines_orig, _ = tile_creator.create_tile(
                 parameters=parameters, closure=closure
             )
-            # Set evaluation points as 4 random spots in the parametric space
+            # Set evaluation points as random spots in the parametric space
             eval_points = np_rng.random(
                 (n_test_points, splines_orig[0].para_dim)
             )