# Notebook 03: Inverse design parameterization

This notebook will introduce a few basic parameterization concepts for inverse design. The same mode converter device concept as in the previous notebook will be used here as an example.

*Parameterization* refers to how we are representing our device. In the previous notebook, we simply used a 2D array of permittivity values, $\epsilon_r[x,y]$.
However, we observed that this led to continuously varying features, which is not ideal for fabricating devices generated by our optimization. 
In this notebook, we will use a modified parameterization of the device to encourage more *binarized* features in the optimized device.
However, as we will see, we will still need to tune various hyperparameters in order to get desirable results.

We begin by importing the necessary python packages:


```python
import numpy as np
import autograd.numpy as npa

import copy

import matplotlib as mpl
mpl.rcParams['figure.dpi']=100

import matplotlib.pylab as plt

from autograd.scipy.signal import convolve as conv
from skimage.draw import circle

import ceviche
from ceviche import fdfd_ez, jacobian
from ceviche.optimizers import adam_optimize
from ceviche.modes import insert_mode

import collections
# Create a container for our slice coords to be used for sources and probes
Slice = collections.namedtuple('Slice', 'x y')
```

## Introduction: device parameterization

Here, our device will be parameterized by a 2D array of density values (a density distribution, referred to as $\rho[x,y]$).
This part of the parameterization actually looks very similar to the 2D array of permitivitty values we used in the previous notebook. 
However, the difference is that we then definite two relatively simple operators that will be applied to this density. They are:

 - **Blur operator**: A 2D convolution blur filter with a configurable radius; can be applied a configurable number of times
 - **Projection operator**: A sigmoid-like nonlinear function for *binarizing* the output materials; has tunable slope and transition region

After these two operations are applied, the resulting density distribution, $\tilde{\rho}[x,y]$, will be a 2D array with values varying between `0.0` and `1.0`. A value of `0.0` represents a *background* material, while a value of `1.0` represents a *foreground* material.
The permitivitty distribution, which represents the optical description of our device, is then constructed as:
\begin{equation}
    \epsilon_r[x,y] = \epsilon_{\text{min}} + \left(\epsilon_{\text{max}} - \epsilon_{\text{min}} \right) \tilde{\rho}{[x,y]}
\end{equation}

We start by defining our parameterization operators:


```python
# Projection that drives rho towards a "binarized" design with values either 0 or 1 
def operator_proj(rho, eta=0.5, beta=100, N=1):
    """Density projection
    eta     : Center of the projection between 0 and 1
    beta    : Strength of the projection
    N       : Number of times to apply the projection
    """
    for i in range(N):
        rho =  npa.divide(npa.tanh(beta * eta) + npa.tanh(beta * (rho - eta)), 
                          npa.tanh(beta * eta) + npa.tanh(beta * (1 - eta)))
    
    return rho

# Blurring filter that results in smooth features of the structure
# First we define a function to create the kernel
def _create_blur_kernel(radius):
    """Helper function used below for creating the conv kernel"""
    rr, cc = circle(radius, radius, radius+1)
    kernel = np.zeros((2*radius+1, 2*radius+1), dtype=np.float)
    kernel[rr, cc] = 1
    return  kernel/kernel.sum()

# Then we define the function to apply the operation
def operator_blur(rho, radius=2, N=1):
    """Blur operator implemented via two-dimensional convolution
    radius    : Radius of the circle for the conv kernel filter
    N         : Number of times to apply the filter
    
    Note that depending on the radius, the kernel is not always a
    perfect circle due to "pixelation" / stair casing
    """
    
    kernel = _create_blur_kernel(radius)
    
    for i in range(N):
        # For whatever reason HIPS autograd doesn't support 'same' mode, so we need to manually crop the output
        rho = conv(rho, kernel, mode='full')[radius:-radius,radius:-radius]
    
    return rho
```

### Visualization of the blur

Below we visualize the blur operator applied to a random 2D array:


```python
# Specify a range of values for the blur radius
blur_radii = [2, 3, 4, 5, 6]

# Number of times to apply the blur filter
N_blur = 1

# Define a random 2D array to test our blur filter on
rho = np.random.rand(50, 50)

# Create a figure with panels to plot into
fig, axs = plt.subplots(2, len(blur_radii)+1,figsize=(12,4), constrained_layout=True)

# First, plot a dummy conv filter
axs[0,0].set_title('Initial random image')
axs[0,0].imshow(np.zeros((1,1)), vmin=0, cmap='Greys')
# And also plot the initial random image
axs[1,0].imshow(rho)

# Now, loop over the blur radii and visualize each result
for i, radius in enumerate(blur_radii):
    kernel = _create_blur_kernel(radius)
    
    kernel_pad = np.pad(kernel, (2+np.max(blur_radii)-radius, 2+np.max(blur_radii)-radius),
                        'constant', constant_values=(0, 0))
    axs[0,i+1].imshow(kernel_pad, vmin=0)
    
    axs[0,i+1].set_title('radius = %d' % radius)

    rho_p = operator_blur(rho, radius=radius, N=N_blur)
    axs[1,i+1].imshow(rho_p)

axs[0,0].set_yticks([])
axs[0,0].set_xticks([])
axs[0,0].set_ylabel('Conv filter')
axs[1,0].set_ylabel('Image')

fig.align_labels()
plt.show()
```

### Visualization of the projection

Below we visualize the projection operator:


```python
rho = np.linspace(-0.5, +1.5, 999)

# Visualize different values of the projection strength
plt.figure()
for beta in [5, 10, 50, 100]:
    plt.plot(rho, operator_proj(rho, beta=beta), label=r"$\beta$ = %d" % beta)

plt.xlabel(r"Input $\rho$")
plt.ylabel(r"Projected $\hat{\rho}$")
plt.legend()
plt.show()
```

### Visualization of blur + projection

Finally, we visualize the combined projection and blur operation. Notice how the larger blur radius leads to larger and smoother features.


```python
# Specify a range of values for the blur radius
blur_radii = [2, 3, 4, 5, 6, 10]

# Number of times to apply the blur filter
N_blur = 1

# Number of times to apply the projection operator
N_proj = 1

# Specify beta value to use for projection
beta = 200

# Specify eta value to use for projection
eta = 0.5

# Define a random 2D array to test our blur filter on
rho = np.random.rand(100, 100)

# Create a figure with panels to plot into
fig, axs = plt.subplots(2, len(blur_radii)+1,figsize=(12,4), constrained_layout=True)

# First, plot a dummy conv filter
axs[0,0].set_title('Initial random image')
axs[0,0].imshow(np.zeros((1,1)), vmin=0, cmap='Greys')
# And also plot the initial random image
axs[1,0].imshow(rho)

# Now, loop over the blur radii and visualize each result
for i, radius in enumerate(blur_radii):
    kernel = _create_blur_kernel(radius)
    
    kernel_pad = np.pad(kernel, (2+np.max(blur_radii)-radius, 2+np.max(blur_radii)-radius),
                        'constant', constant_values=(0, 0))
    axs[0,i+1].imshow(kernel_pad, vmin=0)
    
    axs[0,i+1].set_title('radius = %d' % radius)

    rho_p = operator_blur(rho, radius=radius, N=N_blur)
    rho_p = operator_proj(rho_p, beta=beta, eta=eta, N=N_proj)
    axs[1,i+1].imshow(rho_p)

axs[0,0].set_yticks([])
axs[0,0].set_xticks([])
axs[0,0].set_ylabel('Conv filter')
axs[1,0].set_ylabel('Image')

fig.align_labels()
plt.show()
```

---

## Simulation and optimization parameters


```python
# Angular frequency of the source in 1/s
omega=2*np.pi*200e12
# Spatial resolution in meters
dl=40e-9
# Number of pixels in x-direction
Nx=120
# Number of pixels in y-direction
Ny=120
# Number of pixels in the PMLs in each direction
Npml=20
# Minimum value of the relative permittivity
epsr_min=1.0
# Maximum value of the relative permittivity
epsr_max=12.0
# Radius of the smoothening features
blur_radius=5
# Number of times to apply the blur
N_blur=1
# Strength of the binarizing projection
beta=500.0
# Middle point of the binarizing projection
eta=0.5
# Number of times to apply the blur
N_proj=1
# Space between the PMLs and the design region (in pixels)
space=10
# Width of the waveguide (in pixels)
wg_width=12
# Length in pixels of the source/probe slices on each side of the center point
space_slice=8

# Number of epochs in the optimization 
Nsteps=50
# Step size for the Adam optimizer
step_size=5e-3
```

## Utility functions

As in **Notebook 02**, here we define several utility functions for constructing the simulation domain, feed waveguides, design region, input source, and output probe,


```python
def init_domain(Nx, Ny, Npml, space=10, wg_width=10, space_slice=5):
    """Initializes the domain and design region

    space       : The space between the PML and the structure
    wg_width    : The feed and probe waveguide width
    space_slice : The added space for the probe and source slices
    """
    rho = np.zeros((Nx, Ny))  
    bg_rho = np.zeros((Nx, Ny))  
    design_region = np.zeros((Nx, Ny))
    
    # Input waveguide
    bg_rho[0:int(Npml+space),int(Ny/2-wg_width/2):int(Ny/2+wg_width/2)] = 1

    # Input probe slice
    input_slice = Slice(x=np.array(Npml+1), 
        y=np.arange(int(Ny/2-wg_width/2-space_slice), int(Ny/2+wg_width/2+space_slice)))
    
    # Output waveguide 
    bg_rho[int(Nx-Npml-space)::,int(Ny/2-wg_width/2):int(Ny/2+wg_width/2)] = 1
    
    # Output probe slice
    output_slice = Slice(x=np.array(Nx-Npml-1), 
                         y=np.arange(int(Ny/2-wg_width/2-space_slice), int(Ny/2+wg_width/2+space_slice)))
    
    design_region[Npml+space:Nx-Npml-space, Npml+space:Ny-Npml-space] = 1
    
    # Const init
    rho[Npml+space:Nx-Npml-space, Npml+space:Ny-Npml-space] = 0.5
    # Random init
    # rho = design_region * np.random.rand(Nx, Ny)

    return rho, bg_rho, design_region, input_slice, output_slice

def viz_sim(epsr, source, slices=[]):
    """Solve and visualize a simulation with permittivity 'epsr'
    """
    simulation = fdfd_ez(omega, dl, epsr, [Npml, Npml])
    Hx, Hy, Ez = simulation.solve(source)
    
    fig, ax = plt.subplots(1, 2, constrained_layout=True, figsize=(6,3))
    ceviche.viz.real(Ez, outline=epsr, ax=ax[0], cbar=False)
    for sl in slices:
        ax[0].plot(sl.x*np.ones(len(sl.y)), sl.y, 'b-')
    ceviche.viz.abs(epsr, ax=ax[1], cmap='Greys');
    plt.show()
    
    return (simulation, ax)

def mask_combine_rho(rho, bg_rho, design_region):
    """Utility function for combining the design region rho and the background rho
    """
    return rho*design_region + bg_rho*(design_region==0).astype(np.float)

def epsr_parameterization(rho, bg_rho, design_region, radius=2, N_blur=1, beta=100, eta=0.5, N_proj=1):
    """Defines the parameterization steps for constructing rho
    """
    # Combine rho and bg_rho; Note: this is so the subsequent blur sees the waveguides
    rho = mask_combine_rho(rho, bg_rho, design_region)
    
    rho = operator_blur(rho, radius=radius, N=N_blur)
    rho = operator_proj(rho, beta=beta, eta=eta, N=N_proj)
    
    # Final masking undoes the blurring of the waveguides
    rho = mask_combine_rho(rho, bg_rho, design_region)
    
    return epsr_min + (epsr_max-epsr_min) * rho

def mode_overlap(E1, E2):
    """Defines an overlap integral between the sim field and desired field
    """
    return npa.abs(npa.sum(npa.conj(E1)*E2))
```

## Simulate the initial structure




```python
# Initialize the parametrization rho and the design region
rho, bg_rho, design_region, input_slice, output_slice = \
    init_domain(Nx, Ny, Npml, space=space, wg_width=wg_width, space_slice=space_slice)

# Compute the permittivity from rho_init, including blurring and projection
epsr_init = epsr_parameterization(rho, bg_rho, design_region, \
                                  radius=blur_radius, N_blur=N_blur, beta=beta, eta=eta, N_proj=N_proj)

# Setup source
source = insert_mode(omega, dl, input_slice.x, input_slice.y, epsr_init, m=1)

# Setup probe
probe = insert_mode(omega, dl, output_slice.x, output_slice.y, epsr_init, m=2)

# Simulate initial device
simulation, ax = viz_sim(epsr_init, source, slices=[input_slice, output_slice])

# get normalization factor (field overlap before optimizing)
_, _, Ez = simulation.solve(source)
E0 = mode_overlap(Ez, probe)
```

## Run the mode converter optimization


```python
def objective(rho):
    """Objective function called by optimizer
    
    1) Takes the density distribution as input
    2) Constructs epsr
    2) Runs the simulation
    3) Returns the overlap integral between the output wg field 
       and the desired mode field
    """
    rho = rho.reshape((Nx, Ny))
    epsr = epsr_parameterization(rho, bg_rho, design_region, \
                             radius=blur_radius, N_blur=N_blur, beta=beta, eta=eta, N_proj=N_proj)
    simulation.eps_r = epsr
    _, _, Ez = simulation.solve(source)
    return mode_overlap(Ez, probe) / E0

# Compute the gradient of the objective function using revere-mode differentiation
objective_jac = jacobian(objective, mode='reverse')

# Maximize the objective function using an ADAM optimizer
(rho_optimum, loss) = adam_optimize(objective, rho.flatten(), objective_jac,
                         Nsteps=Nsteps, direction='max', step_size=step_size)

# Simulate optimal device
rho_optimum = rho_optimum.reshape((Nx, Ny))
epsr = epsr_parameterization(rho_optimum, bg_rho, design_region, \
                                  radius=blur_radius, N_blur=N_blur, beta=beta, eta=eta, N_proj=N_proj)

viz_sim(epsr, source, slices=[input_slice, output_slice]);
```

## Visualizing the parameterization steps

We can take a look what each step of our parameterization actually looks like. Below we visualize the raw density distribution, the blurred density distribution, and finally, the projected device density.


```python
fig, axs = plt.subplots(1,3,constrained_layout=True,sharex=True,sharey=True, figsize=(8.5,3))
Z = mask_combine_rho(rho_optimum, bg_rho, design_region)
ceviche.viz.abs(Z, cmap='Greys', ax=axs[0])
axs[0].set_xlabel('')
axs[0].set_xticks([])
axs[0].set_ylabel('')
axs[0].set_yticks([])
axs[0].set_title('Raw input density')

Z = mask_combine_rho(rho_optimum, bg_rho, design_region)
Z = operator_blur(Z, radius=blur_radius, N=N_blur)
# Z = operator_proj(Z, beta=beta, eta=eta)
Z = mask_combine_rho(Z, bg_rho, design_region)
ceviche.viz.abs(Z, cmap='Greys', ax=axs[1])
axs[1].set_xlabel('')
axs[1].set_ylabel('')
axs[1].set_title('Blurred density')

Z = rho_optimum
Z = mask_combine_rho(Z, bg_rho, design_region)
Z = operator_blur(Z, radius=blur_radius, N=N_blur)
Z = operator_proj(Z, beta=beta, eta=eta, N=N_proj)
Z = mask_combine_rho(Z, bg_rho, design_region)
ceviche.viz.abs(Z, cmap='Greys', ax=axs[2])
axs[2].set_xlabel('')
axs[2].set_ylabel('')
axs[2].set_title('Projected density (final structure)');
```

## Penalizing the amount of material

We notice that in the optimized device shown above that there seems to be a lot of unnecessary material in the device. Here, we will try to eliminate some of this extra material by adding a penalty term to the objective function.

To add the penalty term, we only need to modify the objective function:


```python
def objective(rho, penalty_weight=1e8): 
    rho = rho.reshape((Nx, Ny))
    epsr = epsr_parameterization(rho, bg_rho, design_region, \
                             radius=blur_radius, N_blur=N_blur, beta=beta, eta=eta, N_proj=N_proj)
    simulation.eps_r = epsr
    _, _, Ez = simulation.solve(source)
    
    # This penalty term is directly proportional to the material area
    # penalty = penalty_weight * (design_region*(epsr-1)).sum() # penalty_weight = 1e-10
    
    # This penalty term is the L2-norm of the raw density
    penalty = penalty_weight * npa.linalg.norm(rho)
    
    return mode_overlap(Ez, probe) / E0 - penalty
```


```python
# Initialize the parametrization rho and the design region
rho, bg_rho, design_region, input_slice, output_slice = \
    init_domain(Nx, Ny, Npml, space=space, wg_width=wg_width, space_slice=space_slice)

# Compute the permittivity from rho_init, including blurring and projection
epsr_init = epsr_parameterization(rho, bg_rho, design_region, \
                                  radius=blur_radius, N_blur=N_blur, beta=beta, eta=eta, N_proj=N_proj)

# Setup source
source = insert_mode(omega, dl, input_slice.x, input_slice.y, epsr_init, m=1)

# Setup probe
probe = insert_mode(omega, dl, output_slice.x, output_slice.y, epsr_init, m=2)

# Simulate initial device
viz_sim(epsr_init, source, slices=[input_slice, output_slice])

# Run optimization
objective_jac = jacobian(objective, mode='reverse')
(rho_optimum, loss) = adam_optimize(objective, rho.flatten(), objective_jac,
                         Nsteps=Nsteps, direction='max', step_size=step_size)

# Simulate optimal device
rho_optimum = rho_optimum.reshape((Nx, Ny))
epsr_pen = epsr_parameterization(rho_optimum, bg_rho, design_region, \
                                  radius=blur_radius, N_blur=N_blur, beta=beta, eta=eta, N_proj=N_proj)
viz_sim(epsr_pen, source, slices=[input_slice, output_slice]);
```

We can then quantify how much our penalization impacted the amount of material used in the final design. 


```python
# Calculate areas
def calc_design_area(design_region, epsr_max, epsr_min, epsr, dl):
    """Computes the area of material used in the design region
    """
    A = ((epsr-epsr_min)/epsr_max * design_region).sum() * dl**2 * 1e12
    A_design_region = design_region.sum() * (dl)**2 * 1e12
    return A, A_design_region

_, A_design_region = calc_design_area(design_region, epsr_max, epsr_min, epsr, dl)
A_original, _ = calc_design_area(design_region, epsr_max, epsr_min, epsr, dl)
A_pen, _ = calc_design_area(design_region, epsr_max, epsr_min, epsr_pen, dl)

# Print summary
print('Design region area: %.2f um^2' % A_design_region)
print('Unpenalized design area: %.2f um^2' % A_original)
print('Penalized design area: %.2f um^2' % A_pen)
print('---')
print('Improvement: %.2f%%' % (100*(1-A_pen/A_original)))
```

    Design region area: 5.76 um^2
    Unpenalized design area: 1.96 um^2
    Penalized design area: 1.03 um^2
    ---
    Improvement: 47.32%



```python

```


```python

```

e41605cafc168bfa6a48ebfeacfad8a8ee0d112d
ipynb
Jupyter Notebook
03_Invdes_parameterization.ipynb
fancompute/workshop-invdesign
200eaa0abc3f691137e228e98ebb62446015ec38
MIT
2019-11-22T18:22:21.000Z
2022-03-15T15:38:08.000Z
03_Invdes_parameterization.ipynb
Ydeh22/workshop-invdesign
200eaa0abc3f691137e228e98ebb62446015ec38
MIT
2019-12-14T16:57:42.000Z
2021-04-01T05:41:30.000Z
03_Invdes_parameterization.ipynb
Ydeh22/workshop-invdesign
200eaa0abc3f691137e228e98ebb62446015ec38
MIT
2019-11-23T19:37:37.000Z
2022-03-22T22:30:20.000Z
Qwen/Qwen-72B
1. YES
2. YES