Fix off-by-one filtering errors.

python/Makefile.am

@@ -38,6 +38,7 @@ TESTS =                                        \
     $(TEST_DIR)/test_3rd_harm_1d.py            \
     $(TEST_DIR)/test_absorber_1d.py            \
     $(TEST_DIR)/test_adjoint_solver.py         \
+    $(TEST_DIR)/test_adjoint_utils.py          \
     $(TEST_DIR)/test_adjoint_cyl.py            \
     $(TEST_DIR)/test_adjoint_jax.py            \
     $(TEST_DIR)/test_antenna_radiation.py      \

python/adjoint/filters.py

@@ -7,6 +7,19 @@
 import meep as mp
 from scipy import special
 
+def compute_mg_dims(Lx,Ly,resolution):
+    '''Compute the material grid dimensions from
+    the corresponding resolution, x-size, and y-size.
+    The grid dimensions must be odd.    
+    '''
+    Nx = int(Lx * resolution)
+    Ny = int(Ly * resolution)
+    if (Nx % 2) == 0:
+        Nx += 1
+    if (Ny % 2) == 0:
+        Ny += 1
+    return Nx, Ny
+
 
 def _centered(arr, newshape):
     '''Helper function that reformats the padded array of the fft filter operation.
@@ -98,8 +111,7 @@ def simple_2d_filter(x, kernel, Lx, Ly, resolution, symmetries=[]):
         The output of the 2d convolution.
     """
     # Get 2d parameter space shape
-    Nx = int(Lx * resolution) + 1
-    Ny = int(Ly * resolution) + 1
+    Nx, Ny = compute_mg_dims(Lx,Ly,resolution)
     (kx, ky) = kernel.shape
 
     # Adjust parameter space shape for symmetries
@@ -135,8 +147,8 @@ def simple_2d_filter(x, kernel, Lx, Ly, resolution, symmetries=[]):
     # Convolution (multiplication in frequency domain)
     Y = H * X
 
-    # We need to fftshift since we padded both sides if each dimension of our input and kernel.
-    y = npa.fft.fftshift(npa.real(npa.fft.ifft2(Y)))
+    # We need to fftshift since we padded both sides of each dimension of our input and kernel.
+    y = npa.fft.ifftshift(npa.real(npa.fft.ifft2(Y)))
 
     # Remove all the extra padding
     y = _centered(y, (kx, ky))
@@ -179,8 +191,7 @@ def cylindrical_filter(x, radius, Lx, Ly, resolution, symmetries=[]):
     density-based topology optimization. Archive of Applied Mechanics, 86(1-2), 189-218.
     '''
     # Get 2d parameter space shape
-    Nx = int(Lx * resolution) + 1
-    Ny = int(Ly * resolution) + 1
+    Nx, Ny = compute_mg_dims(Lx,Ly,resolution)
 
     # Formulate grid over entire design region
     xv, yv = np.meshgrid(np.linspace(-Lx / 2, Lx / 2, Nx),
@@ -189,7 +200,7 @@ def cylindrical_filter(x, radius, Lx, Ly, resolution, symmetries=[]):
                          indexing='ij')
 
     # Calculate kernel
-    kernel = np.where(np.abs(xv**2 + yv**2) <= radius**2, 1, 0).T
+    kernel = np.where(np.abs(xv**2 + yv**2) <= radius**2, 1, 0)
 
     # Normalize kernel
     kernel = kernel / np.sum(kernel.flatten())  # Normalize the filter
@@ -229,8 +240,7 @@ def conic_filter(x, radius, Lx, Ly, resolution, symmetries=[]):
     density-based topology optimization. Archive of Applied Mechanics, 86(1-2), 189-218.
     '''
     # Get 2d parameter space shape
-    Nx = int(Lx * resolution) + 1
-    Ny = int(Ly * resolution) + 1
+    Nx, Ny = compute_mg_dims(Lx,Ly,resolution)
 
     # Formulate grid over entire design region
     xv, yv = np.meshgrid(np.linspace(-Lx / 2, Lx / 2, Nx),
@@ -279,8 +289,7 @@ def gaussian_filter(x, sigma, Lx, Ly, resolution, symmetries=[]):
     topology-optimized metasurfaces. Optical Materials Express, 9(2), 469-482.
     '''
     # Get 2d parameter space shape
-    Nx = int(Lx * resolution) + 1
-    Ny = int(Ly * resolution) + 1
+    Nx, Ny = compute_mg_dims(Lx,Ly,resolution)
 
     gaussian = lambda x, sigma: np.exp(-x**2 / sigma**2)
 

python/tests/test_adjoint_cyl.py

@@ -26,8 +26,7 @@
 design_region_resolution = int(2*resolution)
 design_r = 4.8
 design_z = 2
-Nr = int(design_region_resolution*design_r) + 1
-Nz = int(design_region_resolution*design_z) + 1
+Nr, Nz = mpa.compute_mg_dims(design_r,design_z,design_region_resolution)
 
 fcen = 1/1.55
 width = 0.2

python/tests/test_adjoint_jax.py

@@ -69,9 +69,7 @@ def build_straight_wg_simulation(
             center=[sx / 2 - pml_width - source_to_pml, 0, 0],
         ),
     ]
-
-    nx = int(design_region_resolution * design_region_shape[0]) + 1
-    ny = int(design_region_resolution * design_region_shape[1]) + 1
+    nx, ny = mpa.compute_mg_dims(design_region_shape[0],design_region_shape[1],design_region_resolution)
     mat_grid = mp.MaterialGrid(
         mp.Vector3(nx, ny),
         sio2,

python/tests/test_adjoint_solver.py

@@ -29,8 +29,7 @@
 
 design_region_size = mp.Vector3(1.5,1.5)
 design_region_resolution = int(2*resolution)
-Nx = int(design_region_resolution*design_region_size.x) + 1
-Ny = int(design_region_resolution*design_region_size.y) + 1
+Nx, Ny = mpa.compute_mg_dims(design_region_size.x,design_region_size.y,design_region_resolution)
 
 ## ensure reproducible results
 rng = np.random.RandomState(9861548)

python/tests/test_adjoint_utils.py

@@ -0,0 +1,48 @@
+'''test_adjoint_utils.py
+
+Check various components of the adjoint solver codebase, like
+filters, which may not need explicit gradient computation
+(i.e. forward and adjoint runs).
+
+'''
+
+import meep as mp
+try:
+    import meep.adjoint as mpa
+except:
+    import adjoint as mpa
+import numpy as np
+from autograd import numpy as npa
+from autograd import tensor_jacobian_product
+import unittest
+from enum import Enum
+from utils import ApproxComparisonTestCase
+import parameterized
+
+_TOL = 1e-6 if mp.is_single_precision() else 1e-14
+
+## ensure reproducible results
+rng = np.random.RandomState(9861548)
+
+class TestAdjointUtils(ApproxComparisonTestCase):
+    @parameterized.parameterized.expand([
+    ('1.0_1.0_20_conic',1.0, 1.0, 20, 0.24, mpa.conic_filter),
+    ('1.0_1.0_23_conic',1.0, 1.0, 23, 0.24, mpa.conic_filter),
+    ('0.887_1.56_conic',0.887, 1.56, 20, 0.24, mpa.conic_filter),
+    ('0.887_1.56_conic',0.887, 1.56, 31, 0.24, mpa.conic_filter),
+    ('0.887_1.56_gaussian',0.887, 1.56, 20, 0.24, mpa.gaussian_filter),
+    ('0.887_1.56_cylindrical',0.887, 1.56, 20, 0.24, mpa.cylindrical_filter)
+    ])
+    def test_filter_offset(self,test_name,Lx,Ly,resolution,radius,filter_func):
+        '''ensure that the filters are indeed zero-phase'''
+        print("Testing ",test_name)
+        Nx, Ny = mpa.compute_mg_dims(Lx,Ly,resolution)
+        x = np.random.rand(Nx,Ny)
+        x = x + np.fliplr(x)
+        x = x + np.flipud(x)
+        y = filter_func(x, radius, Lx, Ly, resolution)
+        self.assertClose(y,np.fliplr(y),epsilon=_TOL)
+        self.assertClose(y,np.flipud(y),epsilon=_TOL)
+
+if __name__ == '__main__':
+    unittest.main()

python/tests/test_get_epsilon_grid.py

@@ -2,6 +2,10 @@
 import parameterized
 import numpy as np
 import meep as mp
+try:
+    import meep.adjoint as mpa
+except:
+    import adjoint as mpa
 from meep.materials import SiN, Co
 
 class TestGetEpsilonGrid(unittest.TestCase):
@@ -12,8 +16,7 @@ def setUp(self):
 
         matgrid_resolution = 200
         matgrid_size = mp.Vector3(1.0,1.0,mp.inf)
-        Nx = int(matgrid_resolution*matgrid_size.x) + 1
-        Ny = int(matgrid_resolution*matgrid_size.y) + 1
+        Nx, Ny = mpa.compute_mg_dims(matgrid_size.x,matgrid_size.y,matgrid_resolution)
         x = np.linspace(-0.5*matgrid_size.x,0.5*matgrid_size.x,Nx)
         y = np.linspace(-0.5*matgrid_size.y,0.5*matgrid_size.y,Ny)
         xv, yv = np.meshgrid(x,y)

python/tests/test_material_grid.py

@@ -1,4 +1,8 @@
 import meep as mp
+try:
+    import meep.adjoint as mpa
+except:
+    import adjoint as mpa
 import numpy as np
 from scipy.ndimage import gaussian_filter
 import unittest
@@ -14,8 +18,7 @@ def compute_transmittance(matgrid_symmetry=False):
         matgrid_size = mp.Vector3(2,2,0)
         matgrid_resolution = 2*resolution
 
-        Nx = int(matgrid_resolution*matgrid_size.x) + 1
-        Ny = int(matgrid_resolution*matgrid_size.y) + 1
+        Nx, Ny = mpa.compute_mg_dims(matgrid_size.x,matgrid_size.y,matgrid_resolution)
 
         # ensure reproducible results
         rng = np.random.RandomState(2069588)
@@ -90,8 +93,7 @@ def compute_resonant_mode(res,default_mat=False):
 
         # for a fixed resolution, compute the number of grid points
         # necessary which are defined on the corners of the voxels
-        Nx = int(matgrid_resolution*matgrid_size.x) + 1
-        Ny = int(matgrid_resolution*matgrid_size.y) + 1
+        Nx, Ny = mpa.compute_mg_dims(matgrid_size.x,matgrid_size.y,matgrid_resolution)
 
         x = np.linspace(-0.5*matgrid_size.x,0.5*matgrid_size.x,Nx)
         y = np.linspace(-0.5*matgrid_size.y,0.5*matgrid_size.y,Ny)