Checks that the constraint function of a periodic design region is the same whether or not a test feature crosses its edge and wraps around back into the design region.

test_constraint_function_periodic.py

"""
Checks that the constraint function of a periodic design region
produces the same value whether or not a test feature crosses its
edge and wraps around back into the design region.
"""

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 = 3.0, 3.0
Nx, Ny = int(sx*res), int(sy*res)

def constraint_solid_void(min_length_scale: float,
                          d: np.ndarray) -> Tuple[float, float]:
    """Evaluates 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 2d array.                                                            
                                                                                                                  
    Returns:                                                                                                      
        Constraint function for solid and void regions of the design weights.                                     
    """
    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,
        (0, 1)
    )

    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, 1, (0, 1))
    M2 = lambda a: mpa.constraint_void(a, c0, eta_d, filter_f, threshold_f, 1, (0, 1))

    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.                                                                                    
        fname: filename of the PNG image for saving.                                                              
    """
    fig, ax = plt.subplots(ncols=2)
    extent = [x[0], x[-1], y[0], y[-1]]
    ax[0].imshow(d, extent=extent, cmap='binary')
    ax[0].set_xlabel('x')
    ax[0].set_ylabel('y')
    ax[0].set_title('solid')
    ax[1].imshow(1-d, extent=extent, cmap='binary')
    ax[1].set_xlabel('x')
    ax[1].set_ylabel('y')
    ax[1].set_title('void')
    fig.subplots_adjust(wspace=0.3)
    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)

    # Case 1: circle does not cross the edge of the design region                                                 
    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,
    )

    # test 1: feasible constraint                                                                                 
    constraint_solid, constraint_void = constraint_solid_void(diam + 0.5, d)
    print(f"test1-solid:, {constraint_solid:.6g} (should be < ~0)")
    print(f"test1-void:,  {constraint_void:.6g} (should be < ~0)")

    # test 2: infeasible constraint                                                                               
    constraint_solid, constraint_void = constraint_solid_void(diam - 0.5, d)
    print(f"test2-solid:, {constraint_solid:.6g} (should be > 0)")
    print(f"test2-void:,  {constraint_void:.6g} (should be > 0)")

    plot_design_weights(d, 'Case1_circle_periodic')


    # Case 2: circle crosses the edge of the design region                                                        
    cxs = (0.5*sx, -0.5*sx)
    cys = (0.5*sy, -0.5*sy)
    d = np.zeros((Nx, Ny))
    for cx in cxs:
        for cy in cys:
            d += np.where(
                np.abs(xv-cx)**2 + np.abs(yv-cy)**2 < (0.5*diam)**2,
                1.0,
                0.0,
            )

    # test 1: feasible constraint                                                                                 
    constraint_solid, constraint_void = constraint_solid_void(diam + 0.5, d)
    print(f"test1-solid:, {constraint_solid:.6g} (should be < ~0)")
    print(f"test1-void:,  {constraint_void:.6g} (should be < ~0)")

    # test 2: infeasible constraint                                                                               
    constraint_solid, constraint_void = constraint_solid_void(diam - 0.5, d)
    print(f"test2-solid:, {constraint_solid:.6g} (should be > 0)")
    print(f"test2-void:,  {constraint_void:.6g} (should be > 0)")

    plot_design_weights(d, 'Case2_circle_periodic')