edoardo/modeconverter
1
0
Fork 0
mirror of https://github.com/edoverg/modeconverter.git synced 2026-10-14 23:06:39 +00:00
Code Issues Packages Projects Releases Wiki Activity

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.

Browse Source
This commit is contained in:
edoardo
2026-07-28 13:18:22 +02:00
parent 9eb682d706
commit ebcaeb5643
1 changed files with 40 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files
metaprop_fibermode_doublet_epigraph.py
+35 -29
View File
@@ -43,7 +43,7 @@ sampling_period = xs[1] - xs[0]
d0 = 345e-6 #propagation distance: source - MS1 d0 = 345e-6 #propagation distance: source - MS1
d1 = 500e-6 #propagation distance: MS1 - MS2 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 d = [d1, d2] #d = [d1,d2] d1:distance MS1-MS2, d2: distance MS2-target
phase_min = 0 phase_min = 0
@@ -53,12 +53,13 @@ phase_max = 2 * np.pi
############################## ##############################
#optimization algorithm #optimization algorithm
epigraph_tolerance = np.array([1e-4]*2) epigraph_tolerance = np.array([1e-4]*2)
#epigraph_tolerance = np.array([1e-4])
############################## ##############################
############################## ##############################
#Source #Source
beam_waist = 9e-6 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 return gaussian_field_normalized
def propagate_source(source_to_propagate): 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.''' source for a distance d0: from its plane of definition up to MS1.'''
fft_input_array[:,:] = source_to_propagate.reshape((Nx, Ny)) fft_input_array[:,:] = source_to_propagate.reshape((Nx, Ny))
fft_operator() 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 #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, 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) 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, target_Efield_2]
#target_field_list = [target_Efield_1]
########################### ###########################
def forward_propagate(input_field:np.ndarray,weights:np.ndarray) -> np.ndarray: 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...") print("Starting propagating...")
intermediate_fields = [None, None] intermediate_fields = [None, None]
start_time = time.time() start_time = time.time()
#phase_mask = phase_given_w(weights) phase_mask = phase_given_w(weights)
fft_input_array[:,:] = (input_field * np.exp(1j * phase_mask[:S])).reshape((Nx, Ny)) fft_input_array[:,:] = np.fft.ifftshift((input_field * np.exp(1j * phase_mask[:S])).reshape((Nx, Ny)))
fft_operator() fft_operator()
#print("BEFORE")
print("BEFORE")
fft_power = np.sum(np.abs(fft_operator.output_array)**2)/S 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.where((wavelength * nu_parallel)**2 > 1, 0, 1)
non_propagating_mask = np.fft.ifftshift(non_propagating_mask) non_propagating_mask = np.fft.ifftshift(non_propagating_mask)
propagating_power = np.sum(np.abs(fft_operator.output_array * non_propagating_mask)**2)/S 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 #FIRST PROPAGATION
ifft_input_array[:,:] = fft_operator.output_array * P_1_nat #element wise product 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 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.where((wavelength * nu_parallel)**2 > 1, 0, 1)
non_propagating_mask = np.fft.ifftshift(non_propagating_mask) non_propagating_mask = np.fft.ifftshift(non_propagating_mask)
propagating_power = np.sum(np.abs(ifft_operator.input_array * non_propagating_mask)**2)/S 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() ifft_operator()
intermediate_fields[0] = ifft_operator.output_array.copy() 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)) #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() fft_operator()
#SECOND PROPAGATION #SECOND PROPAGATION
ifft_input_array[:,:] = fft_operator.output_array * P_2_nat #element wise product ifft_input_array[:,:] = fft_operator.output_array * P_2_nat #element wise product
ifft_operator() 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() end_time = time.time()
print("Propagation finished in {:.6f} seconds.".format(end_time - start_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)) Phi_2_dagger = np.exp(-1j * phase_mask[S:]).reshape((Nx, Ny))
#compute the adjoint source #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() fft_operator()
ifft_input_array[:,:] = fft_operator.output_array * P_2_dagger_nat ifft_input_array[:,:] = fft_operator.output_array * P_2_dagger_nat
ifft_operator() 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() fft_operator()
ifft_input_array[:,:] = fft_operator.output_array * P_1_dagger_nat #element wise ifft_input_array[:,:] = fft_operator.output_array * P_1_dagger_nat #element wise
ifft_operator() 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() end_time = time.time()
print("Adjoint finished in {:.6f} seconds.".format(end_time - start_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): def value_and_grad(weights):
values = [None, None] values = []
gradients = [None, None] gradients = []
phase_mask[:] = phase_given_w(weights) #update phase mask in place phase_mask[:] = phase_given_w(weights) #update phase mask in place
for i, input_field in enumerate(input_field_list): 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_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)) C = np.real(np.sum(C_s))
values[i] = C values.append(-C)
gradients[i] = adjoint_propagate( gradients.append(
adjoint_propagate(
input_field_2d = input_field.reshape((Nx, Ny)), input_field_2d = input_field.reshape((Nx, Ny)),
target_Efield_2d = target_field_list[i].reshape((Nx, Ny)), target_Efield_2d = target_field_list[i].reshape((Nx, Ny)),
intermediate_fields = intermediate_fields) intermediate_fields = intermediate_fields)
)
return np.array(values), np.array(gradients) return np.array(values), np.array(gradients)
@@ -470,10 +475,10 @@ def epigraph_constraint(result,epigraph_and_weights,gradient)->None:
if gradient.size > 0: if gradient.size > 0:
gradient[:,0] = -1 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 of the evaluation of this constraint function
result[:] = -np.real(obj_val) - epigraph result[:] = np.real(obj_val) - epigraph
opt_history = [] opt_history = []
def obj_fun(epigraph_and_weights, grad): def obj_fun(epigraph_and_weights, grad):
@@ -526,7 +531,7 @@ if __name__ == "__main__":
#on metasurface 1, set initial guess on normalized phase elements #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 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 #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) phase_mask = phase_given_w(w_norm)
#make an initial plot of the phase masks #make an initial plot of the phase masks
@@ -563,8 +568,9 @@ if __name__ == "__main__":
epigraph_tolerance, epigraph_tolerance,
) )
solver.set_maxeval(opt_max_eval) solver.set_maxeval(opt_max_eval)
solver.set_param("dual_ftol_rel", 1e-7) solver.set_param("dual_ftol_rel", 1e-5)
solver.set_param("verbosity", 1) solver.set_param("verbosity", 1)
solver.set_param("inner_maxeval", 5)
print("Starting optimization...") print("Starting optimization...")
start = True 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))) np.savetxt("results/modeconv1_optimized_output_field_allSource.txt", verify_out_field_all.reshape((Nx, Ny)))
# Reshape fields to 2D for visualization # Reshape fields to 2D for visualization
phase_mask_1 = phase_given_w(epigraph_and_weights[1: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:]).reshape((Nx, Ny)) phase_mask_2 = phase_given_w(epigraph_and_weights[S+1:]).reshape((Nx, Ny))
phase_mask_zoom = 150 phase_mask_zoom = 150
# Plot: Phase masks # 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_1, choose_quantity="phase_mask", save_name="phase_mask_1", plot_zoom_x=phase_mask_zoom, plot_zoom_y=phase_mask_zoom)