From 762938960ab1d4ff11be02687a3f753d8ced8797 Mon Sep 17 00:00:00 2001 From: Edoardo Vergerio Date: Sun, 24 May 2026 17:47:07 +0200 Subject: [PATCH] modified version of metapeg, using epigraph and linewidth/spacing constraints --- metapeg.ipynb | 7 +- metapeg.py | 19 +- metapeg_v2.py | 731 ++++++++++++++++++++++++++++++++++++++++++++++++++ misc/logo.png | Bin 0 -> 7480 bytes misc/peg.png | Bin 0 -> 4250 bytes 5 files changed, 750 insertions(+), 7 deletions(-) create mode 100644 metapeg_v2.py create mode 100644 misc/logo.png create mode 100644 misc/peg.png diff --git a/metapeg.ipynb b/metapeg.ipynb index 4486bab..cfb30c4 100644 --- a/metapeg.ipynb +++ b/metapeg.ipynb @@ -631,6 +631,11 @@ } ], "metadata": { + "kernelspec": { + "display_name": "mp", + "language": "python", + "name": "python3" + }, "language_info": { "codemirror_mode": { "name": "ipython", @@ -641,7 +646,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": 3 + "version": "3.13.12" } }, "nbformat": 4, diff --git a/metapeg.py b/metapeg.py index b64276b..4d4c51c 100644 --- a/metapeg.py +++ b/metapeg.py @@ -59,7 +59,7 @@ design_region_width = 4.0 design_region_length = 4.0 design_region_height = (wavelength/np.sqrt(subs_index))/2 #half-wavelength in the substrate -resolution = 20 # pixels/um +resolution = 5 # pixels/um design_region_resolution = int(resolution) #number of pixels in the design region @@ -133,6 +133,7 @@ output_plane_xy = mp.Volume( # %% ######SIMULATION###### +sim_run_time = 100 kpoint = mp.Vector3() sim = mp.Simulation( cell_size = cell_size, @@ -261,7 +262,7 @@ norm_near2far = norm_sim.add_near2far(fcen,0,1,*NearRegions) norm_sim.plot2D(output_plane=output_plane,show_sources=True, show_boundaries=True, show_geometry=True) plt.savefig("results/normalization_simulation_geometry.pdf") -norm_sim.run(until=1500) +norm_sim.run(until=sim_run_time) ref_fields = np.array([norm_sim.get_farfield(norm_near2far, point) for point in far_z_plane]) #TODO this todo is linked to the one above. Need to check if the far-field need to be flipped or something #at this time, the two definition are consistent with each other. @@ -424,8 +425,8 @@ ub = np.ones((Nx * Ny,)) cur_beta = 4 #starting beta beta_scale = 2 #beta scaling (update) factor -num_betas = 6 #number of betas to test -update_factor = 12 #number of iterations per beta +num_betas = 3 #number of betas to test +update_factor = 1 #number of iterations per beta #total number of simulations = num_betas * update_factor * 2 (forward + adjoint) ftol = 1e-5 @@ -514,6 +515,12 @@ full_output_plane = mp.Volume( center=mp.Vector3(0,0,0), size=mp.Vector3(Sx,0,full_Sz)) -opt.sim.plot2D(output_plane=full_output_plane,fields=mp.Ex) +design_output_plane = mp.Volume( + center=mp.Vector3(0,0,10), + size=mp.Vector3(design_region_width,design_region_length,0) +) + +#understand how to retrieve the far-field data at the design output plane +opt.sim.plot2D(output_plane=full_output_plane,fields=mp.Ex, show_sources=True, show_boundaries=True, show_geometry=True) plt.savefig("results/full_simulation_geometry.pdf") -opt.sim.run(until=1500) \ No newline at end of file +opt.sim.run(until=sim_run_time) \ No newline at end of file diff --git a/metapeg_v2.py b/metapeg_v2.py new file mode 100644 index 0000000..651091e --- /dev/null +++ b/metapeg_v2.py @@ -0,0 +1,731 @@ +from typing import List, NamedTuple, Tuple + +from autograd import numpy as npa, tensor_jacobian_product, grad +import matplotlib.pyplot as plt +import meep as mp +import meep.adjoint as mpa +import nlopt +import numpy as np +import os +import xarray as xr +from PIL import Image + +mp.verbosity(0) + +if not os.path.exists("results"): + os.makedirs("results") + +###### PHYSICS SETUP ###### +RESOLUTION = 20 #50 #pixels per unit length +MAX_RUN_TIME = 1100 #meep units +WAVELENGTH_MIN_UM = 1.50 +WAVELENGTH_MAX_UM = 1.60 +PML_UM = 1.0 #PML thickness +BUFFER = 2.0 +DESIGN_WAVELENGTHS_UM = [1.50,1.55,1.60] #wavelengths at which to optimize the device +DESIGN_REGION_UM = mp.Vector3(5,5,1) +DESIGN_REGION_RESOLUTION = int(2 * RESOLUTION) +DESIGN_REGION_CENTER = mp.Vector3(0, 0, DESIGN_REGION_UM.z / 2) +NX_DESIGN_GRID = int(DESIGN_REGION_UM.x * DESIGN_REGION_RESOLUTION) + 1 +NY_DESIGN_GRID = int(DESIGN_REGION_UM.y * DESIGN_REGION_RESOLUTION) + 1 +MIN_LENGTH_UM = 0.15 #150 nm minimum length scale +SILICON = mp.Medium(index = 3.5) +SILICON_DIOXIDE = mp.Medium(index = 1.5) +AIR = mp.Medium(index = 1.0) + +SIGMOID_THRESHOLD_INTRINSIC = 0.5 +SIGMOID_THRESHOLD_EROSION = 0.75 +SIGMOID_THRESHOLD_DILATION = 1 - SIGMOID_THRESHOLD_EROSION + +cell_um = mp.Vector3( + DESIGN_REGION_UM.x + 2*BUFFER + 2*PML_UM, + DESIGN_REGION_UM.y + 2*BUFFER + 2*PML_UM, + DESIGN_REGION_UM.z + 2*BUFFER + 2*PML_UM +) + +SUBSTRATE_THICKNESS = cell_um.z/2 +SUBSTRATE_SIZE = mp.Vector3(cell_um.x, cell_um.y, SUBSTRATE_THICKNESS) +SUBSTRATE_CENTER = mp.Vector3(0, 0, - SUBSTRATE_THICKNESS / 2) + +physical_domains_size = mp.Vector3( + DESIGN_REGION_UM.x + 2*BUFFER, + DESIGN_REGION_UM.y + 2*BUFFER, + DESIGN_REGION_UM.z + 2*BUFFER +) + +filter_radius_um = mpa.get_conic_radius_from_eta_e( + MIN_LENGTH_UM, SIGMOID_THRESHOLD_EROSION +) + +frequency_min = 1 / WAVELENGTH_MAX_UM +frequency_max = 1 / WAVELENGTH_MIN_UM +frequency_center = 0.5 * (frequency_min + frequency_max) +wavelength_center = 1 / frequency_center #remember that c = 1 in meep units +frequency_width = frequency_max - frequency_min + +SOURCE_CENTER = mp.Vector3(0,0,-BUFFER) +SOURCE_SIZE = mp.Vector3(physical_domains_size.x,physical_domains_size.y,0) + +NEAR_REGION_MONITOR_CENTER = mp.Vector3(0,0, DESIGN_REGION_UM.z+0.5) +NEAR_REGION_MONITOR_SIZE = mp.Vector3(physical_domains_size.x,physical_domains_size.y,0) + +ff_monitor_center = mp.Vector3(0,0,DESIGN_REGION_UM.z + BUFFER) +ff_monitor_size = mp.Vector3(DESIGN_REGION_UM.x, DESIGN_REGION_UM.y, 0) +xs = ys = np.linspace(-ff_monitor_size.x/2, ff_monitor_size.x/2, NX_DESIGN_GRID) +#collection of points in the far-field monitor +ff_points = [mp.Vector3(x_p,y_p,ff_monitor_center.z) for x_p in xs for y_p in ys] + +stop_cond = mp.stop_when_fields_decayed(50, mp.Ex, ff_monitor_center, 1e-6) + +pml_layers = [mp.PML(PML_UM)] +frequencies = [1/wavelength for wavelength in DESIGN_WAVELENGTHS_UM] +num_wavelengths = len(DESIGN_WAVELENGTHS_UM) + +view_2D_plane = mp.Volume( + center=mp.Vector3(0,0,0), + size=mp.Vector3(cell_um.x,0,cell_um.z) +) + +def str_from_list(list_: List[float]) -> str: + return "[" + ", ".join(f"{val:.4f}" for val in list_) + "]" + +def filter_and_project( + weights: np.ndarray, + sigmoid_threshold: float, + sigmoid_bias: float +) -> np.ndarray: + '''Function to filter and project design variables. Differentiable by autograd + + Args: + weights: design weights as a flattened (1D) array + sigmoid_threshold: erosion/dilation parameter for projection + sigmoid_bias: bias parameter for the projection. 0 is no projection + Returns: + The flattened filtered/projected design weights + ''' + + weights_filtered = mpa.conic_filter( + weights, + filter_radius_um, + DESIGN_REGION_UM.x, + DESIGN_REGION_UM.y, + DESIGN_REGION_RESOLUTION, + ) + + if sigmoid_bias == 0: + return weights_filtered.flatten() + else: + weights_projected = mpa.tanh_projection( + weights_filtered.reshape(NX_DESIGN_GRID,NY_DESIGN_GRID), + sigmoid_bias, + sigmoid_threshold, + ) + return weights_projected.flatten() + +def obj_fun(epigraph_and_weights: np.ndarray, grad: np.ndarray)-> float: + '''Objective function for the epigraph formulation + + Args: + epigraph_and_weights: 1D array containing epigraph variable + (at index 0) and design weigts (remaining elements) + grad: the gradient modified in place of the objective function (flattened 1D array) + Returns: + The scalar epigraph variable (the objective function) + ''' + + epigraph = epigraph_and_weights[0] + + if grad.size > 0: + grad[0] = 1 + grad[1:] = 0 + + return epigraph + +def epigraph_constraint( + result: np.ndarray, + epigraph_and_weights: np.ndarray, + gradient: np.ndarray, + sigmoid_threshold: float, + sigmoid_bias: float, + use_epsavg: bool, +) -> None: + '''Constraint function for the epigraph formulation + + Args: + result: 1D array, modified in place, with the results + of this constraint evaluation + epigraph_and_weights: 1D array of epigraph (first element) + and design weights (other elements) + gradient: the jacobian matrix (backpropagated) with dimensions + (1+NX_DESIGN_GRID*NY_DESIGN_GRID, 2 * num_wavelengths), modified in place. + these gradients are used by the optimization algorithm, + because, being backpropagated, they tell us how the constraint function behaves + with respect to changes in the design parameters. + sigmoid_threshold: erosion/dilation parameter for projection. + sigmoid_bias: bias parameter for projection. + use_epsavg: whether to use subpixel smoothing. + ''' + epigraph = epigraph_and_weights[0] + weights = epigraph_and_weights[1:] + + obj_val, grad = opt( + [ + filter_and_project( + weights, sigmoid_threshold, 0 if use_epsavg else sigmoid_bias) + ] + ) + + #in our case there is only one objective function, + #which needs to be evaluated at multiple wavelengths. + #for this reason, obj_val is a 1-element list (to be verified, + # it might be that the list is directly replaced by the content) + # that contains an array of objective function values at each wavelength. + + #for each wavelength, backpropagate + for k in range(num_wavelengths): + grad[:,k] = tensor_jacobian_product(filter_and_project, 0)( + weights, + sigmoid_threshold, + sigmoid_bias, + grad[:,k], + ) + + #modify in place the gradient result + if gradient.size > 0: + gradient[:,0] = -1 #gradient w.r.t epigraph + gradient[:,1:] = grad.T + + #modify in place the constraint result + result[:] = np.real(obj_val) - epigraph + + print( + f"iteration:, {cur_iter[0]:3d}, sigmoid_bias: {sigmoid_bias:2d}, " + f"epigraph: {epigraph:.5f}, obj. func.: {obj_val}, " + f"epigraph constraint: {str_from_list(result)}" + ) + + cur_iter[0] = cur_iter[0] + 1 + +def line_width_and_spacing_constraint( + result: np.ndarray, + epigraph_and_weights: np.ndarray, + gradient: np.ndarray, + sigmoid_bias: float, +) -> float: + '''Constraint function for minimum line width and spacing + + Args: + result: evaluation of this constraint function, modified in place + epigraph_and_weights: 1D array of epigraph (first element) + and design weights (other elements) + gradient: the Jacobian matrix, modified in place + sigmoid_bias: bias parameter for projection + + Returns: + The constraint value result + ''' + + epigraph = epigraph_and_weights[0] + weights = epigraph_and_weights[1:] + + a1 = 1e-3 + b1 = 0 + gradient[:,0] = -a1 + + filter_func = lambda a: mpa.conic_filter( + a.reshape(NX_DESIGN_GRID,NY_DESIGN_GRID), + filter_radius_um, + DESIGN_REGION_UM.x, + DESIGN_REGION_UM.y, + DESIGN_REGION_RESOLUTION, + ) + + threshold_func = lambda a: mpa.tanh_projection( + a.reshape(NX_DESIGN_GRID,NY_DESIGN_GRID), + sigmoid_bias, + SIGMOID_THRESHOLD_INTRINSIC, + ) + + c0 = 1e7 * (filter_radius_um * 1 / RESOLUTION) ** 4 + + #constraints as lambda functions + M1 = lambda a: mpa.constraint_solid( + a, + c0, + SIGMOID_THRESHOLD_EROSION, + filter_func, + threshold_func, + 1 + ) + + M2 = lambda a: mpa.constraint_void( + a, + c0, + SIGMOID_THRESHOLD_DILATION, + filter_func, + threshold_func, + 1 + ) + + #gradients of the constraints w.r.t design parameters + g1 = grad(M1)(weights) + g2 = grad(M2)(weights) + + result[0] = M1(weights) - a1 * epigraph - b1 + result[1] = M2(weights) - a1 * epigraph - b1 + + gradient[0, 1:] = g1.flatten() + gradient[1, 1:] = g2.flatten() + + #these are the epigraph values for the constraints + t1 = (M1(weights) - b1) / a1 + t2 = (M2(weights) - b1) / a1 + + print(f"line_width_and_spacing_constraint:, {result[0]}, {result[1]}, {t1}, {t2}") + + return max(t1, t2) + +def get_pattern() -> np.ndarray: + '''Imports the intensity pattern from file''' + pattern_name = 'peg.png' + pattern = Image.open('misc/' + pattern_name).convert('L') + pattern = np.array(pattern).astype(float) + + pattern = pattern - np.min(pattern) + pattern = pattern / np.max(pattern) + + return pattern + +def intensity_desired_fn_pattern( + xs: list, + ys: list, + rescale: float = 0.5 +) -> np.ndarray: + '''Returns the pattern intensity value at the specified + (xs,ys) locations, with a scale factor + + Args: + xs: list of x coordinates + ys: list of y coordinates + rescale: scale factor for the pattern + Returns: + The flattened pattern intensity value at the specified (xs,ys) locations, with a scale factor + ''' + pattern_values = get_pattern() + pattern_values = 1 - np.rot90(np.rot90(np.rot90(pattern_values))) + + nx,ny = pattern_values.shape + xs_pattern = np.linspace(rescale*min(xs), rescale*max(xs), nx) + ys_pattern = np.linspace(rescale*min(ys), rescale*max(ys), ny) + pattern_dataArray = xr.DataArray(pattern_values, coords=dict(x=xs_pattern, y=ys_pattern)) + pattern_interp = pattern_dataArray.interp(x=xs, y=ys) + + pattern_final = npa.nan_to_num(pattern_interp.values, nan=npa.min(pattern_interp)) + pattern_to_return = npa.nan_to_num(pattern_interp.values, nan=npa.min(pattern_interp)).flatten() + + #make a plot before returning + fig, ax = plt.subplots() + im = ax.imshow( + pattern_final, + extent=(min(xs), max(xs), min(ys), max(ys)), + origin='lower', + cmap='inferno', + ) + ax.set_xlabel('x (um)') + ax.set_ylabel('y (um)') + ax.set_title('Target Intensity Pattern') + fig.colorbar(im, ax=ax, label='Intensity (a.u.)') + plt.savefig("results/target_intensity_pattern.pdf") + + return pattern_to_return + +def normalization_sim() -> np.ndarray: + ''' Computes the Far-Field pattern at the monitor location and + without the design region. Used for normalization purposes. + + Returns: + A 2D-array [Nx*Ny,num_wavelengths] of reference intensities, for each wavelength, + at the monitor location. + ''' + + sources = [ + mp.Source( + src=mp.GaussianSource(frequency_center, fwidth=frequency_width), + component=mp.Ex, + center=SOURCE_CENTER, + size=SOURCE_SIZE, + ) + ] + + substrate_block = mp.Block( + center = SUBSTRATE_CENTER, + size = SUBSTRATE_SIZE, + material = SILICON_DIOXIDE, + ) + + geometry = [#combine geometries + substrate_block, + ] + + norm_sim = mp.Simulation( + resolution=RESOLUTION, + default_material=AIR, + cell_size=cell_um, + sources=sources, + geometry=geometry, + boundary_layers=pml_layers, + k_point=mp.Vector3(), + ) + + fig, ax = plt.subplot() + norm_sim.plot2D(output_plane=view_2D_plane, ax=ax) + fig.savefig("results/normalization_sim_setup.pdf") + + + NearRegions = [mp.Near2FarRegion( + center=NEAR_REGION_MONITOR_CENTER, + size=NEAR_REGION_MONITOR_SIZE, + weight=+1, + )] + + norm_near2far = norm_sim.add_near2far(frequencies, *NearRegions) + + norm_sim.run(until_after_sources=stop_cond) + + ref_fields = np.array( + [norm_sim.get_farfield(norm_near2far, point) for point in ff_points] + ) + + #select only the Ex component for all frequencies, + # which is the one we will optimize for + ref_fields_x_f1 = ref_fields[:,0] + ref_fields_x_f2 = ref_fields[:,6] + ref_fields_x_f3 = ref_fields[:,12] + #concatenate the Ex component for all frequencies into a single array + #the rows are the points, the columns are the frequencies + ref_fields_x = npa.stack((ref_fields_x_f1, ref_fields_x_f2, ref_fields_x_f3), axis=-1) + + ref_intensity_x = npa.power(npa.abs(ref_fields_x),2) + return ref_intensity_x + +def intensity_from_farfields(FarFields) -> np.ndarray: + '''Computes the intensity pattern at the monitor location from the simulated far-fields + + Args: + FarFields: the simulated far-field patterns at the farfield + monitor location, for each wavelength + Returns: + A 2D-array [Nx*Ny,num_wavelengths] of intensities, for each wavelength, + at the monitor location. + ''' + #select only the Ex component for all frequencies, + # which is the one we will optimize for. + #Reshape into matrix form (Nx,Ny) + #print(FarFields.shape) + #fields_x_f1 = FarFields[:,0,0].reshape(NX_DESIGN_GRID, NY_DESIGN_GRID) + #fields_x_f2 = FarFields[:,1,0].reshape(NX_DESIGN_GRID, NY_DESIGN_GRID) + #fields_x_f3 = FarFields[:,2,0].reshape(NX_DESIGN_GRID, NY_DESIGN_GRID) + #concatenate the Ex component for all frequencies into a single array + #the rows are the points, the columns are the frequencies + #make a 3d arrray with dimensions (Nx,Ny,num_wavelengths) + #fields_stack_3d = npa.stack((fields_x_f1, fields_x_f2, fields_x_f3), axis=-1) + #fields_x = np.stack((fields_x_f1, fields_x_f2, fields_x_f3), axis=-1) + + intensity_x = npa.abs(FarFields[:,:,0]) ** 2 + + #make a plot before returning of the central wavelength + #for the title use cur_iter + intensity_x_f2 = intensity_x[:,1].reshape(NX_DESIGN_GRID, NY_DESIGN_GRID) + fig, ax = plt.subplots() + im = ax.imshow( + intensity_x_f2, + extent=(min(xs), max(xs), min(ys), max(ys)), + origin='lower', + cmap='inferno', + ) + ax.set_xlabel('x (um)') + ax.set_ylabel('y (um)') + ax.set_title('Simulated Intensity Pattern at Monitor (central wavelength) - iteration ' + str(cur_iter[0])) + fig.colorbar(im, ax=ax, label='Intensity (a.u.)') + plt.savefig("results/simulated_intensity_pattern_monitor"+str(cur_iter[0])+".png") + + return intensity_x + +def intensity_optimization( + ref_intensity: np.ndarray, + use_damping:bool, + use_epsavg: bool, + sigmoid_bias: float, +)-> mpa.OptimizationProblem: + '''Sets up the optimization problem for the intensity pattern matching + + Args: + ref_intensity: 2D array of reference intensities at the monitor location, for each wavelength + use_damping: whether to use damping in the optimization algorithm + use_epsavg: whether to use subpixel smoothing in the optimization + sigmoid_bias: bias parameter for projection + Returns: + The optimization problem class instance + ''' + matgrid = mp.MaterialGrid( + grid_size = mp.Vector3(NX_DESIGN_GRID, NY_DESIGN_GRID,0), + medium1 = AIR, + medium2 = SILICON, + weights = np.ones((NX_DESIGN_GRID, NY_DESIGN_GRID)), + beta = sigmoid_bias if use_epsavg else 0, + do_averaging = True if use_epsavg else False, + damping = 0.02*2*np.pi*frequency_center if use_damping else 0, + ) + + matgrid_region = mpa.DesignRegion( + matgrid, + volume=mp.Volume(center=DESIGN_REGION_CENTER, + size=DESIGN_REGION_UM + ) + ) + + matgrid_block = mp.Block( + center=matgrid_region.center, + size=matgrid_region.size, + material=matgrid, + ) + + substrate_block = mp.Block( + center = SUBSTRATE_CENTER, + size = SUBSTRATE_SIZE, + material = SILICON_DIOXIDE, + ) + + geometry = [ + matgrid_block, + substrate_block, + ] + + sources = [ + mp.Source( + src=mp.GaussianSource(frequency_center, fwidth=frequency_width), + component=mp.Ex, + size=SOURCE_SIZE, + center=SOURCE_CENTER, + ) + ] + + sim = mp.Simulation( + resolution=RESOLUTION, + default_material=AIR, + cell_size=cell_um, + sources=sources, + geometry=geometry, + boundary_layers=pml_layers, + k_point=mp.Vector3(), + ) + + fig, ax = plt.subplot() + sim.plot2D(output_plane=view_2D_plane,show_monitors=True,ax=ax) + fig.savefig("results/optimization_sim_setup.pdf") + + NearRegions = [ + mp.Near2FarRegion( + center=NEAR_REGION_MONITOR_CENTER, + size=NEAR_REGION_MONITOR_SIZE, + weight=+1, + ) + ] + + FarFields = mpa.Near2FarFields( + sim=sim, + Near2FarRegions=NearRegions, + far_pts=ff_points, + ) + + obj_list = [FarFields] + + #2D array of target intensity at the monitor location + target_intensity = intensity_desired_fn_pattern(xs, ys) + target_intensity_3d = npa.tile(target_intensity[:,npa.newaxis], (1,num_wavelengths)) + target_intensity_mean = npa.mean(target_intensity) + + def J1(FarFields): + '''Objective function. Computes the distance between simulated and target intensity patterns. + Args: + FarFields: the simulated far-field patterns at the monitor location, for each wavelength + Returns: + The distance between simulated and target intensity patterns + ''' + + measured_intensity = intensity_from_farfields(FarFields) + + measured_intensity_mean = npa.mean(measured_intensity, axis=0) + + norm_measured_intensity = measured_intensity / ref_intensity * target_intensity_mean + norm_measured_intensity = npa.clip(norm_measured_intensity, a_min=0, a_max=1) + + distance = norm_measured_intensity - target_intensity_3d + + return npa.power(npa.sum(distance**2, axis=0),0.5) + + + opt = mpa.OptimizationProblem( + simulation=sim, + objective_functions=[J1], + objective_arguments=obj_list, + design_regions=[matgrid_region], + frequencies=frequencies, + maximum_run_time=MAX_RUN_TIME + ) + + fig, ax = plt.subplot() + opt.plot2D(True,output_plane=view_2D_plane, ax=ax) + fig.savefig("results/OPT_PROB_SETUP.pdf") + + return opt + +if __name__ == "__main__": + print("Hello there! Starting normalization") + ref_intensity = normalization_sim() + + num_weights = NX_DESIGN_GRID * NY_DESIGN_GRID + + epigraph_and_weights = np.ones(shape=num_weights) * 0.5 + epigraph_and_weights = np.insert(epigraph_and_weights, 0, 1.2) + + weights_lower_bound = np.zeros(num_weights) + weights_upper_bound = np.ones(num_weights) + epigraph_and_weights_lower_bound = np.insert(weights_lower_bound, 0, -np.inf) + epigraph_and_weights_upper_bound = np.insert(weights_upper_bound, 0, +np.inf) + + objfunc_history = [] + epivar_history = [] + cur_iter = [0] + + sigmoid_bias_threshold = 64 + + sigmoid_biases = [8, 16, 32, 64, 128, 256] + #max_evals = [80, 80, 100, 120, 120, 100] + max_evals = [10, 10, 10, 10, 10, 10] + epigraph_tolerance = np.array([1e-4]*num_wavelengths) + tolerance_width_and_spacing = np.array([1e-8]*2) + + for sigmoid_bias, max_eval in zip(sigmoid_biases, max_evals): + print(f"Starting optimization epoch with sigmoid bias {sigmoid_bias} and max evals {max_eval}") + #the optimizer need to work with the weights and epigraph (+1) + solver = nlopt.opt(nlopt.LD_MMA, num_weights + 1) + solver.set_lower_bounds(epigraph_and_weights_lower_bound) + solver.set_upper_bounds(epigraph_and_weights_upper_bound) + solver.set_min_objective(obj_fun) + solver.set_maxeval(max_eval) + solver.set_param("dual_ftol_rel", 1e-7) + solver.add_inequality_mconstraint( + lambda result_, epigraph_and_weights_, grad_: epigraph_constraint( + result_, + epigraph_and_weights_, + grad_, + SIGMOID_THRESHOLD_INTRINSIC, + sigmoid_bias, + use_epsavg=False if sigmoid_bias < sigmoid_bias_threshold else True, + ), + epigraph_tolerance, + ) + solver.set_param("verbosity",1) + if sigmoid_bias < sigmoid_bias_threshold: + use_epsavg = False + else: + use_epsavg = True + + opt = intensity_optimization( + ref_intensity, + use_damping=True, + use_epsavg=use_epsavg, + sigmoid_bias=sigmoid_bias, + ) + + fig, ax = plt.subplots() + opt.plot2D(True,output_plane=view_2D_plane, ax=ax) + fig.savefig(f"results/optimization_problem_geometry.pdf") + + #the minimum linewidth and spacing constraint + #is activated only in the last epoch. This is done + #by adding the corresponding constraints to the same solver object. + #constraints are added as lambda functions, which + #are called directly by the nlopt solver. + if sigmoid_bias == sigmoid_biases[-1]: #if we are at the last sigmoid bias... + line_width_and_spacing = np.zeros(2) + grad_line_width_and_spacing = np.zeros((2, num_weights + 1)) + #compute the value of the costraints at this time + linewidth_constraint_val = line_width_and_spacing_constraint( + line_width_and_spacing, + epigraph_and_weights, + grad_line_width_and_spacing, + sigmoid_bias, + ) + solver.add_inequality_mconstraint( + lambda result_, epigraph_and_weights_, grad_: line_width_and_spacing_constraint( + result_, + epigraph_and_weights_, + grad_, + sigmoid_bias, + ), + tolerance_width_and_spacing, + ) + + #execute a SINGLE forward run before the start of each epoch and + #manually set the initial epigraph variable to slightly larger + #than the largest value of the objective function over the + #wavelengths and the lengthscale constraint + #making a call to opt like the following invokes + #the __call__ method of the opt instance and starts the optimization + print("Calibrating epigraph variable before starting optimization...") + epigraph_initial, empty_grads = opt( + [ + filter_and_project( + epigraph_and_weights[1:], + SIGMOID_THRESHOLD_INTRINSIC, + sigmoid_bias if sigmoid_bias < sigmoid_bias_threshold else 0, + ), + ], + need_gradient=False, + ) + #epigraph_initial contains the objective function values + #at each wavelength in the simulation. At each epoch, + #we set the initial epigraph variable taking the max of the objective functions + epigraph_and_weights[0] = np.max(epigraph_initial) + fraction_max_epigraph = 0.05 + if sigmoid_bias == sigmoid_biases[-1]:#if last epoch + epigraph_and_weights[0] = \ + (1 + fraction_max_epigraph) * max(epigraph_and_weights[0], + linewidth_constraint_val) + print( + f"epigraph-calibration:, {sigmoid_bias}, " + f"{str_from_list(epigraph_initial)}, {epigraph_and_weights[0]}" + ) + + epigraph_and_weights[:] = solver.optimize(epigraph_and_weights) + + optimal_design_weights = filter_and_project( + epigraph_and_weights[1:], + SIGMOID_THRESHOLD_INTRINSIC, + sigmoid_bias + ).reshape(NX_DESIGN_GRID, NY_DESIGN_GRID) + + fig, ax = plt.subplots() + ax.imshow( + optimal_design_weights, + cmap="binary", + interpolation="none", + ) + ax.set_axis_off() + if mp.am_master(): + fig.savefig( + f"optimal_design_beta{sigmoid_bias}.png", + dpi=150, + bbox_inches="tight", + ) + # Save the final (unmapped) design as a 2D array in CSV format + np.savetxt( + f"unmapped_design_weights_beta{sigmoid_bias}.csv", + epigraph_and_weights[1:].reshape(NX_DESIGN_GRID, NY_DESIGN_GRID), + fmt="%4.2f", + delimiter=",", + ) + + print("Normalization simulation completed") diff --git a/misc/logo.png b/misc/logo.png new file mode 100644 index 0000000000000000000000000000000000000000..6068339dde460bc769bbe4f66c5bddbf4c8b4716 GIT binary patch literal 7480 zcmds+_dlCo{P$JcPnDvys9mGlQqkI@C>5lN+O1iwRlCF%YE#s%8HyOG+FOj8MTr$# zP%8-%HDg7->G%Ev?w{}bhwGPIkH@*rd7qQ#`}KOAm{&R)^t2qbR8&;-AWb!WDymC& z{yj9*lwWAJgXXEIxX2(i6+_>gZM2`ck=Mnq-JQ1AfKq3`+ML~mU0n_H?z^uKXjcBK zu{8DpG7Xpa8M3BJ_7*M7biZ)nLa7>v=*dG@sD(bs}RF`ft zf1;wf{*3DX`PbjHoyb(*=w%0-JmgZCS|48zzmib&y&XUvJlWfP0b({!dD)k{HHK++ zZL<-l@4wF1DbRxU#cc7mqrEmq8`<_;t$?z+TqH@an<)nF}ML$W3 zuWKKm^!y2ukY>GU#iSY<_0=$AorB4}zBOE&uMej1Vt6b0ZX!^vu z5_gZyMt;Y?vGjXr1Akwdow2o&n~(E+?w*BO7g2N6yf@7@mRt{Q!Jl@>9 zye7_H-Co((W@ybeOT6$*AW%J2*{i$;Y&>MT_gR+Q!81bGdkt%5;S)SIn3Qn}>xe$B z#xi>Wi~sAJUtk6AIKMmw9QA3w#>2((>f?d#y@@VPC4u>W+}vYu)F zRhC@^ve;hAi;B7*sG?x;_PmH^q;*Gr#Q*2bL7xwNG}kXQ-dXk~&u}kU^b6#Ufeb{q zsqyv);)EB+08tvg&Yfi1+MMUKiFNmrE0~* zJG0fyh-P&(0McA1z$+@vxY?i=zei9kUNVmtSDucbM#KMFm|oVwe@a6c0%FnZx!??M z_4`s;cQUtgYu>GW>f~_WIey-iMlT;Pk(-hJWG9c;#$wQ3eC`Yr7&-UDD+k^Es=Gqn zMTjd2Y^rJ*8G!m;W;^{sBpEiq`d(TD!|9DxI|h`~^5GsJU#qG_I?e2@2=C!S!)Buy zpR#&8PKQvw&lyQCziA2_SQ$}r2>EPDAe~3K($yjO!(oNrIX8_VOi+J#VoQN5lwL-D zcx-9}NTx6GmdE%Togl9H028e=?ki68g0nH>dY2;|mf{ErH}z{5S|@*kj`kYrJjuv6 zcMzmAI=`3d0pjAlU4R+u$jHSaBF)WltL3<&+jiy0+ueN%jizVze&#vVt$Um@j%5xa z$sCGNbIu<7(<|&~-ai4bYw+)t^JjRI&W~dzB3dzQuRhQd#C>hTM-)tzEgYydVAo~? zHjcO3vS)V8{QiW<@}dw{(Y1LMM&WRDm>`?Y)jnW;!Ny9f)-NO{7bL*pHn$ zdMb@FDLz^Cr^jV zeS85^ojBE7vi#_7cFIzwYl$xJ;xhn}jF_;~12csWtioL=g3mz28_F6>OPgh7P}9Ln zV~r(9bHv}bviv0fNAhXe8%#`U?fQnXtuy;LFJ#bMHAsubf-y$oYhoUOhbMXjINVqA zV4M2B4&LS6+?;X@#~3b=pcIM};96pn;!GE!r<_!r+o)xUJvbPK;4_+>CV)wclXgK2C*QnoAqmIe%A555SHu{9U6WUiDloF~Z0pB_CPztA>S zbqLnXsfLB#u=TzB1vuR6hhgAW%BmT1T)xSU8w7|Jkocg3+TASA7;Rte1Z{EKC-67` z)+U8F9rc5eN~w^BrgQmp0*&AB(tyoWFz+V3YAm=0x@;1pY09JTVQr^nc+f-ylh(=R zQ>wvNVd`SN8IueW{_L@tW013*ZJ+l(eQd>L##V`ecN4UACGj>-|4jZVKeN|&dY+sC zmKO^P+@Z$2tQ_7NxoAT;u13AAAlmN;rVQ8y@4fF|-%~d87fv6LO<6Y|8fi+^Mcg?w z2~3{tcgc;Fk8SEG-Lu0JQVYvg3S3d8Q&=YLdf6T8w|l8cpt3w?2Nz9twqVnP*=r^u z(ZsQV^@-shVNg3oV8i06BQU1>$hL#0PBZw8Y;ZG?&dmHg!VvWU6p);7O5S?zU6ORP zWToqNc#N~tCCC32&C;NIzVaP$uo1KuG?hk&vIp~QOsVM9z8nz#-ueb{{tWBG=21$1 z-r_C2uq<0vl9Zusdf#)xnENcel^-z7Bi%XXyq3qkX_nj&0d`9cl}>}&P17dFO3NwP`V zn=EZD8K)y5s8TM+TPxBdnxhS8x0*kYPpTi1Nhf%CH-1{~%-b;gbDcf6i7f=*E@5?z z`%3{D@$hp_=y4lz&YbR#UODn)<#oY6pqc zB$Kokc~+&~A2%_SlX~@T&AF?s#Nie854RzFYYLxioQt5u8(IO1HJ5wWYqoZJ2!fT- z`znyL%3F1yGQH#@7*`$}0*Y?eRxC7=2x=j6zQznb?e1@B2q(7PNaQre(7RI?MP0^;cigDg5FSn5^QDSLXlN zrweaDH~thVA7J%1OEY9-`6}C?9|Sq+tU^E3$jsRTk>cpLcMmv9e+Jy@%sHNK z)+m<Gri>a0849@nA{r2FpU)o1JN;Y+Z!I|bt&f?7C6804Tae3=5fbiG2Dj@|%@_=E(k-)x@vlGBcC5T5`Zz2ij6QlG{m`N;c_f-pM}N z%0c<76E(7!D8TkmDHe}w;zwqhZz%PgJ|{+~r4SH;fXT%Qd-DSS%SgnXjAK*zNA=!Y z4^~Dlb;uOb3B!vQu9apb`M;TI-Fd3fpMXx5Tt-#2v=b&$S5=~7@-zpU;nwL?3cNma zoP4~4yvu~?%3n1|ik@N$MDV27yX^NY@XvkNbt zi|7*t#s?N8ZC^k8D0LF~2iM(~v}gs35k#)mjy!}bv%$XrC884M*H$W*j8lr~v`cq* zm)E);+0MM!`)nA%ruu4(w+93B@BMMFGFq!PGGS?1Qln5re*zBe@l~~Ku~&0{au}zn zOR$yp&;F&li3_z$YPy}5t$qPXa#Y_es+z^Fl*CZZ5efjyifl8M{dd~ZU4$z_?7LeqYwr7^eUFs;g_p3qYO z2zL}W{j1z3tI(DO-pd4fw+&t-=KfKc{vJ)t$KkZ48`0kbW5MIY>?wlDvDj$W2c@ND zCu>8y%iT#bnhC(1Y>uo)HC{x>`*&r^p%JQEQK?>V_doW$=H;+OZ?SAwY57rWJIqly zOWL$C8> z=sqD*h)!K@gvVOUF+MlT37R-p?whLst3XaG8}c5)_(93 zvg8`wXa9`I&mDCtTxx)DDt7(W-E@_@(hJN}%+t0=N(w4DWRL*5h=Cprnfe{P-}?+x z_wjTp>N55b)+chFp>ASYOKT(_0OSu)Dw~REu_NCN#rsRuV!}M+B0=P>y$N(-Xl0sS z?ijo&Gky+N-}CC45*w*;imYl}E4X3*Nb!#M(Ro4-%kD#3MB1(5BRv1c>r{%?XRmF) zDxCG>PHm@c1;ykf^5}T{VXR;Z0kZV4S-TH>tXdb63GQZ~Rx0-Fqj=SY!fA66(PHS> z69T2sg|T^{Ur589fE42e>Rxn*!V|lVsST(=d{)|%Xf(jw&!?nF`J)o_OgFRT? z;;U)14!e>LYSlrCLDli=%HVj>wdc45VBo7fwooA*7X752lMqyj1mkwHV>A(f#fjD$ zBkRBpC0hR=wJuHTb(F`BPM|nn=cNn9s5h#))P08b;NKl_e3hDx;GjcoMmjT#3K|1M z$K8*XhI$-#*tmnT0@aBgm6$2H(>Ypyep}U`7Bd_7C%DFn?_yJJLmI)~0=4>o$;C$+ zAg+N1N6h`*r)T%&>R3K1&2|CacBLwXegR%T!Q~>RQ%riAMn<1-LP@WfR8I}ey||%H z$F!#+pPEuK*IA1?LZY=k2mFf}p$Q27HRn=X1Ub7u$J`!U+=;4y&qZd3) z#iJo`0!@p$V&O!JLWev*$|(kEdFdeCT}{9I1@bpI!E%zCMO%+D67)Wde=-gbv_7U= zDk3hmdHSV3J1Rm{*qXVk~@-L(GPo%k<#%z^=A6pOkvfl#Ag0~{kjioReRtE;SxU8 zYT;zIKX>%QWeGgVkFk2(Il}7Z_-g553vd4=NN*^PVCb4PO*pKZ*D{k0kQd)5d=WN4Tr{ZsvgEyZR-MBvgI>{P8X9}O$=m#;xIOg zuKn~WHF~WzgiU?^!`n}kwdQiEUh;~kz9N5&PoW2L4%^%&?Gc6H7eLRz?u7$*l9Kl% zDF4P$lOBpsSG}UBT3k?Jr#m=i^h>Fzq7opd*^oQ~+S?!^Q003UmGjOq! z`$xzzh~!cQ(z+qqQC&$eem=@JW-pfLHZG3+*5R65MRA`V6QK!^mzGR!-~bmG0_ir+ zteTymnRZz3M)uUWb0{P_S^KE~a>0iMDB{$lNfe%ld|DWGX+((iutRT}a*j`aPM(bN zi4qD@X~SO>DUqaQhZ=Zf)%t%!`*jX4u!^Vb+<8w?u3-{(;4TMy%>iZ7W8&&zris{i zJ3U4xOjHIjwDkdHbTrdJr1;kb$%N8xoJOe#Jvn=ijq|K^_4H<;jIo$P7m6*qwJ(A0 zx?HJeH8n6}6Z^z0bw9#+T)47)(R0K1@7Vu=3!vU(*AMceTB$;v1vIZJT=B@t<8JD| zfR<1`2fJ6mk37zm_(ftM*wEAW%oD0Xv7e2#@tdxb*MgXK(IcYO_#kLuWhs64|5m1t zqTcS6m}vOEO<|=MJyA7&*IHU2c^r!!4G&v0p*jqr@Oic&lraA96bOi8)unbN1Z>3J z_QxobtM&u_6v(7_{ZtSa`U!Fec$*Uv$)3(zt;e(6)N>BH`mX5 zb4UKu$#8stMV0<2b|hxL?_HO1=y>$?_Yyb$34VAt15igr&j2s}WB0R@g3^ITG$+o- ze`KIB;oJQ~iN0LteH*dDikE=*=FEp#1I;~ zw1iPVZ;XBu?FQ|B#K|PlvgAHlePHN1m$t_oB333Gt(ZT(G3%;mMVQ}t7zXfI3FbPE zSC346o!IP4KH33nh>z!4I3B#3si@AIorfs|3H<4ld_-|8`YKwc;v=J z#y?IPbJ3h#r~ zPnn%H6k7^^gK)yfdKdh|R&y--&<;y9bxEeWc-d9ZNhTowckFlJ(h@pK_M`BW&Bt{* zI7RBCt zMdt@$OZPb^D$ie@Nt*F=Mf(}K=#jrev4 z^&Sso0H~nG3^)p0=7^iY`{}CcAMwtdZ9K&zS4~CIXtTfI=VF5)pU4jDWgR=3@h2u? zc2)@YvfYNC1^{G+xsZ&*^q<)Z#XS5}txNX(o;K;9kkUh!+f-SqqJ2wLSq%TbUc&uf iZWouym#L`8WWA)BxYjOfE6R;16$q%KR-yVf^nU=-b+DuW literal 0 HcmV?d00001 diff --git a/misc/peg.png b/misc/peg.png new file mode 100644 index 0000000000000000000000000000000000000000..01830ee6f29950791919e0c23e5fd21a0ec24f8c GIT binary patch literal 4250 zcmV;L5M}R)P)1^@s6(C8?-000nSNkl(^j|5U{Yb#osn{e;mJs21mKq4`NsYEj2ukarNq=A6}6c-ktrKN=hh=HLYyuE#( zJ89B{fi>YEJs2EB_RTDGc6RcDpsA^ez@T7E($kv|pe7uo!h(D>H#Mo;JUvN4!6Bf4 z{n9j{Ahi)4l$V<$c+V1|pM3oT1oJ;G#t8)}^+pPM@7xi!7gsko*x64OB{Np$2?MF2 zz8)pTMWWZQuC9*IuyE+;>WY>ft-yqVbR#(leSLkSTDOITCA@rmMkPBMp$P$LkXW?( zx;j)>RieMYe>Cexa7ZZhi7)j_&@TqkU2AOE!*H;}}fqK9Bm zs}~#^#ss_BS(!pGVt67VXX1#&1q1i{v8ts6nns3#WrZ+~8%XJ?H_+45BV-K@2g^egZa|u(U!!hUb4}#=`GCekssP*b3roOAI3r0pp7<@2D z%zh5&wG%^n`uYeVr-{mCkZX5GKx}7hYy!*bTAYc0o|)?DlLukFXr7p9mGV9okh)|t z6y)cVFr!yK*;GeI(D31ARu=K9h3h3cfdBMBbS0-l5*5zqCm}?y+#!7j#)CJtvBdKF z0wH38c*YD8Rm^E5n&(o`4}r?c%wRc^%2r3cX`2TcXtHr;GxFgK^p+0A1}ehD~B?iMYD4>!wkIM5Fo?k`oe2 zbOTyi+Km1ZYVl|Q`}4dN=xu4m!t`@QiKi)@&e_EU<nICcryxAFNJ#fV84aCd}uW_AfYOAOu?@pkVaZCsxN%~{E+q*cK#pKqeIcqP|u5RNLU!YxAsI& z*Q3mXXZ#{W9Nnwxos8bPQ6HtEoNY?Y%*^0A&7Gt|ZV2{)s4Eao^i0lsJ0MxToP-{# zO*Flt1#kWVWqS|9&~XYDJJeIWX*h?J;>l0QkX<-Gno1aALtm*I(So@Xh+>3v|0p zg;d(cDMwX-^qmipci>xv=d*YI363wV7rIOZKBXY33Nq=*pcbegFBffX+>0K$DjQon zsHbM5@!AbIZQZ0YZA)W0SAAnfFHw^CdLs7HcTC=iiwsHJXr5P!^XvaZ;qhdDPt3o0 zU{MXgqK zRTWbr5<=8xx)2)&6AxFoy}t#eCB;g`e29@q`}dbnckVL7=e8~ZF7IsSm7}TP@625Z zh>BM^`h%@-{m)l~EK%5XD=jG@wXs@31JcqC^-<_)5|UA*AxF#K4kMc=@RgOB$g8Qs zT_w7STE*7ghls-em{|`!4J+&Nke|cA1oC^Muk~tY)zw{O^tawg7_KeT}jJv@; z=;&aFDvnN0sD6RyyYp8S)}lqO+}eAj&_JCh5yS}+EkDy!A(W^dh>VKHV9y<773QmK zf}lZfB)Otr>Q8~D-XtW{=Ocx<5_~#Mm(=||t!#0`B}Bd?yLmD2{o6;Jir>=83IqF( zBK3biV{Ku<;a#v_87FjW3?C;*&7>u-kR;K>kcX!i`fg-FKQaXQ0@`VULPGHKtT?o$ zXCvmT&(P=V2HFnDN7Sj@Z{&!9J41;*H4SyFY%D}ZZZ=EAkT|<`P3<^I@=GJh{ps0= z{PD1ut#N}SHZG$xT3g#dBQqZsfqp3Il&Y+?xe-USYWofp?c57%l7S34`i0ULGd(Z} zSC_0o#W0nMW@STR+k)8k0bnF8cbz@VQOWD<&DY>#YiE}laxzB2B|>5=iz{~ zz!>1{hSWE=!jo8$q5|a(rnQX?)Nf?t>iXX@Lp(#`tDIYvyx#jXcXw25`w+!{-7nvR z%bFE%+5R?D)Er5XQBnBXz=7yh#ndrXN9&1RQWMcJ#|ctuN;1>fAeKT%-_S|_SS-uD zfXuucrL8rzBg93|NbV$t?RXa&b0c|!X8Xxg(EG!Aq;C2j`Lc8^%@IeyrRkTEM0+gH zItQV%m95wp9Qqvg3CjgWimEd#00!e5%!Z*WA>i_84BXZXV(2NZunMnQy zE!sNM^4WgmfAl#socHT5!Rpm#d9NuAroVLwW^zjpi3YsU7mqWA(BWigXas{x$+))p zUDi^mE4hB432MY&6GAMIXiv4O%7n6kX=Q{eZxc^91n&8mY2M-^n$lGxX^@l+Z!;`r z(|r+h{wQyRA3%zhm*+o&mb78vh`Jiq_TEO)b;WB%t*v9L45#DPaZbN1976UVcVIXwQpcswyLiW-0r=KaaG{e^7WtEu6v9`K(k84SkN8xS)T4-L&C2GX%4* zo*LbCD=97}o2nBe%*8|_eA8Ou&n#ID&sWyN;jIlyuP@F0d}HZ>{_deR9rqU!;j(oj zii_Bd@KlFsXXMMh^oR3%z4mMPy!$qjNs!4NeaJ&o#V76k(D(OvZ}@sGT^z8aEJM&f zStp711ktqZXo)srIce+RMBp6Q&Yp?L<6ko0kMMQ8O%DpfX}?%>R}Wc1`}wh$bMg@L zqG-S?BMn^-{U%yslGL^^FJL%1EM0)PC%$5Sq8W@<@<|wNRHfJRtaBjoy4h7@BQp~` z(<-c0QE85wR$8etq#+9ZqrF11&JJ{U52aDLh}OP>NKO7Asmbe4u|vYGY$g#$W5G@$ zOWTo$D_CEDe_Wlv60I3SCfwB99n13`6%#2yv_{w7E@kpt#%(hEIAf8*)Rqy!Sdw&_ zb3NZV_@S+aecg3+wXwQa*o_fXcN;l71&_8K_mC^n)6>oTJ3RRnMba93w5%Dib-lIF zASvGAL`&o?EzP7apc7Qjq^|l6YR_I|ShP)ykl4SvanwCrVLan$&csM?d-tb2 z!_ZK~9#2z44bS%Tcf$;yU{oM+y)xQ7MK|6|s=D4JOGy*=_c~-qsL3YPS39M*_wt(G zp!~owInQ%S9KF8QHXImMp`}7P<`S$Sy^>sXE_gO-GKPAVM$+>S*k5XL@jW(>Mhv;>-%ma)F#7_K{mK zQh2wT1#7%@#Ov1J_`EA=*=iJ>_=$7=MDjAUd0xi2T(^;{I2pB=DL^hba~SpuW-A4$ z;M3jm(1bF^7hES{%`{$_Jvz1+yr-BK-uuW2iVj>MED zlKt`~WbE9>`dv#~3v2HR>kk+=kglxTfS%fVEIiNs7Kg98>iCbi5WgO3TAGOXm|~fgWOy?1t4-b5$)%{(Rg(61;*@FeYz!i!laMmp-L3p(trPec)5fcSL(lpu-0wIe+fxx{h?eEc2^-Cf0qGQ9U?UH|$ZiAIG(;9veEJpt)) wARu6jCnxC%NCFTLz@jFT$)5aATq=P2e<*umMK;8f&c&j literal 0 HcmV?d00001