Skip to content

Instantly share code, notes, and snippets.

@oskooi
Last active May 9, 2023 21:12
Show Gist options
  • Save oskooi/fb715bb4faf31c8a02836bc235ebb7ea to your computer and use it in GitHub Desktop.
Save oskooi/fb715bb4faf31c8a02836bc235ebb7ea to your computer and use it in GitHub Desktop.
"""
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}")
@oskooi
Copy link
Author

oskooi commented May 9, 2023

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(

@smartalecH
Copy link

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment