From ebcaeb56433057376284fd24d06f43e2fc20c421 Mon Sep 17 00:00:00 2001 From: Edoardo Vergerio Date: Tue, 28 Jul 2026 13:18:22 +0200 Subject: [PATCH] fixed the epigraph formulation, which consisted in correcting the sign of the gradients to comply with the definition of epigraph minmax problem. See notes on the epigraph formulation for further explanations. --- metaprop_fibermode_doublet_epigraph.py | 74 ++++++++++++++------------ 1 file changed, 40 insertions(+), 34 deletions(-) diff --git a/metaprop_fibermode_doublet_epigraph.py b/metaprop_fibermode_doublet_epigraph.py index 10ed1c6..724cf8c 100644 --- a/metaprop_fibermode_doublet_epigraph.py +++ b/metaprop_fibermode_doublet_epigraph.py @@ -43,7 +43,7 @@ sampling_period = xs[1] - xs[0] d0 = 345e-6 #propagation distance: source - MS1 d1 = 500e-6 #propagation distance: MS1 - MS2 -d2 = 500e-6 #propagation distance: MS2 - target +d2 = 300e-6 #propagation distance: MS2 - target d = [d1, d2] #d = [d1,d2] d1:distance MS1-MS2, d2: distance MS2-target phase_min = 0 @@ -53,12 +53,13 @@ phase_max = 2 * np.pi ############################## #optimization algorithm epigraph_tolerance = np.array([1e-4]*2) +#epigraph_tolerance = np.array([1e-4]) ############################## ############################## #Source beam_waist = 9e-6 -source_spacing = 150e-6 +source_spacing = 300e-6 ############################## ############################## @@ -201,7 +202,7 @@ def get_source_field(beam_waist,center_x=0,center_y=0,): return gaussian_field_normalized def propagate_source(source_to_propagate): - '''Returns the propagated source field just before metasurface 1 (MS1). This involves propagating the + '''Returns the propagated source field (flattened) just before metasurface 1 (MS1). This involves propagating the source for a distance d0: from its plane of definition up to MS1.''' fft_input_array[:,:] = source_to_propagate.reshape((Nx, Ny)) fft_operator() @@ -330,6 +331,7 @@ make_2Dplot_of(input_field_2_2d, choose_quantity="phase", save_name="input_field #list of input fields to be propagated within the optimization loop input_field_list = [input_field_1, input_field_2] +#input_field_list = [input_field_1] ########################### ########################### @@ -351,6 +353,7 @@ make_2Dplot_of(target_Efield_2_2d, choose_quantity="amplitude", save_name="targe make_2Dplot_of(target_Efield_2_2d, choose_quantity="phase", save_name="target_field_2_phase", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) target_field_list = [target_Efield_1, target_Efield_2] +#target_field_list = [target_Efield_1] ########################### def forward_propagate(input_field:np.ndarray,weights:np.ndarray) -> np.ndarray: @@ -361,37 +364,37 @@ def forward_propagate(input_field:np.ndarray,weights:np.ndarray) -> np.ndarray: print("Starting propagating...") intermediate_fields = [None, None] start_time = time.time() - #phase_mask = phase_given_w(weights) - fft_input_array[:,:] = (input_field * np.exp(1j * phase_mask[:S])).reshape((Nx, Ny)) + phase_mask = phase_given_w(weights) + fft_input_array[:,:] = np.fft.ifftshift((input_field * np.exp(1j * phase_mask[:S])).reshape((Nx, Ny))) fft_operator() - - print("BEFORE") + #print("BEFORE") fft_power = np.sum(np.abs(fft_operator.output_array)**2)/S non_propagating_mask = np.where((wavelength * nu_parallel)**2 > 1, 0, 1) non_propagating_mask = np.fft.ifftshift(non_propagating_mask) propagating_power = np.sum(np.abs(fft_operator.output_array * non_propagating_mask)**2)/S - print("FFT power: %.4e, Propagating power: %.4e, Non-propagating power: %.4e" % (fft_power, propagating_power, fft_power - propagating_power)) - + #print("FFT power: %.4e, Propagating power: %.4e, Non-propagating power: %.4e" % (fft_power, propagating_power, fft_power - propagating_power)) + #FIRST PROPAGATION ifft_input_array[:,:] = fft_operator.output_array * P_1_nat #element wise product - print("AFTER") + #print("AFTER") fft_power = np.sum(np.abs(ifft_operator.input_array)**2)/S non_propagating_mask = np.where((wavelength * nu_parallel)**2 > 1, 0, 1) non_propagating_mask = np.fft.ifftshift(non_propagating_mask) propagating_power = np.sum(np.abs(ifft_operator.input_array * non_propagating_mask)**2)/S - print("FFT power: %.4e, Propagating power: %.4e, Non-propagating power: %.4e" % (fft_power, propagating_power, fft_power - propagating_power)) + #print("FFT power: %.4e, Propagating power: %.4e, Non-propagating power: %.4e" % (fft_power, propagating_power, fft_power - propagating_power)) ifft_operator() - intermediate_fields[0] = ifft_operator.output_array.copy() - print("Field intensity after first phase mask and propagation: %.4e" % np.sum(np.abs(intermediate_fields[0])**2)) + intermediate_fields[0] = np.fft.fftshift(ifft_operator.output_array.copy()) + #print("Field intensity after first phase mask and propagation: %.4e" % np.sum(np.abs(intermediate_fields[0])**2)) - fft_input_array[:,:] = (ifft_operator.output_array * np.exp(1j * phase_mask[S:].reshape((Nx, Ny)))) + fft_input_array[:,:] = np.fft.ifftshift((np.fft.fftshift(ifft_operator.output_array) * np.exp(1j * phase_mask[S:].reshape((Nx, Ny))))) fft_operator() + #SECOND PROPAGATION ifft_input_array[:,:] = fft_operator.output_array * P_2_nat #element wise product ifft_operator() - intermediate_fields[1] = ifft_operator.output_array.copy() + intermediate_fields[1] = np.fft.fftshift(ifft_operator.output_array.copy()) end_time = time.time() print("Propagation finished in {:.6f} seconds.".format(end_time - start_time)) @@ -413,19 +416,19 @@ def adjoint_propagate(input_field_2d,target_Efield_2d,intermediate_fields): Phi_2_dagger = np.exp(-1j * phase_mask[S:]).reshape((Nx, Ny)) #compute the adjoint source - fft_input_array[:,:] = np.sum(intermediate_fields[1] * np.conj(target_Efield_2d)) * target_Efield_2d + fft_input_array[:,:] = np.fft.ifftshift(np.sum(intermediate_fields[1] * np.conj(target_Efield_2d)) * target_Efield_2d) fft_operator() ifft_input_array[:,:] = fft_operator.output_array * P_2_dagger_nat ifft_operator() - fft_input_array[:,:] = ifft_operator.output_array * Phi_2_dagger + fft_input_array[:,:] = np.fft.ifftshift(np.fft.fftshift(ifft_operator.output_array) * Phi_2_dagger) - grad_C_phi_2 = 2 * np.real(-1j * intermediate_fields[0].T.conj() * fft_input_array).flatten() + grad_C_phi_2 = 2 * np.real(-1j * intermediate_fields[0].conj() * np.fft.fftshift(fft_input_array)).flatten() fft_operator() ifft_input_array[:,:] = fft_operator.output_array * P_1_dagger_nat #element wise ifft_operator() - grad_C_phi_1 = 2 * np.real(-1j * input_field_2d.T.conj() * ifft_operator.output_array * Phi_1_dagger).flatten() + grad_C_phi_1 = 2 * np.real(-1j * input_field_2d.conj() * np.fft.fftshift(ifft_operator.output_array) * Phi_1_dagger).flatten() end_time = time.time() print("Adjoint finished in {:.6f} seconds.".format(end_time - start_time)) @@ -433,8 +436,8 @@ def adjoint_propagate(input_field_2d,target_Efield_2d,intermediate_fields): def value_and_grad(weights): - values = [None, None] - gradients = [None, None] + values = [] + gradients = [] phase_mask[:] = phase_given_w(weights) #update phase mask in place for i, input_field in enumerate(input_field_list): @@ -445,11 +448,13 @@ def value_and_grad(weights): C_s = np.abs(np.sum(intermediate_fields[-1].flatten() * circular_mask * np.conj(target_field_list[i]))) ** 2 C = np.real(np.sum(C_s)) - values[i] = C - gradients[i] = adjoint_propagate( - input_field_2d = input_field.reshape((Nx, Ny)), - target_Efield_2d = target_field_list[i].reshape((Nx, Ny)), - intermediate_fields = intermediate_fields) + values.append(-C) + gradients.append( + adjoint_propagate( + input_field_2d = input_field.reshape((Nx, Ny)), + target_Efield_2d = target_field_list[i].reshape((Nx, Ny)), + intermediate_fields = intermediate_fields) + ) return np.array(values), np.array(gradients) @@ -470,10 +475,10 @@ def epigraph_constraint(result,epigraph_and_weights,gradient)->None: if gradient.size > 0: gradient[:,0] = -1 - gradient[:,1:] = grad * (2*np.pi) + gradient[:,1:] = - grad * (2*np.pi) - #we take the negative of obj_val because we are minimizing the epigraph - result[:] = -np.real(obj_val) - epigraph + #result of the evaluation of this constraint function + result[:] = np.real(obj_val) - epigraph opt_history = [] def obj_fun(epigraph_and_weights, grad): @@ -526,7 +531,7 @@ if __name__ == "__main__": #on metasurface 1, set initial guess on normalized phase elements w_norm[:S] += circular_mask_1 * target_pattern_1_norm_phase + circular_mask_2 * target_pattern_2_norm_phase #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 + #w_norm[S:] += np.random.rand(S) * circular_mask_1 + np.random.rand(S) * circular_mask_2 phase_mask = phase_given_w(w_norm) #make an initial plot of the phase masks @@ -563,8 +568,9 @@ if __name__ == "__main__": epigraph_tolerance, ) solver.set_maxeval(opt_max_eval) - solver.set_param("dual_ftol_rel", 1e-7) - solver.set_param("verbosity",1) + solver.set_param("dual_ftol_rel", 1e-5) + solver.set_param("verbosity", 1) + solver.set_param("inner_maxeval", 5) print("Starting optimization...") start = True @@ -594,8 +600,8 @@ if __name__ == "__main__": np.savetxt("results/modeconv1_optimized_output_field_allSource.txt", verify_out_field_all.reshape((Nx, Ny))) # Reshape fields to 2D for visualization - phase_mask_1 = phase_given_w(epigraph_and_weights[1:S]).reshape((Nx, Ny)) - phase_mask_2 = phase_given_w(epigraph_and_weights[S:]).reshape((Nx, Ny)) + 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)