plot updates for paper; debugging pyramid sim again

photonics/simulations/command_matrix.py

@@ -10,11 +10,11 @@ def make_command_matrix(
     flat_wavefront: hc.Wavefront,
     probe_amp=1e-10,
     dm_basis="modal",
-    rerun=True
+    rerun=False
 ):
     """
     Make a command matrix for a DM-WFS pair and assigns it to wfs.command_matrix
-    """  
+    """
     cmd_path = join(PROJECT_ROOT, "data", f"secondstage_{wfs.name}", f"cm_{date_now()}_{dm_basis}.npy")
     wfs.image_ref = wfs.readout(flat_wavefront)
     if (not rerun) and os.path.exists(cmd_path):

photonics/simulations/second_stage.py

@@ -21,9 +21,7 @@ def correction(
 		"dm_commands" : [],
 		"dm_shapes" : [],
 		"point_spread_functions" : [],
-		"strehl_ratios" : [],
-		"lantern_zernikes_truth" : [],
-		"lantern_zernikes_measured" : []
+		"strehl_ratios" : []
 	}
 	a = np.exp(-2 * np.pi * f_cutoff / f_loop)
 	pyramid_filter = HighPassFilter(lantern.nmodes, a)

scripts_jl/seal_linearity/interaction_matrix_231025.jl

@@ -1,10 +1,25 @@
 using NPZ
 using Plots
 using PhotonicLantern
+default(fontfamily= "Computer Modern", linewidth=3, framestyle=:box, label=nothing, grid=true)
 
 p_i_all = npzread("data/pl_231025/linearity_responses.npy")# * rad2opd
 amps = npzread("data/pl_231025/linearity_amps.npy")# * rad2opd
 
-command_matrix = make_command_matrix(amps, p_i_all, nmodes=10, push_ref=0.05)
-plot_linearity(amps, p_i_all, command_matrix, nmodes=10)
-plot_linearity(amps, p_i_all, command_matrix, nmodes=6)
+command_matrix = make_command_matrix(amps, p_i_all, nmodes=9, push_ref=0.05)
+# plot_linearity(amps, p_i_all, command_matrix, nmodes=9)
+begin
+    nmodes = 9
+    pl = []
+    for k in 1:nmodes
+        lin_k = ((p_i_all[k,:,:]) * command_matrix)[:,1:nmodes]
+        a = [(k == i ? 1 : 0.2) for i in 1:(nmodes)]'
+        pk = plot(xlabel=mode_names[k], label=nothing, legend=:outertopright, aspect_ratio=:equal, ylim=(-1, 1), xlim=(-1, 1))
+        plot!(amps, amps, ls=:dash, color=:black)
+        plot!(amps, lin_k, alpha=a)
+        push!(pl, pk)
+    end
+    p = plot(pl..., legend=nothing, size=(900,900), dpi=600, suptitle="Photonic lantern linear reconstruction, SEAL")
+    Plots.savefig("figures/linear_seal.pdf")
+    p
+end

scripts_jl/seal_linearity/interaction_svd_230714.jl

@@ -53,7 +53,7 @@ begin
     plot(pl..., layout=(3,6), plot_title="Principal lantern responses", plot_titlefontsize=13)
 end
 
-iA = p \ z
+iA = powers \ z
 
 U, S, Vt = svd(A)
 iU, iS, iVt = svd(iA)

scripts_jl/sim_linearity/interaction_matrix_sim.jl

@@ -2,11 +2,7 @@ using NPZ
 using Plots
 using PhotonicLantern
 
-pgfplotsx()
-
-plot_font = "Computer Modern"
-default(fontfamily=plot_font,
-        linewidth=2, framestyle=:box, label=nothing, grid=false)
+default(fontfamily="Computer Modern", linewidth=3, framestyle=:box, label=nothing, grid=true)
 
 amps = npzread("data/sweep_testset_amplitudes_.npy")
 p_i_all = npzread("data/sweep_testset_lanterns.npy") .|> abs2
@@ -38,12 +34,18 @@ begin
     for k in 1:nmodes
         lin_k = ((p_i_all[k,:,:]) * iA)[:,1:nmodes]
         a = [(k == i ? 1 : 0.2) for i in 1:(nmodes)]'
-        pk = plot(xlabel=mode_names[k], label=nothing, legend=:outertopright)
-        plot!(amps, lin_k, alpha=a)
+        pk = plot(xlabel=mode_names[k], label=nothing, legend=:outertopright, aspect_ratio=:equal, xlim=(-1,1), ylim=(-1,1))
         plot!(amps, amps, ls=:dash, color=:black)
+        plot!(amps, lin_k, alpha=a)
         push!(pl, pk)
     end
-    p = plot(pl..., legend=nothing, size=(750,750), dpi=600, layout=(3, 3), suptitle="Linear reconstruction, simulated data")
+    p = plot(pl..., legend=nothing, size=(900,900), dpi=600, suptitle="Photonic lantern linear reconstruction, simulated")
     Plots.savefig("figures/linear_sim.pdf")
     p
 end
+
+sweep = zeros(nmodes, size(p_i_all, 2), nmodes);
+for k in 1:nmodes
+    sweep[k,:,:] = (p_i_all[k,:,:] * iA)[:,1:nmodes]
+end
+npzwrite("data/linear_sweeps/pl_linear.npy", sweep)

scripts_jl/sim_linearity/linearity_comparison.jl

@@ -0,0 +1,25 @@
+using Plots
+using NPZ
+default(fontfamily= "Computer Modern", linewidth=3, framestyle=:box, label=nothing, grid=true)
+
+amps = npzread("data/sweep_testset_amplitudes_.npy")
+sweep_gs = npzread("data/linear_sweeps/gs.npy")
+sweep_mpwfs = npzread("data/linear_sweeps/mpwfs.npy")
+sweep_linear = npzread("data/linear_sweeps/pl_linear.npy")
+sweep_nn = npzread("data/linear_sweeps/pl_nn.npy")
+
+begin
+    nmodes = 9
+    pl = []
+    for k in 1:nmodes
+        pk = plot(xlabel=mode_names[k], legend=:bottomright, aspect_ratio=:equal, xlim=(-1,1), ylim=(-1,1))
+        plot!(amps, amps, ls=:dash, color=:black)
+        for (i, (s, n)) in enumerate(zip([sweep_mpwfs, sweep_linear, sweep_nn, sweep_gs], ["Modulated pyramid", "PL interaction matrix", "PL neural network", "PL Gerchberg-Saxton"]))
+            plot!(amps, s[k,:,k], color=15-i, label=(k == 3 ? n : nothing))
+        end
+        push!(pl, pk)
+    end
+    p = plot(pl..., size=(1000,1000), dpi=600, suptitle="All PL reconstruction algorithms compared with the modulated pyramid")
+    Plots.savefig("figures/all_recon.pdf")
+    p
+end

scripts_jl/sim_linearity/pl_gs_linearity.jl

@@ -0,0 +1,22 @@
+using Plots
+using NPZ
+default(fontfamily= "Computer Modern", linewidth=3, framestyle=:box, label=nothing, grid=true)
+
+amps = npzread("data/sweep_testset_amplitudes_.npy")
+sweep = npzread("data/linear_sweeps/gs.npy")
+
+begin
+    nmodes = 9
+    pl = []
+    for k in 1:nmodes
+        lin_k = sweep[k,:,:]
+        a = [(k == i ? 1 : 0.2) for i in 1:(nmodes)]'
+        pk = plot(xlabel=mode_names[k], label=nothing, legend=:outertopright, aspect_ratio=:equal, xlim=(-1,1), ylim=(-1,1))
+        plot!(amps, amps, ls=:dash, color=:black)
+        plot!(amps, lin_k, alpha=a)
+        push!(pl, pk)
+    end
+    p = plot(pl..., legend=nothing, size=(900,900), dpi=600, suptitle="Photonic lantern Gerchberg-Saxton reconstruction, simulated")
+    Plots.savefig("figures/gs_linearity_sim.pdf")
+    p
+end

scripts_jl/sim_linearity/pl_nn_linearity.jl

@@ -3,14 +3,10 @@ using NPZ
 using JLD2
 using Flux
 using Plots
+using Plots.PlotMeasures
 using StatsBase: mean
 
-pgfplotsx()
-
-plot_font = "Computer Modern"
-default(fontfamily=plot_font,
-        linewidth=2, framestyle=:box, label=nothing, grid=false)
-
+default(fontfamily="Computer Modern", linewidth=3, framestyle=:box, label=nothing, grid=true)
 
 zero_one_ify(x) = (x .- minimum(x)) ./ (maximum(x) .- minimum(x)), minimum(x), maximum(x)
 rescale(z, zmin, zmax) = z * (zmax - zmin) + zmin
@@ -40,14 +36,21 @@ for k in 1:nmodes
     nn_recon = rescale.(hcat([model(x) for x in eachrow(sweep_lanterns_normalized[k,:,:])]...)', xmin, xmax)
     a = [(k == i ? 1 : 0.2) for i in 1:(nmodes)]'
     pk = plot(xlabel=mode_names[k], label=nothing, legend=:outertopright, ylim=(-1, 1), xlim=(-1,1))
-    plot!(amplitudes, nn_recon, alpha=a)
     plot!(amplitudes, amplitudes, ls=:dash, color=:black)
+    plot!(amplitudes, nn_recon, alpha=a)
     push!(pl, pk)
 end
-p = plot(pl..., legend=nothing, size=(750,750), dpi=600, suptitle="Neural network reconstructor simulation, max train amplitude = $(round(xmax, digits=2)) rad", layout=(3, 3))
+p = plot(pl..., legend=nothing, size=(900,900), dpi=600, suptitle="Neural network reconstructor simulation, max train amplitude = $(round(xmax, digits=2)) rad", layout=(3, 3), left_margin=2Plots.mm)
 Plots.savefig("figures/nn_sim_0.25.pdf")
 p
 
+
+sweep = zeros(nmodes, length(amplitudes), nmodes);
+for k in 1:nmodes
+    sweep[k,:,:] = nn_recon = rescale.(hcat([model(x) for x in eachrow(sweep_lanterns_normalized[k,:,:])]...)', xmin, xmax)
+end
+npzwrite("data/linear_sweeps/pl_nn.npy", sweep)
+
 trainset_amplitudes_scaled_f32 = rescale.(trainset_amplitudes[:,48001:end], xmin, xmax)
 resid_train = rescale.(hcat(model.(eachcol(trainset_lanterns[:,48001:end]))...), xmin, xmax) - trainset_amplitudes_scaled_f32
 mean(sqrt.(sum(resid_train .^ 2, dims=1)))

scripts_py/sim/empirical_identification.ipynb

@@ -93,7 +93,7 @@
     "    ax.set_axis_off()\n",
     "plt.suptitle(\"Phase screens from diagonal basis of interaction matrix\", fontsize=28)\n",
     "plt.tight_layout()\n",
-    "plt.savefig(\"../../figures/diagonal_phase_screens.pdf\")"
+    "plt.savefig(\"../../figures/diagonal_phase_screens.p df\")"
    ]
   },
   {