From 23831d47c453469c456fa47e16f1c436e80208a2 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Thu, 10 Dec 2020 15:10:16 -0500 Subject: [PATCH 01/15] update filter bugs --- python/adjoint/filters.py | 188 +++++++++++++++++++++----------------- 1 file changed, 103 insertions(+), 85 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index b2605f740..e3b9c2056 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -3,7 +3,7 @@ """ import numpy as np -from autograd import numpy as npa +from jax import numpy as npj import meep as mp from scipy import special @@ -26,47 +26,46 @@ def _centered(arr, newshape): def _edge_pad(arr, pad): # fill sides - left = npa.tile(arr[0, :], (pad[0][0], 1)) # left side - right = npa.tile(arr[-1, :], (pad[0][1], 1)) # right side - top = npa.tile(arr[:, 0], (pad[1][0], 1)).transpose() # top side - bottom = npa.tile(arr[:, -1], (pad[1][1], 1)).transpose() # bottom side) - + left = npj.tile(arr[0,:],(pad[0][0],1)) # left side + right = npj.tile(arr[-1,:],(pad[0][1],1)) # right side + top = npj.tile(arr[:,0],(pad[1][0],1)).transpose() # top side + bottom = npj.tile(arr[:,-1],(pad[1][1],1)).transpose() # bottom side) + # fill corners - top_left = npa.tile(arr[0, 0], (pad[0][0], pad[1][0])) # top left - top_right = npa.tile(arr[-1, 0], (pad[0][1], pad[1][0])) # top right - bottom_left = npa.tile(arr[0, -1], (pad[0][0], pad[1][1])) # bottom left - bottom_right = npa.tile(arr[-1, -1], - (pad[0][1], pad[1][1])) # bottom right - - out = npa.concatenate((npa.concatenate( - (top_left, top, top_right)), npa.concatenate((left, arr, right)), - npa.concatenate( - (bottom_left, bottom, bottom_right))), - axis=1) - + top_left = npj.tile(arr[0,0], (pad[0][0],pad[1][0])) # top left + top_right = npj.tile(arr[-1,0], (pad[0][1],pad[1][0])) # top right + bottom_left = npj.tile(arr[0,-1], (pad[0][0],pad[1][1])) # bottom left + bottom_right = npj.tile(arr[-1,-1], (pad[0][1],pad[1][1])) # bottom right + + out = npj.concatenate(( + npj.concatenate((top_left,top,top_right)), + npj.concatenate((left,arr,right)), + npj.concatenate((bottom_left,bottom,bottom_right)) + ),axis=1) + return out def _zero_pad(arr, pad): # fill sides - left = npa.tile(0, (pad[0][0], arr.shape[1])) # left side - right = npa.tile(0, (pad[0][1], arr.shape[1])) # right side - top = npa.tile(0, (arr.shape[0], pad[1][0])) # top side - bottom = npa.tile(0, (arr.shape[0], pad[1][1])) # bottom side - + left = npj.tile(0,(pad[0][0],arr.shape[1])) # left side + right = npj.tile(0,(pad[0][1],arr.shape[1])) # right side + top = npj.tile(0,(arr.shape[0],pad[1][0])) # top side + bottom = npj.tile(0,(arr.shape[0],pad[1][1])) # bottom side + # fill corners - top_left = npa.tile(0, (pad[0][0], pad[1][0])) # top left - top_right = npa.tile(0, (pad[0][1], pad[1][0])) # top right - bottom_left = npa.tile(0, (pad[0][0], pad[1][1])) # bottom left - bottom_right = npa.tile(0, (pad[0][1], pad[1][1])) # bottom right - - out = npa.concatenate((npa.concatenate( - (top_left, top, top_right)), npa.concatenate((left, arr, right)), - npa.concatenate( - (bottom_left, bottom, bottom_right))), - axis=1) - + top_left = npj.tile(0, (pad[0][0],pad[1][0])) # top left + top_right = npj.tile(0, (pad[0][1],pad[1][0])) # top right + bottom_left = npj.tile(0, (pad[0][0],pad[1][1])) # bottom left + bottom_right = npj.tile(0, (pad[0][1],pad[1][1])) # bottom right + + out = npj.concatenate(( + npj.concatenate((top_left,top,top_right)), + npj.concatenate((left,arr,right)), + npj.concatenate((bottom_left,bottom,bottom_right)) + ),axis=1) + return out @@ -114,30 +113,30 @@ def simple_2d_filter(x, kernel, Lx, Ly, resolution, symmetries=[]): # Perform the required reflections for symmetries if mp.X in symmetries: if kx % 2 == 1: - x = npa.concatenate((x, x[-1, :][None, :], x[::-1, :]), axis=0) + x = npj.concatenate((x,x[-1,:][None,:],x[::-1,:]), axis=0) else: - x = npa.concatenate((x, x[::-1, :]), axis=0) + x = npj.concatenate((x,x[::-1,:]), axis=0) if mp.Y in symmetries: if ky % 2 == 1: - x = npa.concatenate((x[:, ::-1], x[:, -1][:, None], x), axis=1) + x = npj.concatenate((x[:,::-1],x[:,-1][:,None],x), axis=1) else: - x = npa.concatenate((x[:, ::-1], x), axis=1) - + x = npj.concatenate((x[:,::-1],x), axis=1) + # pad the kernel and input to avoid circular convolution and # to ensure boundary conditions are met. kernel = _zero_pad(kernel, ((kx, kx), (ky, ky))) x = _edge_pad(x, ((kx, kx), (ky, ky))) # Transform to frequency domain for fast convolution - H = npa.fft.fft2(kernel) - X = npa.fft.fft2(x) - + H = npj.fft.fft2(kernel) + X = npj.fft.fft2(x) + # 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))) - + y = npj.fft.fftshift(npj.real(npj.fft.ifft2(Y))) + # Remove all the extra padding y = _centered(y, (kx, ky)) @@ -338,6 +337,9 @@ def exponential_erosion(x, radius, beta, Lx, Ly, resolution): morphological close-open and open-close filters for topology optimization. Structural and Multidisciplinary Optimization, 54(1), 15-21. ''' + + x_hat = npj.exp(beta*(1-x)) + return 1 - npj.log(cylindrical_filter(x_hat,radius,Lx,Ly,resolution).flatten()) / beta x_hat = npa.exp(beta * (1 - x)) return 1 - npa.log( @@ -375,6 +377,9 @@ def exponential_dilation(x, radius, beta, Lx, Ly, resolution): morphological close-open and open-close filters for topology optimization. Structural and Multidisciplinary Optimization, 54(1), 15-21. ''' + + x_hat = npj.exp(beta*x) + return npj.log(cylindrical_filter(x_hat,radius,Lx,Ly,resolution).flatten()) / beta x_hat = npa.exp(beta * x) return npa.log( @@ -410,6 +415,9 @@ def heaviside_erosion(x, radius, beta, Lx, Ly, resolution): optimization using nodal design variables and projection functions. International journal for numerical methods in engineering, 61(2), 238-254. ''' + + x_hat = cylindrical_filter(x,radius,Lx,Ly,resolution).flatten() + return npj.exp(-beta*(1-x_hat)) + npj.exp(-beta)*(1-x_hat) x_hat = cylindrical_filter(x, radius, Lx, Ly, resolution).flatten() return npa.exp(-beta * (1 - x_hat)) + npa.exp(-beta) * (1 - x_hat) @@ -444,6 +452,9 @@ def heaviside_dilation(x, radius, beta, Lx, Ly, resolution): optimization using nodal design variables and projection functions. International journal for numerical methods in engineering, 61(2), 238-254. ''' + + x_hat = cylindrical_filter(x,radius,Lx,Ly,resolution).flatten() + return 1 - npj.exp(-beta*x_hat) + npj.exp(-beta)*x_hat x_hat = cylindrical_filter(x, radius, Lx, Ly, resolution).flatten() return 1 - npa.exp(-beta * x_hat) + npa.exp(-beta) * x_hat @@ -477,9 +488,8 @@ def geometric_erosion(x, radius, alpha, Lx, Ly, resolution): [1] Svanberg, K., & Svärd, H. (2013). Density filters for topology optimization based on the Pythagorean means. Structural and Multidisciplinary Optimization, 48(5), 859-875. ''' - x_hat = npa.log(x + alpha) - return npa.exp(cylindrical_filter(x_hat, radius, Lx, Ly, - resolution)).flatten() - alpha + x_hat = npj.log(x + alpha) + return npj.exp(cylindrical_filter(x_hat,radius,Lx,Ly,resolution)).flatten() - alpha def geometric_dilation(x, radius, alpha, Lx, Ly, resolution): @@ -511,10 +521,8 @@ def geometric_dilation(x, radius, alpha, Lx, Ly, resolution): Pythagorean means. Structural and Multidisciplinary Optimization, 48(5), 859-875. ''' - x_hat = npa.log(1 - x + alpha) - return -npa.exp(cylindrical_filter(x_hat, radius, Lx, Ly, - resolution)).flatten() + alpha + 1 - + x_hat = npj.log(1 - x + alpha) + return -npj.exp(cylindrical_filter(x_hat,radius,Lx,Ly,resolution)).flatten() + alpha + 1 def harmonic_erosion(x, radius, alpha, Lx, Ly, resolution): ''' Performs a harmonic erosion operation. @@ -612,12 +620,8 @@ def tanh_projection(x, beta, eta): [1] Wang, F., Lazarov, B. S., & Sigmund, O. (2011). On projection methods, convergence and robust formulations in topology optimization. Structural and Multidisciplinary Optimization, 43(6), 767-784. ''' - - return (npa.tanh(beta * eta) + - npa.tanh(beta * - (x - eta))) / (npa.tanh(beta * eta) + npa.tanh(beta * - (1 - eta))) - + + return (npj.tanh(beta*eta) + npj.tanh(beta*(x-eta))) / (npj.tanh(beta*eta) + npj.tanh(beta*(1-eta))) def heaviside_projection(x, beta, eta): '''Projection filter that thresholds the input parameters between 0 and 1. @@ -641,12 +645,10 @@ def heaviside_projection(x, beta, eta): [1] Lazarov, B. S., Wang, F., & Sigmund, O. (2016). Length scale and manufacturability in density-based topology optimization. Archive of Applied Mechanics, 86(1-2), 189-218. ''' - - case1 = eta * npa.exp(-beta * (eta - x) / eta) - (eta - x) * npa.exp(-beta) - case2 = 1 - (1 - eta) * npa.exp(-beta * (x - eta) / - (1 - eta)) - (eta - x) * npa.exp(-beta) - return npa.where(x < eta, case1, case2) - + + case1 = eta*npj.exp(-beta*(eta-x)/eta) - (eta-x)*npj.exp(-beta) + case2 = 1 - (1-eta)*npj.exp(-beta*(x-eta)/(1-eta)) - (eta-x)*npj.exp(-beta) + return npj.where(x < eta,case1,case2) ''' # ------------------------------------------------------------------------------------ # @@ -791,14 +793,11 @@ def indicator_solid(x, c, filter_f, threshold_f, resolution): filtered_field = filter_f(x) design_field = threshold_f(filtered_field) - gradient_filtered_field = npa.gradient(filtered_field) - grad_mag = (gradient_filtered_field[0] * - resolution)**2 + (gradient_filtered_field[1] * resolution)**2 + gradient_filtered_field = npj.gradient(filtered_field) + grad_mag = (gradient_filtered_field[0]*resolution) ** 2 + (gradient_filtered_field[1]*resolution) ** 2 if grad_mag.ndim != 2: - raise ValueError( - "The gradient fields must be 2 dimensional. Check input array and filter functions." - ) - I_s = design_field * npa.exp(-c * grad_mag) + raise ValueError("The gradient fields must be 2 dimensional. Check input array and filter functions.") + I_s = design_field * npj.exp(-c * grad_mag) return I_s @@ -835,9 +834,8 @@ def constraint_solid(x, c, eta_e, filter_f, threshold_f, resolution): ''' filtered_field = filter_f(x) - I_s = indicator_solid(x.reshape(filtered_field.shape), c, filter_f, - threshold_f, resolution).flatten() - return npa.mean(I_s * npa.minimum(filtered_field.flatten() - eta_e, 0)**2) + I_s = indicator_solid(x.reshape(filtered_field.shape),c,filter_f,threshold_f,resolution).flatten() + return npj.mean(I_s * npj.minimum(filtered_field.flatten()-eta_e,0)**2) def indicator_void(x, c, filter_f, threshold_f, resolution): @@ -869,14 +867,11 @@ def indicator_void(x, c, filter_f, threshold_f, resolution): filtered_field = filter_f(x).reshape(x.shape) design_field = threshold_f(filtered_field) - gradient_filtered_field = npa.gradient(filtered_field) - grad_mag = (gradient_filtered_field[0] * - resolution)**2 + (gradient_filtered_field[1] * resolution)**2 + gradient_filtered_field = npj.gradient(filtered_field) + grad_mag = (gradient_filtered_field[0]*resolution) ** 2 + (gradient_filtered_field[1]*resolution) ** 2 if grad_mag.ndim != 2: - raise ValueError( - "The gradient fields must be 2 dimensional. Check input array and filter functions." - ) - return (1 - design_field) * npa.exp(-c * grad_mag) + raise ValueError("The gradient fields must be 2 dimensional. Check input array and filter functions.") + return (1 - design_field) * npj.exp(-c * grad_mag) def constraint_void(x, c, eta_d, filter_f, threshold_f, resolution): @@ -912,10 +907,8 @@ def constraint_void(x, c, eta_d, filter_f, threshold_f, resolution): ''' filtered_field = filter_f(x) - I_v = indicator_void(x.reshape(filtered_field.shape), c, filter_f, - threshold_f, resolution).flatten() - return npa.mean(I_v * npa.minimum(eta_d - filtered_field.flatten(), 0)**2) - + I_v = indicator_void(x.reshape(filtered_field.shape),c,filter_f,threshold_f,resolution).flatten() + return npj.mean(I_v * npj.minimum(eta_d-filtered_field.flatten(),0)**2) def gray_indicator(x): '''Calculates a measure of "grayness" according to [1]. @@ -937,4 +930,29 @@ def gray_indicator(x): [1] Lazarov, B. S., Wang, F., & Sigmund, O. (2016). Length scale and manufacturability in density-based topology optimization. Archive of Applied Mechanics, 86(1-2), 189-218. ''' - return npa.mean(4 * x.flatten() * (1 - x.flatten())) * 100 + return npj.mean(4 * x.flatten() * (1-x.flatten())) * 100 + +def F_diff_open_close(x,f_open,f_close,beta=64): + ''' + + ''' + return npj.mean(tanh_projection(npj.abs(f_close(x.flatten())-f_open(x.flatten()))),beta,0.5) + +def M_diff_open_close(x,f_open,f_close,p=1): + ''' + + ''' + return npj.mean(npj.abs(f_close(x.flatten())-f_open(x.flatten()))**p) + +def M_diff_identity_open(x,f_open,p=1): + ''' + + ''' + return npj.mean(npj.abs(x.flatten()-f_open(x.flatten()))**p) + +def M_diff_identity_close(x,f_close,p=1): + ''' + + ''' + return npj.mean(npj.abs(x.flatten()-f_close(x.flatten()))**p) + From 73803d13d0059788a41996fd46cec1c6f7f29175 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Thu, 10 Dec 2020 15:45:39 -0500 Subject: [PATCH 02/15] rebase --- python/adjoint/optimization_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py index 6ad79b2f0..80485a20d 100644 --- a/python/adjoint/optimization_problem.py +++ b/python/adjoint/optimization_problem.py @@ -1,6 +1,6 @@ import meep as mp import numpy as np -from autograd import grad, jacobian +from jax import grad, jacobian from collections import namedtuple Grid = namedtuple('Grid', ['x', 'y', 'z', 'w']) From e088f711284c2b52052712ab7ba8d35b58e9721a Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Tue, 15 Dec 2020 13:45:59 -0500 Subject: [PATCH 03/15] normalize --- python/adjoint/filters.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index e3b9c2056..bfd108d2f 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -942,17 +942,17 @@ def M_diff_open_close(x,f_open,f_close,p=1): ''' ''' - return npj.mean(npj.abs(f_close(x.flatten())-f_open(x.flatten()))**p) + return npj.mean(npj.abs(f_close(x.flatten())-f_open(x.flatten()))**p) * 100 def M_diff_identity_open(x,f_open,p=1): ''' ''' - return npj.mean(npj.abs(x.flatten()-f_open(x.flatten()))**p) + return npj.mean(npj.abs(x.flatten()-f_open(x.flatten()))**p) * 100 def M_diff_identity_close(x,f_close,p=1): ''' ''' - return npj.mean(npj.abs(x.flatten()-f_close(x.flatten()))**p) + return npj.mean(npj.abs(x.flatten()-f_close(x.flatten()))**p) * 100 From b464432ee65d4a003c52d5e99c21c785dd267328 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Fri, 18 Dec 2020 12:40:48 -0500 Subject: [PATCH 04/15] generalize filters --- python/adjoint/filters.py | 369 ++++++++------------------------------ 1 file changed, 70 insertions(+), 299 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index bfd108d2f..83145a0bd 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -3,9 +3,11 @@ """ import numpy as np +import jax from jax import numpy as npj import meep as mp from scipy import special +from .optimization_problem import atleast_3d def _centered(arr, newshape): @@ -22,291 +24,83 @@ def _centered(arr, newshape): myslice = [slice(startind[k], endind[k]) for k in range(len(endind))] return arr[tuple(myslice)] +def meep_filter(x,kernel): + ''' + General convolution for 1D, 2D, and 3D design spaces. -def _edge_pad(arr, pad): - - # fill sides - left = npj.tile(arr[0,:],(pad[0][0],1)) # left side - right = npj.tile(arr[-1,:],(pad[0][1],1)) # right side - top = npj.tile(arr[:,0],(pad[1][0],1)).transpose() # top side - bottom = npj.tile(arr[:,-1],(pad[1][1],1)).transpose() # bottom side) - - # fill corners - top_left = npj.tile(arr[0,0], (pad[0][0],pad[1][0])) # top left - top_right = npj.tile(arr[-1,0], (pad[0][1],pad[1][0])) # top right - bottom_left = npj.tile(arr[0,-1], (pad[0][0],pad[1][1])) # bottom left - bottom_right = npj.tile(arr[-1,-1], (pad[0][1],pad[1][1])) # bottom right - - out = npj.concatenate(( - npj.concatenate((top_left,top,top_right)), - npj.concatenate((left,arr,right)), - npj.concatenate((bottom_left,bottom,bottom_right)) - ),axis=1) - - return out - - -def _zero_pad(arr, pad): - - # fill sides - left = npj.tile(0,(pad[0][0],arr.shape[1])) # left side - right = npj.tile(0,(pad[0][1],arr.shape[1])) # right side - top = npj.tile(0,(arr.shape[0],pad[1][0])) # top side - bottom = npj.tile(0,(arr.shape[0],pad[1][1])) # bottom side - - # fill corners - top_left = npj.tile(0, (pad[0][0],pad[1][0])) # top left - top_right = npj.tile(0, (pad[0][1],pad[1][0])) # top right - bottom_left = npj.tile(0, (pad[0][0],pad[1][1])) # bottom left - bottom_right = npj.tile(0, (pad[0][1],pad[1][1])) # bottom right - - out = npj.concatenate(( - npj.concatenate((top_left,top,top_right)), - npj.concatenate((left,arr,right)), - npj.concatenate((bottom_left,bottom,bottom_right)) - ),axis=1) - - return out - - -def simple_2d_filter(x, kernel, Lx, Ly, resolution, symmetries=[]): - """A simple 2d filter algorithm that is differentiable with autograd. - Uses a 2D fft approach since it is typically faster and preserves the shape - of the input and output arrays. - - The ffts pad the operation to prevent any circular convolution garbage. + x and kernel must have the shape you wish to filter. + ''' + # store shapes + x_shape = x.shape; k_shape = kernel.shape - Parameters - ---------- - x : array_like (2D) - Input array to be filtered. Must be 2D. - kernel : array_like (2D) - Filter kernel (before the DFT). Must be same size as `x` - Lx : float - Length of design region in X direction (in "meep units") - Ly : float - Length of design region in Y direction (in "meep units") - resolution : int - Resolution of the design grid (not the meep simulation resolution) - symmetries : list - Symmetries to impose on the parameter field (either mp.X or mp.Y) + # edge pad + npad = *((s,s) for s in x_shape), + kernelp = npj.pad(kernel,npad,mode='edge') + xp = npj.pad(x,npad,mode='edge') - Returns - ------- - array_like (2D) - The output of the 2d convolution. - """ - # Get 2d parameter space shape - Nx = int(Lx * resolution) - Ny = int(Ly * resolution) - (kx, ky) = kernel.shape - - # Adjust parameter space shape for symmetries - if mp.X in symmetries: - Nx = int(Nx / 2) - if mp.Y in symmetries: - Ny = int(Ny / 2) - - # Ensure the input is 2D - x = x.reshape(Nx, Ny) - - # Perform the required reflections for symmetries - if mp.X in symmetries: - if kx % 2 == 1: - x = npj.concatenate((x,x[-1,:][None,:],x[::-1,:]), axis=0) - else: - x = npj.concatenate((x,x[::-1,:]), axis=0) - if mp.Y in symmetries: - if ky % 2 == 1: - x = npj.concatenate((x[:,::-1],x[:,-1][:,None],x), axis=1) - else: - x = npj.concatenate((x[:,::-1],x), axis=1) - - # pad the kernel and input to avoid circular convolution and - # to ensure boundary conditions are met. - kernel = _zero_pad(kernel, ((kx, kx), (ky, ky))) - x = _edge_pad(x, ((kx, kx), (ky, ky))) - - # Transform to frequency domain for fast convolution - H = npj.fft.fft2(kernel) - X = npj.fft.fft2(x) - - # 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 = npj.fft.fftshift(npj.real(npj.fft.ifft2(Y))) - - # Remove all the extra padding - y = _centered(y, (kx, ky)) - - # Remove the added symmetry domains - if mp.X in symmetries: - y = y[0:Nx, :] - if mp.Y in symmetries: - y = y[:, -Ny:] - - return y - - -def cylindrical_filter(x, radius, Lx, Ly, resolution, symmetries=[]): - '''A uniform cylindrical filter [1]. Typically allows for sharper transitions. + # convolve + yp = jax.scipy.signal.convolve(xp,kernelp,mode='same') - Parameters - ---------- - x : array_like (2D) - Design parameters - radius : float - Filter radius (in "meep units") - Lx : float - Length of design region in X direction (in "meep units") - Ly : float - Length of design region in Y direction (in "meep units") - resolution : int - Resolution of the design grid (not the meep simulation resolution) - symmetries : list - Symmetries to impose on the parameter field (either mp.X or mp.Y) + # remove paddings + return _centered(yp,x_shape) - Returns - ------- - array_like (2D) - Filtered design parameters. - - References - ---------- - [1] Lazarov, B. S., Wang, F., & Sigmund, O. (2016). Length scale and manufacturability in - density-based topology optimization. Archive of Applied Mechanics, 86(1-2), 189-218. +def symmetric_kernel(kernel,dims): ''' - # Get 2d parameter space shape - Nx = int(Lx * resolution) - Ny = int(Ly * resolution) - - # Formulate grid over entire design region - xv, yv = np.meshgrid(np.linspace(-Lx / 2, Lx / 2, Nx), - np.linspace(-Ly / 2, Ly / 2, Ny), - sparse=True, - indexing='ij') - - # Calculate kernel - kernel = np.where(np.abs(xv**2 + yv**2) <= radius**2, 1, 0).T - - # Normalize kernel - kernel = kernel / np.sum(kernel.flatten()) # Normalize the filter - - # Filter the response - y = simple_2d_filter(x, kernel, Lx, Ly, resolution, symmetries) - - return y - - -def conic_filter(x, radius, Lx, Ly, resolution, symmetries=[]): - '''A linear conic filter, also known as a "Hat" filter in the literature [1]. - - Parameters - ---------- - x : array_like (2D) - Design parameters - radius : float - Filter radius (in "meep units") - Lx : float - Length of design region in X direction (in "meep units") - Ly : float - Length of design region in Y direction (in "meep units") - resolution : int - Resolution of the design grid (not the meep simulation resolution) - symmetries : list - Symmetries to impose on the parameter field (either mp.X or mp.Y) - - Returns - ------- - array_like (2D) - Filtered design parameters. - - References - ---------- - [1] Lazarov, B. S., Wang, F., & Sigmund, O. (2016). Length scale and manufacturability in - density-based topology optimization. Archive of Applied Mechanics, 86(1-2), 189-218. + Given a 1D array for a kernel, do sucessive outer products + to get higher dimensional equivalents. ''' - # Get 2d parameter space shape - Nx = int(Lx * resolution) - Ny = int(Ly * resolution) - - # Formulate grid over entire design region - xv, yv = np.meshgrid(np.linspace(-Lx / 2, Lx / 2, Nx), - np.linspace(-Ly / 2, Ly / 2, Ny), - sparse=True, - indexing='ij') - - # Calculate kernel - kernel = np.where( - np.abs(xv**2 + yv**2) <= radius**2, - (1 - np.sqrt(abs(xv**2 + yv**2)) / radius), 0) - - # Normalize kernel - kernel = kernel / np.sum(kernel.flatten()) # Normalize the filter - - # Filter the response - y = simple_2d_filter(x, kernel, Lx, Ly, resolution, symmetries) - - return y + k = kernel.copy() + for k in range(dims-1): + k = npj.unfunc.outer(k,kernel) + + return k +def cylindrical_filter(x,radius): + x = atleast_3d(x) + xl = npj.linspace(-x.shape[0]/2,x.shape[0]/2,x.shape[0]) + yl = npj.linspace(-x.shape[1]/2,x.shape[1]/2,x.shape[1]) + zl = npj.linspace(-x.shape[2]/2,x.shape[2]/2,x.shape[2]) -def gaussian_filter(x, sigma, Lx, Ly, resolution, symmetries=[]): - '''A simple gaussian filter of the form exp(-x **2 / sigma ** 2) [1]. + X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True) + kernel = X**2+Y**2+Z**2 <= radius**2 + kernel = kernel/npj.sum(kernel.flatten()) # normalize - Parameters - ---------- - x : array_like (2D) - Design parameters - sigma : float - Filter radius (in "meep units") - Lx : float - Length of design region in X direction (in "meep units") - Ly : float - Length of design region in Y direction (in "meep units") - resolution : int - Resolution of the design grid (not the meep simulation resolution) + return meep_filter(npj.squeeze(x),npj.squeeze(kernel)) - Returns - ------- - array_like (2D) - Filtered design parameters. - References - ---------- - [1] Wang, E. W., Sell, D., Phan, T., & Fan, J. A. (2019). Robust design of - topology-optimized metasurfaces. Optical Materials Express, 9(2), 469-482. - ''' - # Get 2d parameter space shape - Nx = int(Lx * resolution) - Ny = int(Ly * resolution) - gaussian = lambda x, sigma: np.exp(-x**2 / sigma**2) +def conic_filter(x,radius): + x = atleast_3d(x) + xl = npj.linspace(-x.shape[0]/2,x.shape[0]/2,x.shape[0]) + yl = npj.linspace(-x.shape[1]/2,x.shape[1]/2,x.shape[1]) + zl = npj.linspace(-x.shape[2]/2,x.shape[2]/2,x.shape[2]) - # Formulate grid over entire design region - xv = np.linspace(-Lx / 2, Lx / 2, Nx) - yv = np.linspace(-Ly / 2, Ly / 2, Ny) + X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True) + kernel = np.where(np.abs(X**2+Y**2+Z**2) <= radius**2,(1-np.sqrt(npj.abs(X**2+Y**2+Z**2))/radius),0) + kernel = kernel/npj.sum(kernel.flatten()) # normalize - # Calculate kernel - kernel = np.outer(gaussian(xv, sigma), - gaussian(yv, sigma)) # Gaussian filter kernel + return meep_filter(npj.squeeze(x),npj.squeeze(kernel)) - # Normalize kernel - kernel = kernel / np.sum(kernel.flatten()) # Normalize the filter +def gaussian_filter(x,sigma): + gaussian = lambda x,sigma: np.exp(-x**2/sigma**2) + x = atleast_3d(x) + xl = npj.linspace(-x.shape[0]/2,x.shape[0]/2,x.shape[0]) + yl = npj.linspace(-x.shape[1]/2,x.shape[1]/2,x.shape[1]) + zl = npj.linspace(-x.shape[2]/2,x.shape[2]/2,x.shape[2]) - # Filter the response - y = simple_2d_filter(x, kernel, Lx, Ly, resolution, symmetries) + X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True) + kernel = np.multiply.outer(np.outer(gaussian(X, sigma), gaussian(Y, sigma)),gaussian(Z, sigma)) # Gaussian filter kernel + kernel = kernel/npj.sum(kernel.flatten()) # normalize - return y + return meep_filter(npj.squeeze(x),npj.squeeze(kernel)) ''' # ------------------------------------------------------------------------------------ # Erosion and dilation operators ''' - - -def exponential_erosion(x, radius, beta, Lx, Ly, resolution): + +def exponential_erosion(x,radius,beta): ''' Performs and exponential erosion operation. Parameters @@ -339,14 +133,9 @@ def exponential_erosion(x, radius, beta, Lx, Ly, resolution): ''' x_hat = npj.exp(beta*(1-x)) - return 1 - npj.log(cylindrical_filter(x_hat,radius,Lx,Ly,resolution).flatten()) / beta - - x_hat = npa.exp(beta * (1 - x)) - return 1 - npa.log( - cylindrical_filter(x_hat, radius, Lx, Ly, resolution).flatten()) / beta - + return 1 - npj.log(cylindrical_filter(x_hat,radius)) / beta -def exponential_dilation(x, radius, beta, Lx, Ly, resolution): +def exponential_dilation(x,radius,beta): ''' Performs a exponential dilation operation. Parameters @@ -379,14 +168,9 @@ def exponential_dilation(x, radius, beta, Lx, Ly, resolution): ''' x_hat = npj.exp(beta*x) - return npj.log(cylindrical_filter(x_hat,radius,Lx,Ly,resolution).flatten()) / beta + return npj.log(cylindrical_filter(x_hat,radius)) / beta - x_hat = npa.exp(beta * x) - return npa.log( - cylindrical_filter(x_hat, radius, Lx, Ly, resolution).flatten()) / beta - - -def heaviside_erosion(x, radius, beta, Lx, Ly, resolution): +def heaviside_erosion(x,radius,beta): ''' Performs a heaviside erosion operation. Parameters @@ -416,14 +200,10 @@ def heaviside_erosion(x, radius, beta, Lx, Ly, resolution): numerical methods in engineering, 61(2), 238-254. ''' - x_hat = cylindrical_filter(x,radius,Lx,Ly,resolution).flatten() + x_hat = cylindrical_filter(x,radius) return npj.exp(-beta*(1-x_hat)) + npj.exp(-beta)*(1-x_hat) - x_hat = cylindrical_filter(x, radius, Lx, Ly, resolution).flatten() - return npa.exp(-beta * (1 - x_hat)) + npa.exp(-beta) * (1 - x_hat) - - -def heaviside_dilation(x, radius, beta, Lx, Ly, resolution): +def heaviside_dilation(x,radius,beta): ''' Performs a heaviside dilation operation. Parameters @@ -453,14 +233,10 @@ def heaviside_dilation(x, radius, beta, Lx, Ly, resolution): numerical methods in engineering, 61(2), 238-254. ''' - x_hat = cylindrical_filter(x,radius,Lx,Ly,resolution).flatten() + x_hat = cylindrical_filter(x,radius) return 1 - npj.exp(-beta*x_hat) + npj.exp(-beta)*x_hat - x_hat = cylindrical_filter(x, radius, Lx, Ly, resolution).flatten() - return 1 - npa.exp(-beta * x_hat) + npa.exp(-beta) * x_hat - - -def geometric_erosion(x, radius, alpha, Lx, Ly, resolution): +def geometric_erosion(x,radius,alpha): ''' Performs a geometric erosion operation. Parameters @@ -489,10 +265,9 @@ def geometric_erosion(x, radius, alpha, Lx, Ly, resolution): Pythagorean means. Structural and Multidisciplinary Optimization, 48(5), 859-875. ''' x_hat = npj.log(x + alpha) - return npj.exp(cylindrical_filter(x_hat,radius,Lx,Ly,resolution)).flatten() - alpha + return npj.exp(cylindrical_filter(x_hat,radius)) - alpha - -def geometric_dilation(x, radius, alpha, Lx, Ly, resolution): +def geometric_dilation(x,radius,alpha): ''' Performs a geometric dilation operation. Parameters @@ -522,9 +297,9 @@ def geometric_dilation(x, radius, alpha, Lx, Ly, resolution): ''' x_hat = npj.log(1 - x + alpha) - return -npj.exp(cylindrical_filter(x_hat,radius,Lx,Ly,resolution)).flatten() + alpha + 1 + return -npj.exp(cylindrical_filter(x_hat,radius)) + alpha + 1 -def harmonic_erosion(x, radius, alpha, Lx, Ly, resolution): +def harmonic_erosion(x,radius,alpha): ''' Performs a harmonic erosion operation. Parameters @@ -554,11 +329,9 @@ def harmonic_erosion(x, radius, alpha, Lx, Ly, resolution): ''' x_hat = 1 / (x + alpha) - return 1 / cylindrical_filter(x_hat, radius, Lx, Ly, - resolution).flatten() - alpha - + return 1 / cylindrical_filter(x_hat,radius) - alpha -def harmonic_dilation(x, radius, alpha, Lx, Ly, resolution): +def harmonic_dilation(x,radius,alpha): ''' Performs a harmonic dilation operation. Parameters @@ -588,9 +361,7 @@ def harmonic_dilation(x, radius, alpha, Lx, Ly, resolution): ''' x_hat = 1 / (1 - x + alpha) - return 1 - 1 / cylindrical_filter(x_hat, radius, Lx, Ly, - resolution).flatten() + alpha - + return 1 - 1 / cylindrical_filter(x_hat,radius) + alpha ''' # ------------------------------------------------------------------------------------ # From 6b5f83fd9218d880ae144b2ba07a0576838f9a19 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Sat, 19 Dec 2020 12:06:13 -0500 Subject: [PATCH 05/15] rever --- python/adjoint/filters.py | 32 +++++++++++++++----------------- 1 file changed, 15 insertions(+), 17 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index 83145a0bd..9221a0318 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -7,7 +7,18 @@ from jax import numpy as npj import meep as mp from scipy import special -from .optimization_problem import atleast_3d + + +def atleast_3d(ary): + if ary.ndim == 0: + result = npj.expand_dims(ary, axis=(0,1,2))#ary.reshape(1, 1, 1) + elif ary.ndim == 1: + result = npj.expand_dims(ary, axis=(1,2)) + elif ary.ndim == 2: + result = npj.expand_dims(ary, axis=(2)) + else: + result = ary + return result def _centered(arr, newshape): @@ -17,8 +28,8 @@ def _centered(arr, newshape): https://github.com/scipy/scipy/blob/v1.4.1/scipy/signal/signaltools.py#L263-L270 ''' # Return the center newshape portion of the array. - newshape = np.asarray(newshape) - currshape = np.array(arr.shape) + newshape = npj.array(newshape) + currshape = npj.array(arr.shape) startind = (currshape - newshape) // 2 endind = startind + newshape myslice = [slice(startind[k], endind[k]) for k in range(len(endind))] @@ -44,17 +55,6 @@ def meep_filter(x,kernel): # remove paddings return _centered(yp,x_shape) -def symmetric_kernel(kernel,dims): - ''' - Given a 1D array for a kernel, do sucessive outer products - to get higher dimensional equivalents. - ''' - k = kernel.copy() - for k in range(dims-1): - k = npj.unfunc.outer(k,kernel) - - return k - def cylindrical_filter(x,radius): x = atleast_3d(x) xl = npj.linspace(-x.shape[0]/2,x.shape[0]/2,x.shape[0]) @@ -64,10 +64,8 @@ def cylindrical_filter(x,radius): X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True) kernel = X**2+Y**2+Z**2 <= radius**2 kernel = kernel/npj.sum(kernel.flatten()) # normalize - - return meep_filter(npj.squeeze(x),npj.squeeze(kernel)) - + return meep_filter(npj.squeeze(x),npj.squeeze(kernel)) def conic_filter(x,radius): x = atleast_3d(x) From 8480e49260a10e974cc3ff064472a668d886793b Mon Sep 17 00:00:00 2001 From: smartalecH Date: Tue, 29 Dec 2020 13:59:44 -0700 Subject: [PATCH 06/15] fix convolution --- python/adjoint/filters.py | 30 +++++++++++++++++++++++++++--- 1 file changed, 27 insertions(+), 3 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index 9221a0318..058d203f4 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -7,6 +7,7 @@ from jax import numpy as npj import meep as mp from scipy import special +import jax.numpy.fft as fft def atleast_3d(ary): @@ -39,18 +40,41 @@ def meep_filter(x,kernel): ''' General convolution for 1D, 2D, and 3D design spaces. - x and kernel must have the shape you wish to filter. + Currently, Jax's XLA convolution kernel *always* uses + the overlap-add method, which is very costly for + the large arrays we typically convolve. + + In the meantime, we can do a manual n-dimensional + FFT (also using jax). Note, however, that Jax + doesn't support automatic padding (like numpy does) + so we need to do the padding ourselves (in addition + to the edge padding we use for boundary conditions). + + Note that we need to fftshift the output since we + padded *both* sides of our signals. + + Parameters + ---------- + x : array_like (1D, 2D, or 3D) + Input array to be filtered. + kernel : array_like (1D, 2D, or 3D) + Filter kernel (before the DFT). Must be same size as `x` + + Returns + ------- + array_like (1D, 2D, or 3D) + The output of the N-D convolution. ''' # store shapes x_shape = x.shape; k_shape = kernel.shape - + # edge pad npad = *((s,s) for s in x_shape), kernelp = npj.pad(kernel,npad,mode='edge') xp = npj.pad(x,npad,mode='edge') # convolve - yp = jax.scipy.signal.convolve(xp,kernelp,mode='same') + yp = fft.fftshift(fft.ifftn(fft.fftn(xp) * fft.fftn(kernelp))).real # remove paddings return _centered(yp,x_shape) From a0766013eb10835cde5bd353254a652c219aedbc Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Fri, 8 Jan 2021 16:04:36 -0500 Subject: [PATCH 07/15] fix filters --- python/adjoint/filters.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index 058d203f4..5b08b68e6 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -9,6 +9,18 @@ from scipy import special import jax.numpy.fft as fft +def centered(arr, newshape): + '''Helper function that reformats the padded array of the fft filter operation. + Borrowed from scipy: + https://github.com/scipy/scipy/blob/v1.4.1/scipy/signal/signaltools.py#L263-L270 + ''' + # Return the center newshape portion of the array. + newshape = np.asarray(newshape) + currshape = np.array(arr.shape) + startind = (currshape - newshape) // 2 + endind = startind + newshape + myslice = [slice(startind[k], endind[k]) for k in range(len(endind))] + return arr[tuple(myslice)] def atleast_3d(ary): if ary.ndim == 0: @@ -80,6 +92,10 @@ def meep_filter(x,kernel): return _centered(yp,x_shape) def cylindrical_filter(x,radius): + # filter radius is in pixels, so make sure radius >= 1 + if radius < 1: + raise ValueError("The filter radius must be at least 1 pixel") + x = atleast_3d(x) xl = npj.linspace(-x.shape[0]/2,x.shape[0]/2,x.shape[0]) yl = npj.linspace(-x.shape[1]/2,x.shape[1]/2,x.shape[1]) From 1220fc49876a9d78e98d7ba038539ab73652aedd Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Mon, 11 Jan 2021 10:10:14 -0500 Subject: [PATCH 08/15] make filters zero-phase --- python/adjoint/filters.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index 5b08b68e6..1e784fbbb 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -86,7 +86,18 @@ def meep_filter(x,kernel): xp = npj.pad(x,npad,mode='edge') # convolve - yp = fft.fftshift(fft.ifftn(fft.fftn(xp) * fft.fftn(kernelp))).real + ''' + As is the case with any convolution, the filter will + introduce a delay in every specified dimension. To + compensate for this delay, we can perform a zero-phase + filter operation. We first convolve normally, and then + "flip" the signal and convolve again, effectively eliminating + the phase delay. This effectively applies the filter kernel + twice, however. We sqrt the kernel in the frequency domain + to somewhat compensate for this. + ''' + yp = fft.fftshift(fft.ifftn(fft.fftn(xp) * npj.sqrt(fft.fftn(kernelp)))).real + yp = fft.fftshift(fft.ifftn(fft.fftn(npj.flip(yp,axis=None)) * npj.sqrt(fft.fftn(kernelp)))).real # remove paddings return _centered(yp,x_shape) @@ -101,7 +112,7 @@ def cylindrical_filter(x,radius): yl = npj.linspace(-x.shape[1]/2,x.shape[1]/2,x.shape[1]) zl = npj.linspace(-x.shape[2]/2,x.shape[2]/2,x.shape[2]) - X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True) + X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True,indexing='ij') kernel = X**2+Y**2+Z**2 <= radius**2 kernel = kernel/npj.sum(kernel.flatten()) # normalize @@ -113,7 +124,7 @@ def conic_filter(x,radius): yl = npj.linspace(-x.shape[1]/2,x.shape[1]/2,x.shape[1]) zl = npj.linspace(-x.shape[2]/2,x.shape[2]/2,x.shape[2]) - X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True) + X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True,indexing='ij') kernel = np.where(np.abs(X**2+Y**2+Z**2) <= radius**2,(1-np.sqrt(npj.abs(X**2+Y**2+Z**2))/radius),0) kernel = kernel/npj.sum(kernel.flatten()) # normalize @@ -126,7 +137,7 @@ def gaussian_filter(x,sigma): yl = npj.linspace(-x.shape[1]/2,x.shape[1]/2,x.shape[1]) zl = npj.linspace(-x.shape[2]/2,x.shape[2]/2,x.shape[2]) - X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True) + X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True,indexing='ij') kernel = np.multiply.outer(np.outer(gaussian(X, sigma), gaussian(Y, sigma)),gaussian(Z, sigma)) # Gaussian filter kernel kernel = kernel/npj.sum(kernel.flatten()) # normalize From 41593691f90d371fd06dc4d4070da4bf1ccb5354 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Tue, 12 Jan 2021 17:12:42 -0500 Subject: [PATCH 09/15] fix small filter bugs --- python/adjoint/filters.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index 1e784fbbb..8e95abbfc 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -84,7 +84,7 @@ def meep_filter(x,kernel): npad = *((s,s) for s in x_shape), kernelp = npj.pad(kernel,npad,mode='edge') xp = npj.pad(x,npad,mode='edge') - + # convolve ''' As is the case with any convolution, the filter will @@ -96,8 +96,9 @@ def meep_filter(x,kernel): twice, however. We sqrt the kernel in the frequency domain to somewhat compensate for this. ''' - yp = fft.fftshift(fft.ifftn(fft.fftn(xp) * npj.sqrt(fft.fftn(kernelp)))).real - yp = fft.fftshift(fft.ifftn(fft.fftn(npj.flip(yp,axis=None)) * npj.sqrt(fft.fftn(kernelp)))).real + + yp = fft.fftshift(fft.ifftn(fft.fftn(xp) * fft.fftn(kernelp))).real + yp = fft.fftshift(fft.ifftn(fft.fftn(npj.flip(yp,axis=None)) * fft.fftn(kernelp))).real # remove paddings return _centered(yp,x_shape) From aedc294711c68a76fb1fa36a57781e8ea901daee Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Wed, 13 Jan 2021 14:21:37 -0500 Subject: [PATCH 10/15] fix flipped filters --- python/adjoint/filters.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index 8e95abbfc..dbeb72b3b 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -98,7 +98,7 @@ def meep_filter(x,kernel): ''' yp = fft.fftshift(fft.ifftn(fft.fftn(xp) * fft.fftn(kernelp))).real - yp = fft.fftshift(fft.ifftn(fft.fftn(npj.flip(yp,axis=None)) * fft.fftn(kernelp))).real + yp = npj.flip(fft.fftshift(fft.ifftn(fft.fftn(npj.flip(yp,axis=None)) * fft.fftn(kernelp))).real,axis=None) # remove paddings return _centered(yp,x_shape) @@ -751,29 +751,29 @@ def gray_indicator(x): [1] Lazarov, B. S., Wang, F., & Sigmund, O. (2016). Length scale and manufacturability in density-based topology optimization. Archive of Applied Mechanics, 86(1-2), 189-218. ''' - return npj.mean(4 * x.flatten() * (1-x.flatten())) * 100 + return npj.mean(4 * x.flatten() * (1-x.flatten())) -def F_diff_open_close(x,f_open,f_close,beta=64): +def F_diff_open_close(x,f_open,f_close,beta=64,p=2): ''' ''' - return npj.mean(tanh_projection(npj.abs(f_close(x.flatten())-f_open(x.flatten()))),beta,0.5) + return npj.mean(tanh_projection(npj.abs(f_close(x)-f_open(x)).flatten()**2),beta,0.5) * 100 -def M_diff_open_close(x,f_open,f_close,p=1): +def M_diff_open_close(x,f_open,f_close,p=2): ''' ''' - return npj.mean(npj.abs(f_close(x.flatten())-f_open(x.flatten()))**p) * 100 + return npj.mean(npj.abs(f_close(x)-f_open(x)).flatten()**p) * 100 -def M_diff_identity_open(x,f_open,p=1): +def M_diff_identity_open(x,f_open,p=2): ''' ''' - return npj.mean(npj.abs(x.flatten()-f_open(x.flatten()))**p) * 100 + return npj.mean(npj.abs(x-f_open(x)).flatten()**p) * 100 -def M_diff_identity_close(x,f_close,p=1): +def M_diff_identity_close(x,f_close,p=2): ''' ''' - return npj.mean(npj.abs(x.flatten()-f_close(x.flatten()))**p) * 100 + return npj.mean(npj.abs(x-f_close(x)).flatten()**p) * 100 From 76e1c5f83632a39ef08d4b0fe5c7522e887a6803 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Wed, 20 Jan 2021 10:09:56 -0500 Subject: [PATCH 11/15] update geometric constraints --- python/adjoint/filters.py | 35 +++++++++++++++-------------------- 1 file changed, 15 insertions(+), 20 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index dbeb72b3b..53ccfb05c 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -614,15 +614,12 @@ def indicator_solid(x, c, filter_f, threshold_f, resolution): filtered_field = filter_f(x) design_field = threshold_f(filtered_field) - gradient_filtered_field = npj.gradient(filtered_field) - grad_mag = (gradient_filtered_field[0]*resolution) ** 2 + (gradient_filtered_field[1]*resolution) ** 2 - if grad_mag.ndim != 2: - raise ValueError("The gradient fields must be 2 dimensional. Check input array and filter functions.") - I_s = design_field * npj.exp(-c * grad_mag) - return I_s + gradient_filtered_field = npj.gradient(filtered_field,resolution) + #gradient_filtered_field = gradient_filtered_field if isinstance(gradient_filtered_field,list) else [gradient_filtered_field] + grad_mag = npj.sum(npj.array(gradient_filtered_field)**2,axis=0) + return design_field * npj.exp(-c * grad_mag) - -def constraint_solid(x, c, eta_e, filter_f, threshold_f, resolution): +def constraint_solid(x,c,eta_e,filter_f,threshold_f,resolution=1): '''Calculates the constraint function of the solid phase needed for minimum length optimization [1]. Parameters @@ -655,8 +652,8 @@ def constraint_solid(x, c, eta_e, filter_f, threshold_f, resolution): ''' filtered_field = filter_f(x) - I_s = indicator_solid(x.reshape(filtered_field.shape),c,filter_f,threshold_f,resolution).flatten() - return npj.mean(I_s * npj.minimum(filtered_field.flatten()-eta_e,0)**2) + I_s = indicator_solid(x,c,filter_f,threshold_f,resolution) + return npj.mean(I_s.flatten() * npj.minimum(filtered_field.flatten()-eta_e,0)**2) def indicator_void(x, c, filter_f, threshold_f, resolution): @@ -685,17 +682,15 @@ def indicator_void(x, c, filter_f, threshold_f, resolution): [1] Zhou, M., Lazarov, B. S., Wang, F., & Sigmund, O. (2015). Minimum length scale in topology optimization by geometric constraints. Computer Methods in Applied Mechanics and Engineering, 293, 266-282. ''' - - filtered_field = filter_f(x).reshape(x.shape) + + filtered_field = filter_f(x) design_field = threshold_f(filtered_field) - gradient_filtered_field = npj.gradient(filtered_field) - grad_mag = (gradient_filtered_field[0]*resolution) ** 2 + (gradient_filtered_field[1]*resolution) ** 2 - if grad_mag.ndim != 2: - raise ValueError("The gradient fields must be 2 dimensional. Check input array and filter functions.") + gradient_filtered_field = npj.gradient(filtered_field,resolution) + #gradient_filtered_field = gradient_filtered_field if isinstance(list, gradient_filtered_field) else [gradient_filtered_field] + grad_mag = npj.sum(npj.array(gradient_filtered_field)**2,axis=0) return (1 - design_field) * npj.exp(-c * grad_mag) - -def constraint_void(x, c, eta_d, filter_f, threshold_f, resolution): +def constraint_void(x,c,eta_d,filter_f,threshold_f,resolution=1): '''Calculates the constraint function of the void phase needed for minimum length optimization [1]. Parameters @@ -728,8 +723,8 @@ def constraint_void(x, c, eta_d, filter_f, threshold_f, resolution): ''' filtered_field = filter_f(x) - I_v = indicator_void(x.reshape(filtered_field.shape),c,filter_f,threshold_f,resolution).flatten() - return npj.mean(I_v * npj.minimum(eta_d-filtered_field.flatten(),0)**2) + I_v = indicator_void(x,c,filter_f,threshold_f,resolution) + return npj.mean(I_v.flatten() * npj.minimum(eta_d-filtered_field.flatten(),0)**2) def gray_indicator(x): '''Calculates a measure of "grayness" according to [1]. From d83e2709412866a85971450d9ce290ad952b5bc8 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Wed, 24 Feb 2021 14:58:18 -0500 Subject: [PATCH 12/15] revert filters --- python/adjoint/filters.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index 53ccfb05c..8d81ba71b 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -97,7 +97,7 @@ def meep_filter(x,kernel): to somewhat compensate for this. ''' - yp = fft.fftshift(fft.ifftn(fft.fftn(xp) * fft.fftn(kernelp))).real + yp = fft.fftshift(fft.ifftn(fft.fftn(xp) * (fft.fftn(kernelp)))).real yp = npj.flip(fft.fftshift(fft.ifftn(fft.fftn(npj.flip(yp,axis=None)) * fft.fftn(kernelp))).real,axis=None) # remove paddings From 3ecbec1037f2f636251d48b2cead8c399ce80bb2 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Wed, 16 Jun 2021 12:57:19 -0400 Subject: [PATCH 13/15] reformat --- python/adjoint/basis.py | 25 ++- python/adjoint/filter_source.py | 14 +- python/adjoint/filters.py | 217 +++++++++++++++---------- python/adjoint/objective.py | 7 +- python/adjoint/optimization_problem.py | 176 ++++++++------------ python/adjoint/utils.py | 52 +++--- python/adjoint/wrapper.py | 146 ++++++++--------- 7 files changed, 308 insertions(+), 329 deletions(-) diff --git a/python/adjoint/basis.py b/python/adjoint/basis.py index 3e138ee42..25370f272 100644 --- a/python/adjoint/basis.py +++ b/python/adjoint/basis.py @@ -14,13 +14,11 @@ class Basis(ABC): """ """ - def __init__( - self, - rho_vector=None, - volume=None, - size=None, - center=mp.Vector3(), - ): + def __init__(self, + rho_vector=None, + volume=None, + size=None, + center=mp.Vector3()): self.volume = volume if volume else mp.Volume(center=center, size=size) self.rho_vector = rho_vector @@ -54,6 +52,9 @@ class BilinearInterpolationBasis(Basis): Simple bilinear interpolation basis set. ''' def __init__(self, resolution, symmetry=None, **kwargs): + ''' + + ''' self.dim = 2 super(BilinearInterpolationBasis, self).__init__(**kwargs) @@ -194,14 +195,8 @@ def get_bilinear_row(self, rx, ry, rho_x, rho_y): return weights, interp_idx - def gen_interpolation_matrix( - self, - rho_x, - rho_y, - rho_x_interp, - rho_y_interp, - rho_z_interp, - ): + def gen_interpolation_matrix(self, rho_x, rho_y, rho_x_interp, + rho_y_interp, rho_z_interp): ''' Generates a bilinear interpolation matrix. diff --git a/python/adjoint/filter_source.py b/python/adjoint/filter_source.py index cce5fc881..b69ef2d80 100644 --- a/python/adjoint/filter_source.py +++ b/python/adjoint/filter_source.py @@ -4,14 +4,12 @@ class FilteredSource(CustomSource): - def __init__( - self, - center_frequency, - frequencies, - frequency_response, - dt, - time_src=None, - ): + def __init__(self, + center_frequency, + frequencies, + frequency_response, + dt, + time_src=None): dt = dt / 2 # divide by two to compensate for staggered E,H time interval self.dt = dt self.frequencies = frequencies diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index 8d81ba71b..833cf2df5 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -9,6 +9,7 @@ from scipy import special import jax.numpy.fft as fft + def centered(arr, newshape): '''Helper function that reformats the padded array of the fft filter operation. Borrowed from scipy: @@ -22,11 +23,12 @@ def centered(arr, newshape): myslice = [slice(startind[k], endind[k]) for k in range(len(endind))] return arr[tuple(myslice)] + def atleast_3d(ary): if ary.ndim == 0: - result = npj.expand_dims(ary, axis=(0,1,2))#ary.reshape(1, 1, 1) + result = npj.expand_dims(ary, axis=(0, 1, 2)) #ary.reshape(1, 1, 1) elif ary.ndim == 1: - result = npj.expand_dims(ary, axis=(1,2)) + result = npj.expand_dims(ary, axis=(1, 2)) elif ary.ndim == 2: result = npj.expand_dims(ary, axis=(2)) else: @@ -48,7 +50,8 @@ def _centered(arr, newshape): myslice = [slice(startind[k], endind[k]) for k in range(len(endind))] return arr[tuple(myslice)] -def meep_filter(x,kernel): + +def meep_filter(x, kernel): ''' General convolution for 1D, 2D, and 3D design spaces. @@ -78,12 +81,13 @@ def meep_filter(x,kernel): The output of the N-D convolution. ''' # store shapes - x_shape = x.shape; k_shape = kernel.shape - + x_shape = x.shape + k_shape = kernel.shape + # edge pad - npad = *((s,s) for s in x_shape), - kernelp = npj.pad(kernel,npad,mode='edge') - xp = npj.pad(x,npad,mode='edge') + npad = *((s, s) for s in x_shape), + kernelp = npj.pad(kernel, npad, mode='edge') + xp = npj.pad(x, npad, mode='edge') # convolve ''' @@ -98,59 +102,70 @@ def meep_filter(x,kernel): ''' yp = fft.fftshift(fft.ifftn(fft.fftn(xp) * (fft.fftn(kernelp)))).real - yp = npj.flip(fft.fftshift(fft.ifftn(fft.fftn(npj.flip(yp,axis=None)) * fft.fftn(kernelp))).real,axis=None) - + yp = npj.flip(fft.fftshift( + fft.ifftn(fft.fftn(npj.flip(yp, axis=None)) * fft.fftn(kernelp))).real, + axis=None) + # remove paddings - return _centered(yp,x_shape) + return _centered(yp, x_shape) + -def cylindrical_filter(x,radius): +def cylindrical_filter(x, radius): # filter radius is in pixels, so make sure radius >= 1 if radius < 1: - raise ValueError("The filter radius must be at least 1 pixel") + raise ValueError("The filter radius must be at least 1 pixel") x = atleast_3d(x) - xl = npj.linspace(-x.shape[0]/2,x.shape[0]/2,x.shape[0]) - yl = npj.linspace(-x.shape[1]/2,x.shape[1]/2,x.shape[1]) - zl = npj.linspace(-x.shape[2]/2,x.shape[2]/2,x.shape[2]) + xl = npj.linspace(-x.shape[0] / 2, x.shape[0] / 2, x.shape[0]) + yl = npj.linspace(-x.shape[1] / 2, x.shape[1] / 2, x.shape[1]) + zl = npj.linspace(-x.shape[2] / 2, x.shape[2] / 2, x.shape[2]) + + X, Y, Z = npj.meshgrid(xl, yl, zl, sparse=True, indexing='ij') + kernel = X**2 + Y**2 + Z**2 <= radius**2 + kernel = kernel / npj.sum(kernel.flatten()) # normalize - X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True,indexing='ij') - kernel = X**2+Y**2+Z**2 <= radius**2 - kernel = kernel/npj.sum(kernel.flatten()) # normalize + return meep_filter(npj.squeeze(x), npj.squeeze(kernel)) - return meep_filter(npj.squeeze(x),npj.squeeze(kernel)) -def conic_filter(x,radius): +def conic_filter(x, radius): x = atleast_3d(x) - xl = npj.linspace(-x.shape[0]/2,x.shape[0]/2,x.shape[0]) - yl = npj.linspace(-x.shape[1]/2,x.shape[1]/2,x.shape[1]) - zl = npj.linspace(-x.shape[2]/2,x.shape[2]/2,x.shape[2]) + xl = npj.linspace(-x.shape[0] / 2, x.shape[0] / 2, x.shape[0]) + yl = npj.linspace(-x.shape[1] / 2, x.shape[1] / 2, x.shape[1]) + zl = npj.linspace(-x.shape[2] / 2, x.shape[2] / 2, x.shape[2]) - X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True,indexing='ij') - kernel = np.where(np.abs(X**2+Y**2+Z**2) <= radius**2,(1-np.sqrt(npj.abs(X**2+Y**2+Z**2))/radius),0) - kernel = kernel/npj.sum(kernel.flatten()) # normalize + X, Y, Z = npj.meshgrid(xl, yl, zl, sparse=True, indexing='ij') + kernel = np.where( + np.abs(X**2 + Y**2 + Z**2) <= radius**2, + (1 - np.sqrt(npj.abs(X**2 + Y**2 + Z**2)) / radius), 0) + kernel = kernel / npj.sum(kernel.flatten()) # normalize - return meep_filter(npj.squeeze(x),npj.squeeze(kernel)) + return meep_filter(npj.squeeze(x), npj.squeeze(kernel)) -def gaussian_filter(x,sigma): - gaussian = lambda x,sigma: np.exp(-x**2/sigma**2) + +def gaussian_filter(x, sigma): + gaussian = lambda x, sigma: np.exp(-x**2 / sigma**2) x = atleast_3d(x) - xl = npj.linspace(-x.shape[0]/2,x.shape[0]/2,x.shape[0]) - yl = npj.linspace(-x.shape[1]/2,x.shape[1]/2,x.shape[1]) - zl = npj.linspace(-x.shape[2]/2,x.shape[2]/2,x.shape[2]) + xl = npj.linspace(-x.shape[0] / 2, x.shape[0] / 2, x.shape[0]) + yl = npj.linspace(-x.shape[1] / 2, x.shape[1] / 2, x.shape[1]) + zl = npj.linspace(-x.shape[2] / 2, x.shape[2] / 2, x.shape[2]) + + X, Y, Z = npj.meshgrid(xl, yl, zl, sparse=True, indexing='ij') + kernel = np.multiply.outer(np.outer(gaussian(X, sigma), gaussian(Y, + sigma)), + gaussian(Z, sigma)) # Gaussian filter kernel + kernel = kernel / npj.sum(kernel.flatten()) # normalize - X,Y,Z = npj.meshgrid(xl,yl,zl,sparse=True,indexing='ij') - kernel = np.multiply.outer(np.outer(gaussian(X, sigma), gaussian(Y, sigma)),gaussian(Z, sigma)) # Gaussian filter kernel - kernel = kernel/npj.sum(kernel.flatten()) # normalize + return meep_filter(npj.squeeze(x), npj.squeeze(kernel)) - return meep_filter(npj.squeeze(x),npj.squeeze(kernel)) ''' # ------------------------------------------------------------------------------------ # Erosion and dilation operators ''' - -def exponential_erosion(x,radius,beta): + + +def exponential_erosion(x, radius, beta): ''' Performs and exponential erosion operation. Parameters @@ -181,11 +196,12 @@ def exponential_erosion(x,radius,beta): morphological close-open and open-close filters for topology optimization. Structural and Multidisciplinary Optimization, 54(1), 15-21. ''' - - x_hat = npj.exp(beta*(1-x)) - return 1 - npj.log(cylindrical_filter(x_hat,radius)) / beta -def exponential_dilation(x,radius,beta): + x_hat = npj.exp(beta * (1 - x)) + return 1 - npj.log(cylindrical_filter(x_hat, radius)) / beta + + +def exponential_dilation(x, radius, beta): ''' Performs a exponential dilation operation. Parameters @@ -216,11 +232,12 @@ def exponential_dilation(x,radius,beta): morphological close-open and open-close filters for topology optimization. Structural and Multidisciplinary Optimization, 54(1), 15-21. ''' - - x_hat = npj.exp(beta*x) - return npj.log(cylindrical_filter(x_hat,radius)) / beta -def heaviside_erosion(x,radius,beta): + x_hat = npj.exp(beta * x) + return npj.log(cylindrical_filter(x_hat, radius)) / beta + + +def heaviside_erosion(x, radius, beta): ''' Performs a heaviside erosion operation. Parameters @@ -249,11 +266,12 @@ def heaviside_erosion(x,radius,beta): optimization using nodal design variables and projection functions. International journal for numerical methods in engineering, 61(2), 238-254. ''' - - x_hat = cylindrical_filter(x,radius) - return npj.exp(-beta*(1-x_hat)) + npj.exp(-beta)*(1-x_hat) -def heaviside_dilation(x,radius,beta): + x_hat = cylindrical_filter(x, radius) + return npj.exp(-beta * (1 - x_hat)) + npj.exp(-beta) * (1 - x_hat) + + +def heaviside_dilation(x, radius, beta): ''' Performs a heaviside dilation operation. Parameters @@ -282,11 +300,12 @@ def heaviside_dilation(x,radius,beta): optimization using nodal design variables and projection functions. International journal for numerical methods in engineering, 61(2), 238-254. ''' - - x_hat = cylindrical_filter(x,radius) - return 1 - npj.exp(-beta*x_hat) + npj.exp(-beta)*x_hat -def geometric_erosion(x,radius,alpha): + x_hat = cylindrical_filter(x, radius) + return 1 - npj.exp(-beta * x_hat) + npj.exp(-beta) * x_hat + + +def geometric_erosion(x, radius, alpha): ''' Performs a geometric erosion operation. Parameters @@ -315,9 +334,10 @@ def geometric_erosion(x,radius,alpha): Pythagorean means. Structural and Multidisciplinary Optimization, 48(5), 859-875. ''' x_hat = npj.log(x + alpha) - return npj.exp(cylindrical_filter(x_hat,radius)) - alpha + return npj.exp(cylindrical_filter(x_hat, radius)) - alpha -def geometric_dilation(x,radius,alpha): + +def geometric_dilation(x, radius, alpha): ''' Performs a geometric dilation operation. Parameters @@ -347,9 +367,10 @@ def geometric_dilation(x,radius,alpha): ''' x_hat = npj.log(1 - x + alpha) - return -npj.exp(cylindrical_filter(x_hat,radius)) + alpha + 1 + return -npj.exp(cylindrical_filter(x_hat, radius)) + alpha + 1 + -def harmonic_erosion(x,radius,alpha): +def harmonic_erosion(x, radius, alpha): ''' Performs a harmonic erosion operation. Parameters @@ -379,9 +400,10 @@ def harmonic_erosion(x,radius,alpha): ''' x_hat = 1 / (x + alpha) - return 1 / cylindrical_filter(x_hat,radius) - alpha + return 1 / cylindrical_filter(x_hat, radius) - alpha + -def harmonic_dilation(x,radius,alpha): +def harmonic_dilation(x, radius, alpha): ''' Performs a harmonic dilation operation. Parameters @@ -411,7 +433,8 @@ def harmonic_dilation(x,radius,alpha): ''' x_hat = 1 / (1 - x + alpha) - return 1 - 1 / cylindrical_filter(x_hat,radius) + alpha + return 1 - 1 / cylindrical_filter(x_hat, radius) + alpha + ''' # ------------------------------------------------------------------------------------ # @@ -441,8 +464,12 @@ def tanh_projection(x, beta, eta): [1] Wang, F., Lazarov, B. S., & Sigmund, O. (2011). On projection methods, convergence and robust formulations in topology optimization. Structural and Multidisciplinary Optimization, 43(6), 767-784. ''' - - return (npj.tanh(beta*eta) + npj.tanh(beta*(x-eta))) / (npj.tanh(beta*eta) + npj.tanh(beta*(1-eta))) + + return (npj.tanh(beta * eta) + + npj.tanh(beta * + (x - eta))) / (npj.tanh(beta * eta) + npj.tanh(beta * + (1 - eta))) + def heaviside_projection(x, beta, eta): '''Projection filter that thresholds the input parameters between 0 and 1. @@ -466,10 +493,12 @@ def heaviside_projection(x, beta, eta): [1] Lazarov, B. S., Wang, F., & Sigmund, O. (2016). Length scale and manufacturability in density-based topology optimization. Archive of Applied Mechanics, 86(1-2), 189-218. ''' - - case1 = eta*npj.exp(-beta*(eta-x)/eta) - (eta-x)*npj.exp(-beta) - case2 = 1 - (1-eta)*npj.exp(-beta*(x-eta)/(1-eta)) - (eta-x)*npj.exp(-beta) - return npj.where(x < eta,case1,case2) + + case1 = eta * npj.exp(-beta * (eta - x) / eta) - (eta - x) * npj.exp(-beta) + case2 = 1 - (1 - eta) * npj.exp(-beta * (x - eta) / + (1 - eta)) - (eta - x) * npj.exp(-beta) + return npj.where(x < eta, case1, case2) + ''' # ------------------------------------------------------------------------------------ # @@ -614,12 +643,13 @@ def indicator_solid(x, c, filter_f, threshold_f, resolution): filtered_field = filter_f(x) design_field = threshold_f(filtered_field) - gradient_filtered_field = npj.gradient(filtered_field,resolution) + gradient_filtered_field = npj.gradient(filtered_field, resolution) #gradient_filtered_field = gradient_filtered_field if isinstance(gradient_filtered_field,list) else [gradient_filtered_field] - grad_mag = npj.sum(npj.array(gradient_filtered_field)**2,axis=0) + grad_mag = npj.sum(npj.array(gradient_filtered_field)**2, axis=0) return design_field * npj.exp(-c * grad_mag) -def constraint_solid(x,c,eta_e,filter_f,threshold_f,resolution=1): + +def constraint_solid(x, c, eta_e, filter_f, threshold_f, resolution=1): '''Calculates the constraint function of the solid phase needed for minimum length optimization [1]. Parameters @@ -652,8 +682,9 @@ def constraint_solid(x,c,eta_e,filter_f,threshold_f,resolution=1): ''' filtered_field = filter_f(x) - I_s = indicator_solid(x,c,filter_f,threshold_f,resolution) - return npj.mean(I_s.flatten() * npj.minimum(filtered_field.flatten()-eta_e,0)**2) + I_s = indicator_solid(x, c, filter_f, threshold_f, resolution) + return npj.mean(I_s.flatten() * + npj.minimum(filtered_field.flatten() - eta_e, 0)**2) def indicator_void(x, c, filter_f, threshold_f, resolution): @@ -682,15 +713,16 @@ def indicator_void(x, c, filter_f, threshold_f, resolution): [1] Zhou, M., Lazarov, B. S., Wang, F., & Sigmund, O. (2015). Minimum length scale in topology optimization by geometric constraints. Computer Methods in Applied Mechanics and Engineering, 293, 266-282. ''' - + filtered_field = filter_f(x) design_field = threshold_f(filtered_field) - gradient_filtered_field = npj.gradient(filtered_field,resolution) + gradient_filtered_field = npj.gradient(filtered_field, resolution) #gradient_filtered_field = gradient_filtered_field if isinstance(list, gradient_filtered_field) else [gradient_filtered_field] - grad_mag = npj.sum(npj.array(gradient_filtered_field)**2,axis=0) + grad_mag = npj.sum(npj.array(gradient_filtered_field)**2, axis=0) return (1 - design_field) * npj.exp(-c * grad_mag) -def constraint_void(x,c,eta_d,filter_f,threshold_f,resolution=1): + +def constraint_void(x, c, eta_d, filter_f, threshold_f, resolution=1): '''Calculates the constraint function of the void phase needed for minimum length optimization [1]. Parameters @@ -723,8 +755,10 @@ def constraint_void(x,c,eta_d,filter_f,threshold_f,resolution=1): ''' filtered_field = filter_f(x) - I_v = indicator_void(x,c,filter_f,threshold_f,resolution) - return npj.mean(I_v.flatten() * npj.minimum(eta_d-filtered_field.flatten(),0)**2) + I_v = indicator_void(x, c, filter_f, threshold_f, resolution) + return npj.mean(I_v.flatten() * + npj.minimum(eta_d - filtered_field.flatten(), 0)**2) + def gray_indicator(x): '''Calculates a measure of "grayness" according to [1]. @@ -746,29 +780,34 @@ def gray_indicator(x): [1] Lazarov, B. S., Wang, F., & Sigmund, O. (2016). Length scale and manufacturability in density-based topology optimization. Archive of Applied Mechanics, 86(1-2), 189-218. ''' - return npj.mean(4 * x.flatten() * (1-x.flatten())) + return npj.mean(4 * x.flatten() * (1 - x.flatten())) -def F_diff_open_close(x,f_open,f_close,beta=64,p=2): - ''' +def F_diff_open_close(x, f_open, f_close, beta=64, p=2): ''' - return npj.mean(tanh_projection(npj.abs(f_close(x)-f_open(x)).flatten()**2),beta,0.5) * 100 -def M_diff_open_close(x,f_open,f_close,p=2): ''' + return npj.mean( + tanh_projection(npj.abs(f_close(x) - f_open(x)).flatten()**2), beta, + 0.5) * 100 + +def M_diff_open_close(x, f_open, f_close, p=2): ''' - return npj.mean(npj.abs(f_close(x)-f_open(x)).flatten()**p) * 100 -def M_diff_identity_open(x,f_open,p=2): ''' + return npj.mean(npj.abs(f_close(x) - f_open(x)).flatten()**p) * 100 + +def M_diff_identity_open(x, f_open, p=2): ''' - return npj.mean(npj.abs(x-f_open(x)).flatten()**p) * 100 -def M_diff_identity_close(x,f_close,p=2): ''' + return npj.mean(npj.abs(x - f_open(x)).flatten()**p) * 100 + +def M_diff_identity_close(x, f_close, p=2): ''' - return npj.mean(npj.abs(x-f_close(x)).flatten()**p) * 100 + ''' + return npj.mean(npj.abs(x - f_close(x)).flatten()**p) * 100 diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index ecbf2cedd..804f783ad 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -320,11 +320,8 @@ def place_adjoint_source(self, dJ): dJ = dJ.flatten() farpt_list = np.array([list(pi) for pi in self.far_pts]).flatten() far_pt0 = self.far_pts[0] - far_pt_vec = py_v3_to_vec( - self.sim.dimensions, - far_pt0, - self.sim.is_cylindrical, - ) + far_pt_vec = py_v3_to_vec(self.sim.dimensions, far_pt0, + self.sim.is_cylindrical) all_nearsrcdata = self._monitor.swigobj.near_sourcedata( far_pt_vec, farpt_list, self._nfar_pts, dJ) diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py index 80485a20d..14a80119a 100644 --- a/python/adjoint/optimization_problem.py +++ b/python/adjoint/optimization_problem.py @@ -8,14 +8,12 @@ class DesignRegion(object): - def __init__( - self, - design_parameters, - volume=None, - size=None, - center=mp.Vector3(), - MaterialGrid=None, - ): + def __init__(self, + design_parameters, + volume=None, + size=None, + center=mp.Vector3(), + MaterialGrid=None): self.volume = volume if volume else mp.Volume(center=center, size=size) self.size = self.volume.size self.center = self.volume.center @@ -59,22 +57,20 @@ class OptimizationProblem(object): This is done by the __call__ method. """ - def __init__( - self, - simulation, - objective_functions, - objective_arguments, - design_regions, - frequencies=None, - fcen=None, - df=None, - nf=None, - decay_dt=50, - decay_fields=[mp.Ez], - decay_by=1e-6, - minimum_run_time=0, - maximum_run_time=None, - ): + def __init__(self, + simulation, + objective_functions, + objective_arguments, + design_regions, + frequencies=None, + fcen=None, + df=None, + nf=None, + decay_dt=50, + decay_fields=[mp.Ez], + decay_by=1e-6, + minimum_run_time=0, + maximum_run_time=None): self.sim = simulation @@ -107,12 +103,10 @@ def __init__( fmax = fcen + 0.5 * df fmin = fcen - 0.5 * df dfreq = (fmax - fmin) / (nf - 1) - self.frequencies = np.linspace( - fmin, - fmin + dfreq * nf, - num=nf, - endpoint=False, - ) + self.frequencies = np.linspace(fmin, + fmin + dfreq * nf, + num=nf, + endpoint=False) self.nf = nf if self.nf == 1: @@ -210,12 +204,11 @@ def prepare_forward_run(self): # register design region self.design_region_monitors = [ - self.sim.add_dft_fields( - [mp.Ex, mp.Ey, mp.Ez], - self.frequencies, - where=dr.volume, - yee_grid=True, - ) for dr in self.design_regions + self.sim.add_dft_fields([mp.Ex, mp.Ey, mp.Ez], + self.frequencies, + where=dr.volume, + yee_grid=True) + for dr in self.design_regions ] # store design region voxel parameters @@ -233,16 +226,9 @@ def forward_run(self): # Forward run self.sim.run(until_after_sources=stop_when_dft_decayed( - self.sim, - self.design_region_monitors, - self.decay_dt, - self.decay_fields, - self.fcen_idx, - self.decay_by, - True, - self.minimum_run_time, - self.maximum_run_time, - )) + self.sim, self.design_region_monitors, self.decay_dt, + self.decay_fields, self.fcen_idx, self.decay_by, True, + self.minimum_run_time, self.maximum_run_time)) # record objective quantities from user specified monitors self.results_list = [] @@ -278,8 +264,8 @@ def prepare_adjoint_run(self): for ar in range(len(self.objective_functions)): for mi, m in enumerate(self.objective_arguments): dJ = jacobian(self.objective_functions[ar], - mi)(*self.results_list) - # get gradient of objective w.r.t. monitor + mi)(*self.results_list + ) # get gradient of objective w.r.t. monitor if np.any(dJ): self.adjoint_sources[ar] += m.place_adjoint_source( dJ) # place the appropriate adjoint sources @@ -296,26 +282,18 @@ def adjoint_run(self): # register design flux self.design_region_monitors = [ - self.sim.add_dft_fields( - [mp.Ex, mp.Ey, mp.Ez], - self.frequencies, - where=dr.volume, - yee_grid=True, - ) for dr in self.design_regions + self.sim.add_dft_fields([mp.Ex, mp.Ey, mp.Ez], + self.frequencies, + where=dr.volume, + yee_grid=True) + for dr in self.design_regions ] # Adjoint run self.sim.run(until_after_sources=stop_when_dft_decayed( - self.sim, - self.design_region_monitors, - self.decay_dt, - self.decay_fields, - self.fcen_idx, - self.decay_by, - True, - self.minimum_run_time, - self.maximum_run_time, - )) + self.sim, self.design_region_monitors, self.decay_dt, + self.decay_fields, self.fcen_idx, self.decay_by, True, + self.minimum_run_time, self.maximum_run_time)) # Store adjoint fields for each design set of design variables in array (x,y,z,field_components,frequencies) self.a_E.append([[ @@ -334,12 +312,9 @@ def adjoint_run(self): def calculate_gradient(self): # Iterate through all design regions and calculate gradient self.gradient = [[ - dr.get_gradient( - self.sim, - self.a_E[ar][dri], - self.d_E[dri], - self.frequencies, - ) for dri, dr in enumerate(self.design_regions) + dr.get_gradient(self.sim, self.a_E[ar][dri], self.d_E[dri], + self.frequencies) + for dri, dr in enumerate(self.design_regions) ] for ar in range(len(self.objective_functions))] # Cleanup list of lists @@ -356,13 +331,11 @@ def calculate_gradient(self): # Return optimizer's state to initialization self.current_state = "INIT" - def calculate_fd_gradient( - self, - num_gradients=1, - db=1e-4, - design_variables_idx=0, - filter=None, - ): + def calculate_fd_gradient(self, + num_gradients=1, + db=1e-4, + design_variables_idx=0, + filter=None): ''' Estimate central difference gradients. @@ -399,8 +372,7 @@ def calculate_fd_gradient( fd_gradient_idx = np.random.choice( self.num_design_params[design_variables_idx], num_gradients, - replace=False, - ) + replace=False) for k in fd_gradient_idx: @@ -424,16 +396,9 @@ def calculate_fd_gradient( m.register_monitors(self.frequencies)) self.sim.run(until_after_sources=stop_when_dft_decayed( - self.sim, - self.forward_monitors, - self.decay_dt, - self.decay_fields, - self.fcen_idx, - self.decay_by, - True, - self.minimum_run_time, - self.maximum_run_time, - )) + self.sim, self.forward_monitors, self.decay_dt, + self.decay_fields, self.fcen_idx, self.decay_by, True, + self.minimum_run_time, self.maximum_run_time)) # record final objective function value results_list = [] @@ -459,16 +424,9 @@ def calculate_fd_gradient( # add monitor used to track dft convergence self.sim.run(until_after_sources=stop_when_dft_decayed( - self.sim, - self.forward_monitors, - self.decay_dt, - self.decay_fields, - self.fcen_idx, - self.decay_by, - True, - self.minimum_run_time, - self.maximum_run_time, - )) + self.sim, self.forward_monitors, self.decay_dt, + self.decay_fields, self.fcen_idx, self.decay_by, True, + self.minimum_run_time, self.maximum_run_time)) # record final objective function value results_list = [] @@ -525,17 +483,15 @@ def plot2D(self, init_opt=False, **kwargs): self.sim.plot2D(**kwargs) -def stop_when_dft_decayed( - simob, - mon, - dt, - c, - fcen_idx, - decay_by, - yee_grid=False, - minimum_run_time=0, - maximum_run_time=None, -): +def stop_when_dft_decayed(simob, + mon, + dt, + c, + fcen_idx, + decay_by, + yee_grid=False, + minimum_run_time=0, + maximum_run_time=None): '''Step function that monitors the relative change in DFT fields for a list of monitors. mon ............. a list of monitors diff --git a/python/adjoint/utils.py b/python/adjoint/utils.py index 110a4c4c2..600ca8fa8 100644 --- a/python/adjoint/utils.py +++ b/python/adjoint/utils.py @@ -80,11 +80,11 @@ def gather_monitor_values(monitors: List[ObjectiveQuantity]) -> onp.ndarray: Args: monitors: the mode monitors. - Returns: - a rank-2 ndarray, where the dimensions are (monitor, frequency), of dtype - complex128. Note that these values refer to the mode as oriented (i.e. they - are unidirectional). - """ + Returns: + a rank-2 ndarray, where the dimensions are (monitor, frequency), of dtype + complex128. Note that these values refer to the mode as oriented (i.e. they + are unidirectional). + """ monitor_values = [] for monitor in monitors: monitor_values.append(monitor()) @@ -101,21 +101,21 @@ def gather_design_region_fields( ) -> List[List[onp.ndarray]]: """Collects the design region DFT fields from the simulation. - Args: - simulation: the simulation object. - design_region_monitors: the installed design region monitors. - frequencies: the frequencies to monitor. - - Returns: - A list of lists. Each entry (list) in the overall list corresponds one-to- - one with a declared design region. For each such contained list, the - entries correspond to the field components that are monitored. The entries - are ndarrays of rank 4 with dimensions (freq, x, y, (z-or-pad)). - - The design region fields are sampled on the *Yee grid*. This makes them - fairly awkward to inspect directly. Their primary use case is supporting - gradient calculations. - """ + Args: + simulation: the simulation object. + design_region_monitors: the installed design region monitors. + frequencies: the frequencies to monitor. + + Returns: + A list of lists. Each entry (list) in the overall list corresponds one-to- + one with a declared design region. For each such contained list, the + entries correspond to the field components that are monitored. The entries + are ndarrays of rank 4 with dimensions (freq, x, y, (z-or-pad)). + + The design region fields are sampled on the *Yee grid*. This makes them + fairly awkward to inspect directly. Their primary use case is supporting + gradient calculations. + """ fwd_fields = [] for monitor in design_region_monitors: fields_by_component = [] @@ -134,17 +134,17 @@ def validate_and_update_design( design_variables: Iterable[onp.ndarray]) -> None: """Validate the design regions and variables. - In particular the design variable should be 1,2,3-D and the design region - shape should match the design variable shape after dimension expansion. - The arguments are modified in place. + In particular the design variable should be 1,2,3-D and the design region + shape should match the design variable shape after dimension expansion. + The arguments are modified in place. Args: design_regions: List of mpa.DesignRegion, design_variables: Iterable with numpy arrays representing design variables. - Raises: - ValueError if the validation of dimensions fails. - """ + Raises: + ValueError if the validation of dimensions fails. + """ for i, (design_region, design_variable) in enumerate(zip(design_regions, design_variables)): diff --git a/python/adjoint/wrapper.py b/python/adjoint/wrapper.py index d36d1b89d..26d37e3df 100644 --- a/python/adjoint/wrapper.py +++ b/python/adjoint/wrapper.py @@ -64,48 +64,46 @@ def loss(x): class MeepJaxWrapper: """Wraps a Meep simulation object into a JAX-differentiable callable. - Attributes: - simulation: the pre-configured Meep simulation object. - sources: a list of Meep sources for the forward simulation. - monitors: a list of eigenmode coefficient monitors from the `meep.adjoint` - module. - design_regions: a list of design regions from the `meep.adjoint` module. - frequencies: the list of frequencies, in normalized Meep units. - measurement_interval: the time interval between DFT field convergence - measurements, in Meep time units. The default value is 50. - dft_field_components: a list of Meep field components, such as `mp.Ex`, - `mp.Hy`, etc, whose DFT will be monitored for convergence to stop the - simulation. The default is `mp.Ez`. - dft_threshold: the threshold for DFT field convergence. Once the norm of the - change in the fields (the maximum over all design regions and field - components) is less than this value, the simulation will be stopped. The - default value is 1e-6. - minimum_run_time: the minimum run time of the simulation, in Meep time - units. The default value is 0. - maximum_run_time: the maximum run time of the simulation, in Meep time - units. The default value is infinity. - until_after_sources: whether `maximum_run_time` should be ignored until the - sources have turned off. This parameter specifies whether `until` or - `until_after_sources` is used. See - https://meep.readthedocs.io/en/latest/Python_User_Interface/#Simulation + Attributes: + simulation: the pre-configured Meep simulation object. + sources: a list of Meep sources for the forward simulation. + monitors: a list of eigenmode coefficient monitors from the `meep.adjoint` + module. + design_regions: a list of design regions from the `meep.adjoint` module. + frequencies: the list of frequencies, in normalized Meep units. + measurement_interval: the time interval between DFT field convergence + measurements, in Meep time units. The default value is 50. + dft_field_components: a list of Meep field components, such as `mp.Ex`, + `mp.Hy`, etc, whose DFT will be monitored for convergence to stop the + simulation. The default is `mp.Ez`. + dft_threshold: the threshold for DFT field convergence. Once the norm of the + change in the fields (the maximum over all design regions and field + components) is less than this value, the simulation will be stopped. The + default value is 1e-6. + minimum_run_time: the minimum run time of the simulation, in Meep time + units. The default value is 0. + maximum_run_time: the maximum run time of the simulation, in Meep time + units. The default value is infinity. + until_after_sources: whether `maximum_run_time` should be ignored until the + sources have turned off. This parameter specifies whether `until` or + `until_after_sources` is used. See + https://meep.readthedocs.io/en/latest/Python_User_Interface/#Simulation for more information. The default is true. - """ + """ _log_fn = print - def __init__( - self, - simulation: mp.Simulation, - sources: List[mp.Source], - monitors: List[EigenmodeCoefficient], - design_regions: List[DesignRegion], - frequencies: List[float], - measurement_interval: float = 50.0, - dft_field_components: Tuple[int, ...] = (mp.Ez, ), - dft_threshold: float = 1e-6, - minimum_run_time: float = 0, - maximum_run_time: float = onp.inf, - until_after_sources: bool = True, - ): + def __init__(self, + simulation: mp.Simulation, + sources: List[mp.Source], + monitors: List[EigenmodeCoefficient], + design_regions: List[DesignRegion], + frequencies: List[float], + measurement_interval: float = 50.0, + dft_field_components: Tuple[int, ...] = (mp.Ez, ), + dft_threshold: float = 1e-6, + minimum_run_time: float = 0, + maximum_run_time: float = onp.inf, + until_after_sources: bool = True): self.simulation = simulation self.sources = sources self.monitors = monitors @@ -124,20 +122,20 @@ def __init__( def __call__(self, designs: List[jnp.ndarray]) -> jnp.ndarray: """Performs a Meep simulation, taking a list of designs and returning mode overlaps. - Args: - designs: a list of design variables as 1D, 2D, or 3D JAX arrays. Valid shapes for - design variables are (Nx, Ny, Nz) where Nx{y,z} match the elements of the - `grid_size` constructor argument of Meep's `MaterialGrid` used for the - corresponding design region. Singleton dimensions of the `grid_size` may be - omitted from the corresponding design variable. For example, a design variable - with a shape of either (10, 20) or (10, 20, 1) would be compatible with a - `grid_size` of (10, 20, 1). Similarly, a design variable with shapes of (25,), - (25, 1), or (25, 1, 1) would be compatible with a `grid_size` of (25, 1, 1). - - Returns: - a complex-valued JAX ndarray of differentiable mode monitor overlaps values with - a shape of (num monitors, num frequencies). - """ + Args: + designs: a list of design variables as 1D, 2D, or 3D JAX arrays. Valid shapes for + design variables are (Nx, Ny, Nz) where Nx{y,z} match the elements of the + `grid_size` constructor argument of Meep's `MaterialGrid` used for the + corresponding design region. Singleton dimensions of the `grid_size` may be + omitted from the corresponding design variable. For example, a design variable + with a shape of either (10, 20) or (10, 20, 1) would be compatible with a + `grid_size` of (10, 20, 1). Similarly, a design variable with shapes of (25,), + (25, 1), or (25, 1, 1) would be compatible with a `grid_size` of (25, 1, 1). + + Returns: + a complex-valued JAX ndarray of differentiable mode monitor overlaps values with + a shape of (num monitors, num frequencies). + """ return self._simulate_fn(designs) def _reset_convergence_measurement(self, @@ -202,19 +200,15 @@ def _run_adjoint_simulation(self, monitor_values_grad): self._reset_convergence_measurement(design_region_monitors) self.simulation.run(**sim_run_args) - return utils.gather_design_region_fields( - self.simulation, - design_region_monitors, - self.frequencies, - ) + return utils.gather_design_region_fields(self.simulation, + design_region_monitors, + self.frequencies) - def _calculate_vjps( - self, - fwd_fields, - adj_fields, - design_variable_shapes, - sum_freq_partials=True, - ): + def _calculate_vjps(self, + fwd_fields, + adj_fields, + design_variable_shapes, + sum_freq_partials=True): """Calculates the VJP for a given set of forward and adjoint fields.""" return utils.calculate_vjps( self.simulation, @@ -298,19 +292,19 @@ def _are_dfts_converged(self, sim: mp.Simulation) -> bool: def _simulation_run_callback(self, sim: mp.Simulation) -> bool: """A callback function returning `True` when the simulation should stop. - This is a step function that gets called at each time step of the Meep - simulation, taking a Meep simulation object as an input and returning `True` - when the simulation should be terminated and returning `False` when the - simulation should continue. The resulting step function is designed to be - used with the `until` and `until_after_sources` arguments of the Meep - simulation `run()` routine. + This is a step function that gets called at each time step of the Meep + simulation, taking a Meep simulation object as an input and returning `True` + when the simulation should be terminated and returning `False` when the + simulation should continue. The resulting step function is designed to be + used with the `until` and `until_after_sources` arguments of the Meep + simulation `run()` routine. - Args: - sim: a Meep simulation object. + Args: + sim: a Meep simulation object. - Returns: - a boolean indicating whether the simulation should be terminated. - """ + Returns: + a boolean indicating whether the simulation should be terminated. + """ current_meep_time = sim.round_time() if current_meep_time <= self._last_measurement_meep_time + self.measurement_interval: return False From 1b0d5462a85389e559ddfc72b925aac6fde4666d Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Wed, 25 Aug 2021 10:49:59 -0400 Subject: [PATCH 14/15] rebase cont. --- python/adjoint/optimization_problem.py | 33 -------------------------- 1 file changed, 33 deletions(-) diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py index d6267d3a3..fe5e20da0 100644 --- a/python/adjoint/optimization_problem.py +++ b/python/adjoint/optimization_problem.py @@ -57,22 +57,6 @@ class OptimizationProblem(object): This is done by the __call__ method. """ -<<<<<<< HEAD - def __init__(self, - simulation, - objective_functions, - objective_arguments, - design_regions, - frequencies=None, - fcen=None, - df=None, - nf=None, - decay_dt=50, - decay_fields=[mp.Ez], - decay_by=1e-6, - minimum_run_time=0, - maximum_run_time=None): -======= def __init__( self, simulation, @@ -90,7 +74,6 @@ def __init__( minimum_run_time=0, maximum_run_time=None, ): ->>>>>>> 21c94c0d25ac205168819fa98b10de2dce9dc0e4 self.sim = simulation @@ -225,13 +208,6 @@ def prepare_forward_run(self): # register design region self.design_region_monitors = [ -<<<<<<< HEAD - self.sim.add_dft_fields([mp.Ex, mp.Ey, mp.Ez], - self.frequencies, - where=dr.volume, - yee_grid=True) - for dr in self.design_regions -======= self.sim.add_dft_fields( [mp.Ex, mp.Ey, mp.Ez], self.frequencies, @@ -239,7 +215,6 @@ def prepare_forward_run(self): yee_grid=True, decimation_factor=self.decimation_factor, ) for dr in self.design_regions ->>>>>>> 21c94c0d25ac205168819fa98b10de2dce9dc0e4 ] # store design region voxel parameters @@ -314,13 +289,6 @@ def adjoint_run(self): # register design flux self.design_region_monitors = [ -<<<<<<< HEAD - self.sim.add_dft_fields([mp.Ex, mp.Ey, mp.Ez], - self.frequencies, - where=dr.volume, - yee_grid=True) - for dr in self.design_regions -======= self.sim.add_dft_fields( [mp.Ex, mp.Ey, mp.Ez], self.frequencies, @@ -328,7 +296,6 @@ def adjoint_run(self): yee_grid=True, decimation_factor=self.decimation_factor, ) for dr in self.design_regions ->>>>>>> 21c94c0d25ac205168819fa98b10de2dce9dc0e4 ] # Adjoint run From 94ae1508b72861d5c5adbe6b84bd8c4b57e6bc03 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Wed, 25 Aug 2021 12:05:09 -0400 Subject: [PATCH 15/15] fix None bug --- python/adjoint/objective.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index dac95575e..3a90692dc 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -43,7 +43,7 @@ def place_adjoint_source(self, dJ): def get_evaluation(self): """Evaluates the objective quantity.""" - if self._eval: + if self._eval is not None: return self._eval else: raise RuntimeError(