test_constraint_function.py

"""                                                                             
Checks that the constraint function with constraint L applied to a circle
with diameter d has a return value based on the convention that a violation
(L > d) is >=0 and a non-violation is ~0. 
"""

from typing import Tuple

import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
import meep.adjoint as mpa
import numpy as np

# properties of design region                                                   
res = 100
sx, sy = 5.0, 5.0
Nx, Ny = int(sx*res), int(sy*res)

def constraint_solid_void(min_length_scale: float,
                          d: np.ndarray) -> Tuple[float, float]:
    """Returns the value of the constraint function for the solid and void  
       regions of the design weights.                                                          
                                                                                
    Args:                                                                       
        min_length_scale: the minimum length scale (line width and spacing) of  
            the design region to check for a violation.                         
        d: weights of the design region as a flattened 1d array.
    """
    beta = 32  # projection bias                                                
    eta_i = 0.5  # blueprint design field thresholding point in [0, 1]          
    eta_e = 0.75  # erosion design field thresholding point in [0, 1]           
    eta_d = 1 - eta_e  # dilation design field thresholding point in [0, 1]     

    filter_radius = mpa.get_conic_radius_from_eta_e(min_length_scale, eta_e)

    # filter radius should be equivalent to minimum length scale                
    assert filter_radius == min_length_scale

    threshold_f = lambda a: mpa.tanh_projection(a, beta, eta_i)

    filter_f = lambda a: mpa.conic_filter(
        a.reshape(Nx, Ny),
        filter_radius,
        sx,
        sy,
        res,
    )

    prefac = 1e7  # hyperparameter 1                                            
    pwr = 4.0  # hyperparameter 2                                               
    c0 = prefac * (filter_radius * 1 / res) ** pwr

    M1 = lambda a: mpa.constraint_solid(
        a, c0, eta_e, filter_f, threshold_f, res,
    )
    M2 = lambda a: mpa.constraint_void(
        a, c0, eta_d, filter_f, threshold_f, res,
    )

    return M1(d), M2(1-d)
  
  def plot_design_weights(d: np.ndarray, fname: str):
    """Saves a plot of the design weights.                                      
                                                                                
    Args:                                                                       
        d: the design weights as a flattened 1d array.                                                  
        fname: filename of the PNG image for saving.                            
    """
    fig, ax = plt.subplots(ncols=2)
    extent = [x[0], x[-1], y[0], y[-1]]
    ax.imshow(d.reshape(Nx, Ny), extent=extent, cmap='binary')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_title('solid')
    fig.savefig(fname + '.png', bbox_inches='tight', dpi=150)

if __name__ == "__main__":
    # grid coordinates of design region                                         
    x = np.linspace(-0.5*sx, 0.5*sx, Nx)
    y = np.linspace(-0.5*sy, 0.5*sy, Ny)
    xv, yv = np.meshgrid(x, y)

    diam = 1.0
    cx, cy = 0.0, 0.0
    d = np.where(
        np.abs(xv-cx)**2 + np.abs(yv-cy)**2 < (0.5*diam)**2,
        1.0,
        0.0,
    ).flatten()

    plot_design_weights(d, f'circle_diam{diam}')

    lc = 0.8367
    c_solid, c_void = constraint_solid_void(lc, d)
    print(f"data:, {lc:.4f}, {c_solid:.6f}, {c_void:.6f}")
    
    length_constraint = np.linspace(0.1, 2.0, 50)
    for lc in length_constraint:
        c_solid, c_void = constraint_solid_void(lc, d)
        print(f"data:, {lc:.4f}, {c_solid:.6f}, {c_void:.6f}")

diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py
index 2637a747..812b5fe3 100644
--- a/python/adjoint/filters.py
+++ b/python/adjoint/filters.py
@@ -6,6 +6,10 @@ import numpy as np
 from autograd import numpy as npa
 from scipy import signal, special
 
+import matplotlib
+matplotlib.use('agg')
+import matplotlib.pyplot as plt
+
 
 def _centered(arr, newshape):
     """Helper function that reformats the padded array of the fft filter operation.
@@ -845,6 +849,14 @@ def indicator_solid(x, c, filter_f, threshold_f, resolution, periodic_axes=None)
     filtered_field = filter_f(x)
     design_field = threshold_f(filtered_field)
 
+    fig, ax = plt.subplots()
+    extent = [-2.5, 2.5, -2.5, 2.5]
+    ax.imshow(design_field.reshape(500, 500), extent=extent, cmap='binary')
+    ax.set_xlabel('x')
+    ax.set_ylabel('y')
+    ax.set_title('design_field from mpa.indicator_solid')
+    fig.savefig('design_field_indicator_solid.png', bbox_inches='tight', dpi=150)
+
     if periodic_axes is None:
         gradient_filtered_field = npa.gradient(filtered_field)
     else:
@@ -864,7 +876,18 @@ def indicator_solid(x, c, filter_f, threshold_f, resolution, periodic_axes=None)
         raise ValueError(
             "The gradient fields must be 2 dimensional. Check input array and filter functions."
         )
-    return design_field * npa.exp(-c * grad_mag)
+
+    indicator_func = design_field * npa.exp(-c * grad_mag)
+    fig, ax = plt.subplots()
+    extent = [-2.5, 2.5, -2.5, 2.5]
+    pos = ax.imshow(indicator_func.reshape(500, 500), extent=extent, cmap='inferno_r')
+    ax.set_xlabel('x')
+    ax.set_ylabel('y')
+    ax.set_title('indicator function from mpa.indicator_solid')
+    fig.colorbar(pos, ax=ax)
+    fig.savefig('indicator_func_indicator_solid.png', bbox_inches='tight', dpi=150)
+
+    return indicator_func
 
 
 def constraint_solid(

@oskooi see my fork for changes as described in the original issue.