test(benchmarks): PureFockSimulator benchmarks (Budapest-Quantum-Comp…

…uting-Group#217)

The benchmark `purefock_beamsplitter_increasing_cutoff_benchmark.py` is
not a realistic scenario, since only one squeezing is applied, and most
beamsplitters after that do nothing. Instead, an occupation number has
been prepared.

Moreover, a benchmark has been created testing displacement, squeezing,
interferometer and kerr together in `purefock_general_benchmark.py`.

benchmarks/purefock_beamsplitter_increasing_cutoff_benchmark.py

@@ -36,13 +36,15 @@ def d():
     return 5
 
 
-@pytest.mark.parametrize("cutoff", (3, 4, 5, 6, 7, 8))
+@pytest.mark.parametrize("cutoff", range(3, 14))
 def piquasso_benchmark(benchmark, d, cutoff, theta):
     @benchmark
     def func():
+        state_vector = [cutoff // d] * d
+        state_vector[0] += cutoff % d - 1
+
         with pq.Program() as program:
-            pq.Q() | pq.Vacuum()
-            pq.Q(1) | pq.Squeezing(0.1)
+            pq.Q(all) | pq.StateVector(state_vector)
             for i in range(d - 1):
                 pq.Q(i, i + 1) | pq.Beamsplitter(theta)
 
@@ -51,16 +53,21 @@ def func():
         simulator_fock.execute(program)
 
 
-@pytest.mark.parametrize("cutoff", (3, 4, 5, 6, 7, 8))
+@pytest.mark.parametrize("cutoff", (3, 4, 5))
 def strawberryfields_benchmark(benchmark, d, cutoff, theta):
     @benchmark
     def func():
         eng = sf.Engine(backend="fock", backend_options={"cutoff_dim": cutoff})
 
         circuit = sf.Program(d)
 
+        state_vector = [cutoff // d] * d
+        state_vector[0] += cutoff % d - 1
+
         with circuit.context as q:
-            sf.ops.Sgate(0.1) | q[1]
+            for i, n in enumerate(state_vector):
+                sf.ops.Fock(n) | q[i]
+
             for w in range(d - 1):
                 sf.ops.BSgate(theta) | (q[w], q[w + 1])
 

benchmarks/purefock_general_benchmark.py

@@ -0,0 +1,84 @@
+#
+# Copyright 2021-2023 Budapest Quantum Computing Group
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import pytest
+
+import piquasso as pq
+import strawberryfields as sf
+
+from scipy.stats import unitary_group
+
+
+pytestmark = pytest.mark.benchmark(
+    group="pure-fock-general",
+)
+
+
+@pytest.fixture
+def alpha():
+    return 0.1
+
+
+@pytest.fixture
+def r():
+    return 0.2
+
+
+@pytest.fixture
+def xi():
+    return 0.3
+
+
+parameters = [(d, unitary_group.rvs(d)) for d in range(3, 7)]
+
+
+@pytest.mark.parametrize("d, interferometer", parameters)
+def piquasso_benchmark(benchmark, d, interferometer, r, alpha, xi):
+    @benchmark
+    def func():
+        with pq.Program() as program:
+            pq.Q(all) | pq.Vacuum()
+
+            pq.Q(all) | pq.Displacement(alpha=alpha)
+            pq.Q(all) | pq.Squeezing(r)
+
+            pq.Q(all) | pq.Interferometer(interferometer)
+
+            pq.Q(all) | pq.Kerr(xi)
+
+        simulator_fock = pq.PureFockSimulator(d=d, config=pq.Config(cutoff=d))
+
+        simulator_fock.execute(program)
+
+
+@pytest.mark.parametrize("d, interferometer", parameters)
+def strawberryfields_benchmark(benchmark, d, interferometer, r, alpha, xi):
+    @benchmark
+    def func():
+        eng = sf.Engine(backend="fock", backend_options={"cutoff_dim": d})
+
+        circuit = sf.Program(d)
+
+        with circuit.context as q:
+            for i in range(d):
+                sf.ops.Dgate(alpha) | q[i]
+                sf.ops.Sgate(r) | q[i]
+
+            sf.ops.Interferometer(interferometer) | tuple(q[i] for i in range(d))
+
+            for i in range(d):
+                sf.ops.Kgate(xi) | q[i]
+
+        eng.run(circuit)

benchmarks/tf_general_benchmark.py

@@ -0,0 +1,100 @@
+#
+# Copyright 2021-2023 Budapest Quantum Computing Group
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import pytest
+
+import piquasso as pq
+import strawberryfields as sf
+import tensorflow as tf
+
+from scipy.stats import unitary_group
+
+
+pytestmark = pytest.mark.benchmark(
+    group="tf-general",
+)
+
+
+@pytest.fixture
+def alpha():
+    return 0.01
+
+
+@pytest.fixture
+def r():
+    return 0.01
+
+
+@pytest.fixture
+def xi():
+    return 0.3
+
+
+parameters = [(d, unitary_group.rvs(d)) for d in range(3, 5)]
+
+
+@pytest.mark.parametrize("d, interferometer", parameters)
+def piquasso_benchmark(benchmark, d, interferometer, alpha, r, xi):
+    @benchmark
+    def func():
+        alpha_ = tf.Variable(alpha)
+
+        with pq.Program() as program:
+            pq.Q(all) | pq.Vacuum()
+
+            pq.Q(all) | pq.Displacement(alpha=alpha_)
+            pq.Q(all) | pq.Squeezing(r)
+            pq.Q(all) | pq.Interferometer(interferometer)
+            pq.Q(all) | pq.Kerr(xi)
+
+        simulator_fock = pq.TensorflowPureFockSimulator(d=d, config=pq.Config(cutoff=d))
+
+        with tf.GradientTape() as tape:
+            state = simulator_fock.execute(program).state
+            mean_photon_number = state.mean_photon_number()
+
+        tape.gradient(mean_photon_number, [alpha_])
+
+
+@pytest.mark.parametrize("d, interferometer", parameters)
+def strawberryfields_benchmark(benchmark, d, interferometer, alpha, r, xi):
+    @benchmark
+    def func():
+        program = sf.Program(d)
+
+        mapping = {}
+
+        alpha_ = tf.Variable(alpha)
+        param = program.params("alpha")
+        mapping["alpha"] = alpha_
+
+        engine = sf.Engine(backend="tf", backend_options={"cutoff_dim": d})
+
+        with program.context as q:
+            for i in range(d):
+                sf.ops.Dgate(param) | q[i]
+                sf.ops.Sgate(r) | q[i]
+
+            sf.ops.Interferometer(interferometer) | tuple(q[i] for i in range(d))
+
+            for i in range(d):
+                sf.ops.Kgate(xi) | q[i]
+
+        with tf.GradientTape() as tape:
+            result = engine.run(program, args=mapping)
+            state = result.state
+            mean = sum([state.mean_photon(mode)[0] for mode in range(d)])
+
+        tape.gradient(mean, [alpha_])