From 5cf9fcfea38b24dc91b1c8162dab5e33856137d7 Mon Sep 17 00:00:00 2001 From: benoitseron Date: Tue, 21 Mar 2023 09:49:10 +0100 Subject: [PATCH 1/4] first commit --- src/certification/gaussian_partition.jl | 129 ++++++++++++++++++++++++ 1 file changed, 129 insertions(+) create mode 100644 src/certification/gaussian_partition.jl diff --git a/src/certification/gaussian_partition.jl b/src/certification/gaussian_partition.jl new file mode 100644 index 00000000..b4f62b79 --- /dev/null +++ b/src/certification/gaussian_partition.jl @@ -0,0 +1,129 @@ +begin + using Revise + + using BosonSampling + using Plots + using ProgressMeter + using Distributions + using Random + using Test + using ArgCheck + using StatsBase + using ColorSchemes + using Interpolations + using Dierckx + using LinearAlgebra + using PrettyTables + using LaTeXStrings + using JLD + using AutoHashEquals + using LinearRegression + using DataStructures + using Parameters + using UnPack + using Optim + + + using Permanents + using Plots + using Test + using Combinatorics + using Random + using IterTools + using Statistics + using LinearAlgebra #removed so as to be able to use generic types such as BigFloats, can be put back if needed + using PolynomialRoots + using StatsBase + #using JLD + using CSV + using DataFrames + using Tables + using Plots#; plotly() #plotly allows to make "dynamic" plots where you can point the mouse and see the values, they can also be saved like that but note that this makes them heavy (and html instead of png) + using PrettyTables #to display large matrices + using Roots + using BenchmarkTools + using Optim + using ProgressMeter + using ProgressBars + using Parameters + using ArgCheck + using Distributions + using Luxor + using AutoHashEquals + using LinearRegression + using HypothesisTests + using SimpleTraits + using Parameters + using UnPack + using Dates + using JLD + using DelimitedFiles + using BosonSampling + using Plots + using ProgressMeter + using Distributions + using Random + using Test + using ArgCheck + using StatsBase + using ColorSchemes + using Interpolations + using Dierckx + using LinearAlgebra + using PrettyTables + using LaTeXStrings + using JLD + using AutoHashEquals + using LinearRegression + using DataStructures + using Parameters + using UnPack + using Combinatorics + using BosonSampling + using Plots + using ProgressMeter + using Distributions + using Random + using Test + using ArgCheck + using ColorSchemes + using Interpolations + using Dierckx + using LinearAlgebra + using PrettyTables + using LaTeXStrings + using JLD + using AutoHashEquals + using LinearRegression + using DataStructures + using Parameters + using UnPack + using Optim + using DelimitedFiles + using Dates + using DelimitedFiles + + using ITensors +end + + +# r - squeezing +# x - generatinc funtion +# η - vector of fourier coefficients +# m - modes +# Q - +# Λ - +# δ - +# λ - thermal parameters + +# the lambdas are created page 24 +# Q page 23 +# C just above + +m = 4 +r = 1.5 * ones(m) +λ = 1 * ones(m) + +C_array = [0.5 - 1/(1+ λ[j] * exp(2*r[j])) for j in 1:m] + +C = diagm(C_array) \ No newline at end of file From 4c2f61f5432fefd1b53ef40bec6287b2ac516567 Mon Sep 17 00:00:00 2001 From: benoitseron Date: Tue, 21 Mar 2023 10:18:28 +0100 Subject: [PATCH 2/4] virtual interferometer matrix --- .../gaussian_partition.jl | 66 ++++++++++++++++++- src/types/measurements.jl | 5 +- 2 files changed, 66 insertions(+), 5 deletions(-) rename src/{certification => partitions}/gaussian_partition.jl (59%) diff --git a/src/certification/gaussian_partition.jl b/src/partitions/gaussian_partition.jl similarity index 59% rename from src/certification/gaussian_partition.jl rename to src/partitions/gaussian_partition.jl index b4f62b79..01f89373 100644 --- a/src/certification/gaussian_partition.jl +++ b/src/partitions/gaussian_partition.jl @@ -115,15 +115,77 @@ end # Λ - # δ - # λ - thermal parameters +# V = U_round # the lambdas are created page 24 # Q page 23 # C just above -m = 4 +### interferometer ### +m = 2 +physical_interferometer = RandHaar(m) +U = physical_interferometer.U + +### squeezing ### r = 1.5 * ones(m) λ = 1 * ones(m) +### useful matrices ### + C_array = [0.5 - 1/(1+ λ[j] * exp(2*r[j])) for j in 1:m] +C = diagm(C_array) + +### partition ### + +part = equilibrated_partition(m, 2) + +### cutoff ### + +n_max = 10 +@warn "arbitrary cutoff" + + +### virtual interferometer matrix ### + +fourier_indexes = all_mode_configurations(n_max,part.n_subset, only_photon_number_conserving = false) +# probas_fourier = Array{ComplexF64}(undef, length(fourier_indexes)) +virtual_interferometer_matrix = similar(U) + +index_fourier_array = 1 +fourier_index = fourier_indexes[index_fourier_array] + +# for (index_fourier_array, fourier_index) in enumerate(fourier_indexes) + + # for each fourier index, we recompute the virtual interferometer + virtual_interferometer_matrix = conj(physical_interferometer.U) + + diag = [0.0 + 0im for i in 1:m] # previously 1.0 + 0im + + for (i,fourier_element) in enumerate(fourier_index) + + this_phase = exp(2*pi*1im/(n_max+1) * fourier_element) + + for j in 1:length(diag) + + if part.subsets[i].subset[j] == 1 + + diag[j] = this_phase - 1 # previously multiply by phase + + end + + end + + end + + virtual_interferometer_matrix *= Diagonal(diag) + virtual_interferometer_matrix *= conj(physical_interferometer.U') # more practical than the transpose function + + + # beware, only the modes corresponding to the + # virtual_interferometer_matrix[input_config,input_config] + # must be taken into account ! + #probas_fourier[index_fourier_array] = permanent(virtual_interferometer_matrix[mode_occupation_list,mode_occupation_list] .* S) +# end + +@warn "check if this is actually what we want or the definition of Gabrielle in U_round" -C = diagm(C_array) \ No newline at end of file diff --git a/src/types/measurements.jl b/src/types/measurements.jl index 2bb7fdb1..01154ebe 100644 --- a/src/types/measurements.jl +++ b/src/types/measurements.jl @@ -92,7 +92,7 @@ function all_threshold_detections(n, m; only_photon_number_conserving = true) return [ThresholdFockDetection(x) for x in array] - +end """ @@ -222,8 +222,7 @@ mutable struct PartitionSample <: OutputMeasurementType part_occ::Union{PartitionOccupancy, Nothing} PartitionSample() = new(nothing) PartitionSample(p::PartitionOccupancy) = new(p) -endregards - +end StateMeasurement(::Type{PartitionSample}) = PartitionMeasurement() From 805fadbda165ec51bb8f9f7f8d9529cc08ad1da2 Mon Sep 17 00:00:00 2001 From: benoitseron Date: Tue, 21 Mar 2023 12:06:04 +0100 Subject: [PATCH 3/4] most quantities --- src/partitions/gaussian_partition.jl | 44 ++++++++++++++++++++++++++-- 1 file changed, 42 insertions(+), 2 deletions(-) diff --git a/src/partitions/gaussian_partition.jl b/src/partitions/gaussian_partition.jl index 01f89373..13068f67 100644 --- a/src/partitions/gaussian_partition.jl +++ b/src/partitions/gaussian_partition.jl @@ -129,6 +129,9 @@ U = physical_interferometer.U ### squeezing ### r = 1.5 * ones(m) λ = 1 * ones(m) +displacement = (0.5 + 0.1im) * ones(m) +delta_x = real.(displacement) +delta_y = imag.(displacement) ### useful matrices ### @@ -151,7 +154,7 @@ fourier_indexes = all_mode_configurations(n_max,part.n_subset, only_photon_numbe # probas_fourier = Array{ComplexF64}(undef, length(fourier_indexes)) virtual_interferometer_matrix = similar(U) -index_fourier_array = 1 +index_fourier_array = 3 fourier_index = fourier_indexes[index_fourier_array] # for (index_fourier_array, fourier_index) in enumerate(fourier_indexes) @@ -187,5 +190,42 @@ fourier_index = fourier_indexes[index_fourier_array] #probas_fourier[index_fourier_array] = permanent(virtual_interferometer_matrix[mode_occupation_list,mode_occupation_list] .* S) # end -@warn "check if this is actually what we want or the definition of Gabrielle in U_round" +@warn "check if this is actually what we want and the definition of Gabrielle in U_round" + +### matrix Q ### + +Q = zeros(ComplexF64, 4 .* size(C)) + +Q[1:m, 1:m] = I - C +Q[1:m, 2m+1:3m] = - C - virtual_interferometer_matrix +Q[1:m, 3m+1:4m] = -1im .* virtual_interferometer_matrix + +Q[m+1:2m, m+1:2m] = I + C +Q[m+1:2m, 2m+1:3m] = -1im .* virtual_interferometer_matrix +Q[m+1:2m, 3m+1:4m] = - C + virtual_interferometer_matrix + +Q[2m+1:3m, 1:m] = - C - virtual_interferometer_matrix +Q[2m+1:3m, m+1:2m] = -1im .* virtual_interferometer_matrix +Q[2m+1:3m, 2m+1:3m] = I - C + +Q[3m+1:4m, 1:m] = -1im .* virtual_interferometer_matrix +Q[3m+1:4m, m+1:2m] = - C + virtual_interferometer_matrix +Q[3m+1:4m, 3m+1:4m] = I + C + +Q + +### Lambda matrices ### + +Λ_plus = [2 * delta_x[j] /(1+ λ[j] * exp(2*r[j])) for j in 1:m] +Λ_minus = [2 * delta_y[j] /(1+ λ[j] * exp(-2*r[j])) for j in 1:m] + +Λ = vcat(Λ_plus , Λ_minus, Λ_plus, - Λ_minus) + +### + +prod([2 / sqrt((1+ λ[j] * exp(2*r[j]))*(1+ λ[j] * exp(-2*r[j]))) for j in 1:m]) * (det(Q))^(-0.5) * 0.5 * dot(Λ, Q^(-1) * Λ) - dot(Λ_plus, delta_x) - dot(Λ_minus, delta_y) + + + + From 0f0774e0e669c8582206c79bf8e1e25f96ecfdf1 Mon Sep 17 00:00:00 2001 From: benoitseron Date: Tue, 21 Mar 2023 12:32:58 +0100 Subject: [PATCH 4/4] negative probabilities to be fixed --- src/partitions/gaussian_partition.jl | 72 +++++++++++++--------------- 1 file changed, 32 insertions(+), 40 deletions(-) diff --git a/src/partitions/gaussian_partition.jl b/src/partitions/gaussian_partition.jl index 13068f67..e559be0d 100644 --- a/src/partitions/gaussian_partition.jl +++ b/src/partitions/gaussian_partition.jl @@ -129,7 +129,7 @@ U = physical_interferometer.U ### squeezing ### r = 1.5 * ones(m) λ = 1 * ones(m) -displacement = (0.5 + 0.1im) * ones(m) +displacement = (0.1 + 0.1im) * ones(m) delta_x = real.(displacement) delta_y = imag.(displacement) @@ -138,6 +138,13 @@ delta_y = imag.(displacement) C_array = [0.5 - 1/(1+ λ[j] * exp(2*r[j])) for j in 1:m] C = diagm(C_array) +### Lambda matrices ### + +Λ_plus = [2 * delta_x[j] /(1+ λ[j] * exp(2*r[j])) for j in 1:m] +Λ_minus = [2 * delta_y[j] /(1+ λ[j] * exp(-2*r[j])) for j in 1:m] + +Λ = vcat(Λ_plus , Λ_minus, Λ_plus, - Λ_minus) + ### partition ### part = equilibrated_partition(m, 2) @@ -151,13 +158,10 @@ n_max = 10 ### virtual interferometer matrix ### fourier_indexes = all_mode_configurations(n_max,part.n_subset, only_photon_number_conserving = false) -# probas_fourier = Array{ComplexF64}(undef, length(fourier_indexes)) +probas_fourier = Array{ComplexF64}(undef, length(fourier_indexes)) virtual_interferometer_matrix = similar(U) -index_fourier_array = 3 -fourier_index = fourier_indexes[index_fourier_array] - -# for (index_fourier_array, fourier_index) in enumerate(fourier_indexes) +for (index_fourier_array, fourier_index) in enumerate(fourier_indexes) # for each fourier index, we recompute the virtual interferometer virtual_interferometer_matrix = conj(physical_interferometer.U) @@ -183,49 +187,37 @@ fourier_index = fourier_indexes[index_fourier_array] virtual_interferometer_matrix *= Diagonal(diag) virtual_interferometer_matrix *= conj(physical_interferometer.U') # more practical than the transpose function + ### matrix Q ### - # beware, only the modes corresponding to the - # virtual_interferometer_matrix[input_config,input_config] - # must be taken into account ! - #probas_fourier[index_fourier_array] = permanent(virtual_interferometer_matrix[mode_occupation_list,mode_occupation_list] .* S) -# end - -@warn "check if this is actually what we want and the definition of Gabrielle in U_round" - -### matrix Q ### - -Q = zeros(ComplexF64, 4 .* size(C)) + Q = zeros(ComplexF64, 4 .* size(C)) -Q[1:m, 1:m] = I - C -Q[1:m, 2m+1:3m] = - C - virtual_interferometer_matrix -Q[1:m, 3m+1:4m] = -1im .* virtual_interferometer_matrix + Q[1:m, 1:m] = I - C + Q[1:m, 2m+1:3m] = - C - virtual_interferometer_matrix + Q[1:m, 3m+1:4m] = -1im .* virtual_interferometer_matrix -Q[m+1:2m, m+1:2m] = I + C -Q[m+1:2m, 2m+1:3m] = -1im .* virtual_interferometer_matrix -Q[m+1:2m, 3m+1:4m] = - C + virtual_interferometer_matrix + Q[m+1:2m, m+1:2m] = I + C + Q[m+1:2m, 2m+1:3m] = -1im .* virtual_interferometer_matrix + Q[m+1:2m, 3m+1:4m] = - C + virtual_interferometer_matrix -Q[2m+1:3m, 1:m] = - C - virtual_interferometer_matrix -Q[2m+1:3m, m+1:2m] = -1im .* virtual_interferometer_matrix -Q[2m+1:3m, 2m+1:3m] = I - C + Q[2m+1:3m, 1:m] = - C - virtual_interferometer_matrix + Q[2m+1:3m, m+1:2m] = -1im .* virtual_interferometer_matrix + Q[2m+1:3m, 2m+1:3m] = I - C -Q[3m+1:4m, 1:m] = -1im .* virtual_interferometer_matrix -Q[3m+1:4m, m+1:2m] = - C + virtual_interferometer_matrix -Q[3m+1:4m, 3m+1:4m] = I + C + Q[3m+1:4m, 1:m] = -1im .* virtual_interferometer_matrix + Q[3m+1:4m, m+1:2m] = - C + virtual_interferometer_matrix + Q[3m+1:4m, 3m+1:4m] = I + C -Q - -### Lambda matrices ### - -Λ_plus = [2 * delta_x[j] /(1+ λ[j] * exp(2*r[j])) for j in 1:m] -Λ_minus = [2 * delta_y[j] /(1+ λ[j] * exp(-2*r[j])) for j in 1:m] - -Λ = vcat(Λ_plus , Λ_minus, Λ_plus, - Λ_minus) - -### + probas_fourier[index_fourier_array] = prod([2 / sqrt((1+ λ[j] * exp(2*r[j]))*(1+ λ[j] * exp(-2*r[j]))) for j in 1:m]) * (det(Q))^(-0.5) * 0.5 * dot(Λ, Q^(-1) * Λ) - dot(Λ_plus, delta_x) - dot(Λ_minus, delta_y) + +end -prod([2 / sqrt((1+ λ[j] * exp(2*r[j]))*(1+ λ[j] * exp(-2*r[j]))) for j in 1:m]) * (det(Q))^(-0.5) * 0.5 * dot(Λ, Q^(-1) * Λ) - dot(Λ_plus, delta_x) - dot(Λ_minus, delta_y) +physical_indexes = copy(fourier_indexes) +probas_physical(physical_index) = 1/(n_max+1)^(part.n_subset) * sum(probas_fourier[i] * exp(-2pi*1im/(n_max+1) * dot(physical_index, fourier_index)) for (i,fourier_index) in enumerate(fourier_indexes)) +pdf = [probas_physical(physical_index) for physical_index in physical_indexes] +# pdf = clean_pdf(pdf) +bar(real(pdf)) \ No newline at end of file