pardon - added remaining task: probabilistic non-stabilizer simulatio…

…n as mentiond in earlier conversation comment

test/Project.toml

@@ -25,5 +25,6 @@ SIMD = "fdea26ae-647d-5447-a871-4b548cad5224"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StridedViews = "4db3bf67-4bd7-4b4e-b153-31dc3fb37143"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/test_nonclifford_quantumoptics.jl

@@ -6,6 +6,7 @@ using LinearAlgebra
 using Test
 using InteractiveUtils
 using Random
+using Statistics
 
 qo_basis = SpinBasis(1//2)
 qo_tgate = sparse(identityoperator(qo_basis))
@@ -87,6 +88,18 @@ end
     end
 end
 
+function multiqubit_projrand(τ,p)
+    qo_state = Operator(τ)
+    projectrand!(τ, p)[1]
+    qo_state_after_proj = Operator(τ)
+    qo_pauli = Operator(p)
+    qo_proj1 = (identityoperator(qo_pauli) - qo_pauli)/2
+    qo_proj2 = (identityoperator(qo_pauli) + qo_pauli)/2
+    result1 = qo_proj1*qo_state*qo_proj1'
+    result2 = qo_proj2*qo_state*qo_proj2'
+    return qo_state_after_proj, result1, result2
+end
+
 @testset "Multi-qubit projections using GeneralizedStabilizer for stabilizer states" begin
     for n in 1:10
         for repetition in 1:5
@@ -95,14 +108,7 @@ end
                 p = random_pauli(n)
                 τ = state(s)
                 apply!(τ, random_clifford(n))
-                qo_state = Operator(τ)
-                projectrand!(τ, p)[1]
-                qo_state_after_proj = Operator(τ)
-                qo_pauli = Operator(p)
-                qo_proj1 = (identityoperator(qo_pauli) - qo_pauli)/2
-                qo_proj2 = (identityoperator(qo_pauli) + qo_pauli)/2
-                result1 = qo_proj1*qo_state*qo_proj1'
-                result2 = qo_proj2*qo_state*qo_proj2'
+                qo_state_after_proj, result1, result2 = multiqubit_projrand(τ,p)
                 # Normalize to ensure consistent comparison of the projected state, independent of scaling factors
                 norm_qo_state_after_proj = qo_state_after_proj/tr(qo_state_after_proj)
                 norm_result1 = result1/tr(result1)
@@ -112,3 +118,46 @@ end
         end
     end
 end
+
+function non_stabilizer_simulator(num_trials,n)
+    count = 0
+    for n in 1:n # exponential cost in this term
+        for _ in 1:num_trials
+            s = random_stabilizer(n)
+            p = random_pauli(n)
+            gs = GeneralizedStabilizer(s)
+            i = rand(1:n)
+            nc = embed(n, i, pcT) # multi-qubit random non-Clifford gate
+            apply!(gs, nc)
+            qo_state_after_proj, result1, result2 = multiqubit_projrand(gs,p)
+            # Normalize to ensure consistent comparison of the projected state, independent of scaling factors
+            norm_qo_state_after_proj = qo_state_after_proj/tr(qo_state_after_proj)
+            norm_result1 = result1/tr(result1)
+            norm_result2 = result2/tr(result2)
+            if norm_qo_state_after_proj ≈ norm_result2 || norm_qo_state_after_proj ≈ norm_result1
+                count += 1
+            end
+        end
+    end
+    prob = count/(num_trials*n)
+    return prob
+end
+
+@testset "probabilistic non-stabilizer simulator" begin
+# As described in https://www.scottaaronson.com/showcase2/report/ted-yoder.pdf,
+# "The set of states that are stabilizer circuit efficient is related to the
+# set of magic states, those states that, when combined with stabilizer circuits,
+# allow universal quantum computation. The problems are essentially complementary;
+# no magic state will be simulatable through every stabilizer circuit (assuming
+# BQP is larger than BPP, that is). With qudits of odd prime dimension, it was
+# recently found in https://arxiv.org/pdf/1201.1256 that a sufficient condition
+# for a state to be efficiently (weakly) simulated through stabilizer circuits is
+# that a certain quasi-probability representation of the state be positive."
+# Therefore, non-Clifford simulators generally require probabilistic sampling
+# techniques.
+
+    num_qubits = 10
+    num_trials = 10
+    prob = non_stabilizer_simulator(num_trials, num_qubits)
+    @test prob > 0.5
+end