From ca9bcee144819942de2d3ba64be57b5320fd8e8d Mon Sep 17 00:00:00 2001
From: Fe-r-oz <feroz.ahmad.email@gmail.com>
Date: Mon, 11 Nov 2024 11:32:08 +0500
Subject: [PATCH] pardon: added remaining task: probabilistic non-stabilizer
 simulation as mentiond in earlier conversation comment

---
 test/test_nonclifford_quantumoptics.jl | 69 +++++++++++++++++++++++---
 1 file changed, 61 insertions(+), 8 deletions(-)

diff --git a/test/test_nonclifford_quantumoptics.jl b/test/test_nonclifford_quantumoptics.jl
index 5dddd705f..cdbb0d850 100644
--- a/test/test_nonclifford_quantumoptics.jl
+++ b/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,50 @@ 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
+    probs = []
+    for i in 1:10
+        prob = non_stabilizer_simulator(num_trials, num_qubits)
+        push!(probs, prob)
+    end
+    @test mean(probs) > 0.5
+end