From 36c1db0df21e4bafcf09d976867b154dec651c21 Mon Sep 17 00:00:00 2001 From: Edoardo Vergerio Date: Tue, 28 Jul 2026 22:41:42 +0200 Subject: [PATCH] code alignment between the two versions. There are now two possibilities: 1)optimize the average objective function, this generally leads to different performance of the converted modes; 2) optimize using an epigraph formulation. This assures that whatever the final epigraph value, the performance of the converted modes will be greater (in terms of mode overlap) than the epigraph itself. The epigraph is implemented as a constrained optimization, which is stricter than optimizing the average --- metaprop_fibermode_doublet_epigraph.py | 159 ++++++++++++++++--------- metaprop_fibermode_doublet_multi.py | 11 +- 2 files changed, 109 insertions(+), 61 deletions(-) diff --git a/metaprop_fibermode_doublet_epigraph.py b/metaprop_fibermode_doublet_epigraph.py index 724cf8c..06bc55a 100644 --- a/metaprop_fibermode_doublet_epigraph.py +++ b/metaprop_fibermode_doublet_epigraph.py @@ -12,6 +12,7 @@ import nlopt import numpy as np import matplotlib.pyplot as plt +import matplotlib.patches as patches #set the global font size for all plots plt.rcParams.update({'font.size': 14}) @@ -101,6 +102,7 @@ P_2_dagger_nat = np.fft.ifftshift(P_2).T.conj() #Plot zoom parameters [um] zoom_x = 60 zoom_y = 60 +zoom_core_diameter = 2 * 6 full_view_x = size_x * 1e6 / 2 full_view_y = size_y * 1e6 / 2 ############################## @@ -121,14 +123,23 @@ def get_pattern() -> np.ndarray: return pattern.flatten() -def get_fiber_mode_pattern(n_mode_sel): - '''Returns the target fiber mode pattern''' - print("Computing target fiber mode pattern for LP %d1" % n_mode_sel) +def get_fiber_mode_pattern(mode_list:list)->list: + '''Returns the fiber mode pattern for the selected n index. + Args: + mode_list: the list of the selected fiber mode indices + Returns: + The fiber mode patterns corresponding to the selected indices returned + as a list of 2D arrays with the same order as the input mode_list + ''' + #check if the modes in mode_list have already been compute by checking + #the saved files in the 'results' folder. If they have, load them from the files. If not, compute them and save them to the files. + + print("Computing fiber patterns for mode indices: ", mode_list) lda = 1.55e-6 - n_co = 1.4630 - n_cl = 1.4585 - a = 6e-6 + n_co = 1.4630 #core refractive index + n_cl = 1.4585 #cladding refractive index + a = 6e-6 #fiber core radius #find the v parameter range that corresponds to the 1550 nm wavelength #v_1550 = 2 * np.pi * a / lda * np.sqrt(n_co**2 - n_cl**2) @@ -138,15 +149,13 @@ def get_fiber_mode_pattern(n_mode_sel): V_list = np.linspace(2, 5.1, 20) V_obs = v_1550 #b, V = np.meshgrid(b_list, V_list) - n_list = [0,1,2] b_solutions_at_V_obs = [] - for n in n_list: + for n in mode_list: b_solutions = [] for V in V_list: LHS = np.sqrt(1-b) * jv(n+1, V * np.sqrt(1-b)) / jv(n, V * np.sqrt(1-b)) RHS = np.sqrt(b) * kv(n+1, V * np.sqrt(b)) / kv(n, V * np.sqrt(b)) - intersection_index = np.argmax(np.where(np.isclose(LHS, RHS, atol=0.000001),1,0).flatten()) @@ -156,29 +165,32 @@ def get_fiber_mode_pattern(n_mode_sel): RHS = np.sqrt(b) * kv(n+1, V_obs * np.sqrt(b)) / kv(n, V_obs * np.sqrt(b)) b_solutions_at_V_obs.append(b[np.argmax(np.where(np.isclose(LHS, RHS, atol=0.00001),1,0).flatten())]) - b_mode_sel = b_solutions_at_V_obs[n_mode_sel] + output_patterns = [] + for enum_index,n_mode_sel in enumerate(mode_list): + b_mode_sel = b_solutions_at_V_obs[enum_index] - k_co = 2 * np.pi / lda * n_co - k_cl = 2 * np.pi / lda * n_cl + k_co = 2 * np.pi / lda * n_co + k_cl = 2 * np.pi / lda * n_cl - beta = np.sqrt(b_mode_sel * (k_co**2 - k_cl**2) + k_cl**2) + beta = np.sqrt(b_mode_sel * (k_co**2 - k_cl**2) + k_cl**2) - chi_co = v_1550 * np.sqrt(1 - b_mode_sel) / a - chi_cl = v_1550 * np.sqrt(b_mode_sel) / a + chi_co = v_1550 * np.sqrt(1 - b_mode_sel) / a + chi_cl = v_1550 * np.sqrt(b_mode_sel) / a - E_field = np.zeros_like(X, dtype=complex) + E_field = np.zeros_like(X, dtype=complex) - for i in range(Nx): - for j in range(Ny): - r = np.sqrt(X[i,j]**2 + Y[i,j]**2) - if r < 10*a: - phi = np.arctan2(Y[i,j], X[i,j]) - if r < a:#inside core - E_field[i,j] = jv(n_mode_sel,chi_co * r) / jv(n_mode_sel, chi_co * a) * np.cos(n_mode_sel * phi) - else:#outside core - E_field[i,j] = kv(n_mode_sel, chi_cl * r) / kv(n_mode_sel, chi_cl * a) * np.cos(n_mode_sel * phi) + for i in range(Nx): + for j in range(Ny): + r = np.sqrt(X[i,j]**2 + Y[i,j]**2) + if r < 10*a: + phi = np.arctan2(Y[i,j], X[i,j]) + if r < a:#inside core + E_field[i,j] = jv(n_mode_sel,chi_co * r) / jv(n_mode_sel, chi_co * a) * np.cos(n_mode_sel * phi) + else:#outside core + E_field[i,j] = kv(n_mode_sel, chi_cl * r) / kv(n_mode_sel, chi_cl * a) * np.cos(n_mode_sel * phi) - return E_field + output_patterns.append(E_field) + return output_patterns def get_source_field(beam_waist,center_x=0,center_y=0,): '''Returns a normalized Gaussian source field for a given beam waist. Normalization consists in scaling @@ -211,7 +223,7 @@ def propagate_source(source_to_propagate): propagated_source_field = ifft_operator.output_array.flatten() return propagated_source_field -def make_1Dplot_of(x_axis,y_axis,plot_zoom_x=zoom_x,save_name="plot_1D",) -> None: +def make_1Dplot_of(x_axis,y_axis,legend_labels,plot_zoom_x=zoom_x,save_name="plot_1D",) -> None: '''Makes a 1D plot of the given data and saves it to a PDF file. Args: x_axis: list of x-axis data (list of 1D arrays) @@ -225,8 +237,12 @@ def make_1Dplot_of(x_axis,y_axis,plot_zoom_x=zoom_x,save_name="plot_1D",) -> Non raise ValueError("Lists x_axis and y_axis must have the same length.") plt.figure() for i in range(len(x_axis)): - plt.plot(x_axis[i], y_axis[i]) - plt.xlabel("x [um]") + if len(legend_labels) != len(x_axis): + raise ValueError("Length of legend_labels must match the number of plots.") + else: + plt.plot(x_axis[i], y_axis[i], label=legend_labels[i]) + + plt.legend() plt.tight_layout() plt.savefig("results/" + save_name + ".pdf") plt.xlim(-plot_zoom_x, plot_zoom_x) @@ -264,6 +280,37 @@ def make_2Dplot_of(given_field,choose_quantity="amplitude",save_name="plot",plot plt.savefig("results/" + save_name + ".pdf") plt.close() +def make_polarPlot_of(given_field,choose_quantity="amplitude",save_name="polar_plot",plot_zoom_r=zoom_x) -> None: + '''Makes a polar plot of the given 2D field in a zoomed region. + Saves the plot with the chosen filename to a PDF file. + Args: + given_field: the 2D field to be plotted (flattened) + choose_quantity: either "amplitude" or "phase" to select quantity to plot + save_name: name of the file to save the plot (without extension) + plot_zoom_r: zoom level for the radial axis (in micrometers) + Returns: + None + ''' + _, ax = plt.subplots(figsize=(6, 6)) + if choose_quantity == "amplitude": + given_field_2d = np.abs(given_field.reshape((Nx, Ny))) #assure the field is 2D before plotting + im = ax.imshow(given_field_2d,origin='lower') + elif choose_quantity == "phase": + given_field_2d = (np.angle(given_field.reshape((Nx, Ny))) + 2*np.pi) % (2*np.pi)#assure the field is 2D before plotting + im = ax.imshow(given_field_2d, cmap='hsv', vmin=0, vmax=2*np.pi) + else: + raise ValueError("choose_quantity must be either 'amplitude' or 'phase'") + + center_x, center_y = Nx//2, Ny//2 + radius = plot_zoom_r * 1e-6 / (size_x/2) * (Nx/2) #convert zoom in um to pixels + circle = patches.Circle((center_x, center_y), radius, transform=ax.transData) + im.set_clip_path(circle) + ax.set_xlim(center_x - radius, center_x + radius) + ax.set_ylim(center_y - radius, center_y + radius) + ax.axis('off') + cbar = plt.colorbar(im, ax=ax) + plt.savefig("results/" + save_name + ".pdf") + def shift_spatial_grid(shift_amount_x,shift_amount_y): '''Shifts the spatial grid by a given amount. This is useful for centering the source field on the metasurface. Args: @@ -337,8 +384,7 @@ input_field_list = [input_field_1, input_field_2] ########################### #initialize target field target_mode_indices = [1,2] -fiber_mode_pattern_1 = get_fiber_mode_pattern(target_mode_indices[0]) -fiber_mode_pattern_2 = get_fiber_mode_pattern(target_mode_indices[1]) +fiber_mode_pattern_1, fiber_mode_pattern_2 = get_fiber_mode_pattern(target_mode_indices) target_Efield_1_2d = fiber_mode_pattern_1 / np.sqrt(np.sum(np.abs(fiber_mode_pattern_1)**2)) target_Efield_1 = target_Efield_1_2d.flatten() #normalized target field (intensity = 1) @@ -458,7 +504,8 @@ def value_and_grad(weights): return np.array(values), np.array(gradients) - +obj_values1 = [] +obj_values2 = [] def epigraph_constraint(result,epigraph_and_weights,gradient)->None: '''Constraint function for the epigraph formulation Args: @@ -472,6 +519,8 @@ def epigraph_constraint(result,epigraph_and_weights,gradient)->None: weights = epigraph_and_weights[1:] obj_val, grad = value_and_grad(weights) + obj_values1.append(obj_val[0]) + obj_values2.append(obj_val[1]) if gradient.size > 0: gradient[:,0] = -1 @@ -498,23 +547,12 @@ def obj_fun(epigraph_and_weights, grad): grad[1:] = 0 opt_history.append(epigraph) - - plt.figure() - plt.plot(opt_history) - plt.xlabel('Iteration') - plt.ylabel('Epigraph Value') - plt.title('Optimization History') - plt.tight_layout() - plt.savefig("results/modeconv1_optimization_history.pdf") - plt.close() + #make a 1d plot of the optimization history together with the objective values + make_1Dplot_of(x_axis=[np.arange(len(obj_values1)), np.arange(len(obj_values2)), np.arange(len(opt_history))], y_axis=[obj_values1, obj_values2, opt_history], legend_labels=["ob1","ob2","t"], save_name="optimization_history", plot_zoom_x=len(opt_history)) return epigraph -if __name__ == "__main__": - ########################### - #initialize the normalized phase masks and global phase mask - w_norm = np.zeros((2*S)) #concatenated normalized parameters for both masks - +def initialize_masks(): X1,Y1 = shift_spatial_grid(-source_spacing/2,0) rho1 = np.sqrt(X1**2 + Y1**2) circular_mask_1 = np.where(rho1.flatten() < 10*6e-6, 1, 0) @@ -525,7 +563,7 @@ if __name__ == "__main__": target_pattern_1_norm_phase = (np.angle(target_Efield_1) + 2 * np.pi) % (2 * np.pi) / (2*np.pi) target_pattern_1_norm_phase = np.roll(target_pattern_1_norm_phase, shift=int(source_spacing/sampling_period/2), axis=0) - + target_pattern_2_norm_phase = (np.angle(target_Efield_2) + 2 * np.pi) % (2 * np.pi) / (2*np.pi) target_pattern_2_norm_phase = np.roll(target_pattern_2_norm_phase, shift=-int(source_spacing/sampling_period/2), axis=0) #on metasurface 1, set initial guess on normalized phase elements @@ -533,7 +571,13 @@ if __name__ == "__main__": #on metasurface 2, generate random initial guess for normalized phase elements #w_norm[S:] += np.random.rand(S) * circular_mask_1 + np.random.rand(S) * circular_mask_2 +if __name__ == "__main__": + ########################### + #initialize the normalized phase masks and global phase mask + w_norm = np.zeros((2*S)) #concatenated normalized parameters for both masks + initialize_masks() phase_mask = phase_given_w(w_norm) + #make an initial plot of the phase masks make_2Dplot_of(phase_mask[:S].reshape((Nx, Ny)), choose_quantity="phase_mask", save_name="initial_phase_mask_1", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) make_2Dplot_of(phase_mask[S:].reshape((Nx, Ny)), choose_quantity="phase_mask", save_name="initial_phase_mask_2", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) @@ -570,7 +614,7 @@ if __name__ == "__main__": solver.set_maxeval(opt_max_eval) solver.set_param("dual_ftol_rel", 1e-5) solver.set_param("verbosity", 1) - solver.set_param("inner_maxeval", 5) + solver.set_param("inner_maxeval", 10) print("Starting optimization...") start = True @@ -585,12 +629,12 @@ if __name__ == "__main__": all_source_on_field = input_field_list[0] + input_field_list[1] verify_out_field_all = forward_propagate(all_source_on_field, opt_weights)[-1] - make_2Dplot_of(verify_out_field_1, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_source1_amplitude", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) - make_2Dplot_of(verify_out_field_1, choose_quantity="phase", save_name="modeconv1_optimized_output_field_source1_phase", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) - make_2Dplot_of(verify_out_field_2, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_source2_amplitude", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) - make_2Dplot_of(verify_out_field_2, choose_quantity="phase", save_name="modeconv1_optimized_output_field_source2_phase", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) - make_2Dplot_of(verify_out_field_all, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_allSource_amplitude", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) - make_2Dplot_of(verify_out_field_all, choose_quantity="phase", save_name="modeconv1_optimized_output_field_allSource_phase", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) + make_polarPlot_of(verify_out_field_1, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_source1_amplitude", plot_zoom_r=zoom_core_diameter) + make_polarPlot_of(verify_out_field_1, choose_quantity="phase", save_name="modeconv1_optimized_output_field_source1_phase", plot_zoom_r=zoom_core_diameter) + make_polarPlot_of(verify_out_field_2, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_source2_amplitude", plot_zoom_r=zoom_core_diameter) + make_polarPlot_of(verify_out_field_2, choose_quantity="phase", save_name="modeconv1_optimized_output_field_source2_phase", plot_zoom_r=zoom_core_diameter) + make_polarPlot_of(verify_out_field_all, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_allSource_amplitude", plot_zoom_r=zoom_core_diameter) + make_polarPlot_of(verify_out_field_all, choose_quantity="phase", save_name="modeconv1_optimized_output_field_allSource_phase", plot_zoom_r=zoom_core_diameter) #save the output fields for later use saveFields = False @@ -602,13 +646,12 @@ if __name__ == "__main__": # Reshape fields to 2D for visualization phase_mask_1 = phase_given_w(epigraph_and_weights[1:S+1]).reshape((Nx, Ny)) phase_mask_2 = phase_given_w(epigraph_and_weights[S+1:]).reshape((Nx, Ny)) - phase_mask_zoom = 150 # Plot: Phase masks - make_2Dplot_of(phase_mask_1, choose_quantity="phase_mask", save_name="phase_mask_1", plot_zoom_x=phase_mask_zoom, plot_zoom_y=phase_mask_zoom) - make_2Dplot_of(phase_mask_2, choose_quantity="phase_mask", save_name="phase_mask_2", plot_zoom_x=phase_mask_zoom, plot_zoom_y=phase_mask_zoom) + make_2Dplot_of(phase_mask_1, choose_quantity="phase_mask", save_name="phase_mask_1", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) + make_2Dplot_of(phase_mask_2, choose_quantity="phase_mask", save_name="phase_mask_2", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) #save the phase masks for later use - saveMasks = True + saveMasks = False if saveMasks: np.savetxt("results/optimized_phase_mask_1.txt", phase_mask_1) np.savetxt("results/optimized_phase_mask_2.txt", phase_mask_2) diff --git a/metaprop_fibermode_doublet_multi.py b/metaprop_fibermode_doublet_multi.py index 1232b5f..8ca4730 100644 --- a/metaprop_fibermode_doublet_multi.py +++ b/metaprop_fibermode_doublet_multi.py @@ -52,7 +52,7 @@ rho = np.sqrt(X**2 + Y**2) sampling_period = xs[1] - xs[0] d0 = 345e-6 #propagation distance: source - MS1 -d1 = 500e-6 #propagation distance: MS1 - MS2 +d1 = 100e-6 #propagation distance: MS1 - MS2 d2 = 300e-6 #propagation distance: MS2 - target d = [d1, d2] #d = [d1,d2] d1:distance MS1-MS2, d2: distance MS2-target @@ -510,6 +510,8 @@ def value_and_grad(weights): return np.array(values), np.array(gradients) opt_history = [] +opt_history1 = [] +opt_history2 = [] def obj_fun(weights, grad): ''' Objective function for the optimization. Computes the overall cost function by taking the mean value @@ -525,6 +527,7 @@ def obj_fun(weights, grad): obj_val, grads = value_and_grad(weights) + average_obj_val = np.mean(obj_val) if grads.size > 0: @@ -540,7 +543,9 @@ def obj_fun(weights, grad): make_2Dplot_of(np.abs(grad[S:].reshape((Nx, Ny))), choose_quantity="amplitude", save_name="gradient_magnitude", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) opt_history.append(average_obj_val) - make_1Dplot_of([np.arange(len(opt_history))], [opt_history], plot_zoom_x=len(opt_history), save_name="optimization_history") + opt_history1.append(obj_val[0]) + opt_history2.append(obj_val[1]) + make_1Dplot_of([np.arange(len(opt_history))]*3, [opt_history,opt_history1,opt_history2], plot_zoom_x=len(opt_history), save_name="optimization_history") return average_obj_val @@ -619,7 +624,7 @@ if __name__ == "__main__": make_polarPlot_of(verify_out_field_all, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_allSource_amplitude", plot_zoom_r=zoom_core_diameter) make_polarPlot_of(verify_out_field_all, choose_quantity="phase", save_name="modeconv1_optimized_output_field_allSource_phase", plot_zoom_r=zoom_core_diameter) - #make 1D slices of target/outpu fields to compare + #make 1D slices of target/output fields to compare slice_y_index = Ny//2 slice_x = xs * 1e6 slice_target_1 = target_Efield_1.reshape((Nx, Ny))[slice_y_index,:]