From 3ac3c2c73519391c528cb32a41db92dd9eafbace Mon Sep 17 00:00:00 2001 From: Fernando Date: Mon, 6 May 2019 08:50:03 +0200 Subject: [PATCH 1/7] Phase Estimation as a Gate in ops (#260) * First testing version of Phase Estimation with a lot of hardcoding * Adapt to All(Measure) * Adding operators for more than 1 quibit, first version * Adding operators for more than 1 quibit, first versioni: testing * Work in progress: create a PhaseX gate to tests via class * Work in progress: create a PhaseX gate to tests via class. Clean garbaje files * Work in progress: create a PhaseX gate to tests via class. Some enhanement * Work in progress: create a PhaseX gate to tests via class. PhaseX testing * Work in progress: Debugging algorithm * Work in progress: Debugging algorithm * Adding 2qubit example * adding 2 qubit Gate * Initial version * Create Phase Estimation as a new Gate in operations * Solving travis checks * python 2 compatibility + error in StatePreparation normalization * test coverage includes now the string * Improve the check test for no eigenvector test * Start modifying to decomposition * QPE as decomposition and gate * QPE as decomposition and gate: correct a detail in the test * try to get the travis-ci freeze solved * Solve a name not defined error in the phaseestimation tests * Solve coverage in tests * Address comments in review + change how to assert the tests * Enhance statistis in the tests bi more executions * Correct bad calculation in tests * Refine test * Address Andi comments: add detail in the examples and atributes and removing code in the test that is never executed --- docs/projectq.ops.rst | 1 + docs/projectq.setups.decompositions.rst | 8 + projectq/ops/__init__.py | 1 + projectq/ops/_qpegate.py | 29 ++++ projectq/ops/_qpegate_test.py | 23 +++ projectq/setups/decompositions/__init__.py | 6 +- .../setups/decompositions/phaseestimation.py | 129 ++++++++++++++ .../decompositions/phaseestimation_test.py | 162 ++++++++++++++++++ 8 files changed, 357 insertions(+), 2 deletions(-) create mode 100755 projectq/ops/_qpegate.py create mode 100755 projectq/ops/_qpegate_test.py create mode 100644 projectq/setups/decompositions/phaseestimation.py create mode 100644 projectq/setups/decompositions/phaseestimation_test.py diff --git a/docs/projectq.ops.rst b/docs/projectq.ops.rst index 8fbe2e287..58d66ee16 100755 --- a/docs/projectq.ops.rst +++ b/docs/projectq.ops.rst @@ -53,6 +53,7 @@ The operations collection consists of various default gates and is a work-in-pro projectq.ops.UniformlyControlledRy projectq.ops.UniformlyControlledRz projectq.ops.StatePreparation + projectq.ops.QPE projectq.ops.FlipBits diff --git a/docs/projectq.setups.decompositions.rst b/docs/projectq.setups.decompositions.rst index 883395db6..d900f5bfc 100755 --- a/docs/projectq.setups.decompositions.rst +++ b/docs/projectq.setups.decompositions.rst @@ -26,6 +26,7 @@ The decomposition package is a collection of gate decomposition / replacement ru projectq.setups.decompositions.time_evolution projectq.setups.decompositions.toffoli2cnotandtgate projectq.setups.decompositions.uniformlycontrolledr2cnot + projectq.setups.decompositions.phaseestimation Submodules @@ -172,6 +173,13 @@ projectq.setups.decompositions.uniformlycontrolledr2cnot module :members: :undoc-members: +projectq.setups.decompositions.phaseestimation module +--------------------------------------------------------------- + +.. automodule:: projectq.setups.decompositions.phaseestimation + :members: + :undoc-members: + Module contents --------------- diff --git a/projectq/ops/__init__.py b/projectq/ops/__init__.py index db4a38b79..388350259 100755 --- a/projectq/ops/__init__.py +++ b/projectq/ops/__init__.py @@ -37,3 +37,4 @@ from ._uniformly_controlled_rotation import (UniformlyControlledRy, UniformlyControlledRz) from ._state_prep import StatePreparation +from ._qpegate import QPE diff --git a/projectq/ops/_qpegate.py b/projectq/ops/_qpegate.py new file mode 100755 index 000000000..08beee743 --- /dev/null +++ b/projectq/ops/_qpegate.py @@ -0,0 +1,29 @@ +# Copyright 2017 ProjectQ-Framework (www.projectq.ch) +# +# 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. + +from ._basics import BasicGate + + +class QPE(BasicGate): + """ + Quantum Phase Estimation gate. + + See setups.decompositions for the complete implementation + """ + def __init__(self, unitary): + BasicGate.__init__(self) + self.unitary = unitary + + def __str__(self): + return 'QPE({})'.format(str(self.unitary)) diff --git a/projectq/ops/_qpegate_test.py b/projectq/ops/_qpegate_test.py new file mode 100755 index 000000000..5ffcbf185 --- /dev/null +++ b/projectq/ops/_qpegate_test.py @@ -0,0 +1,23 @@ +# Copyright 2017 ProjectQ-Framework (www.projectq.ch) +# +# 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. + +"""Tests for projectq.ops._qpegate.""" + +from projectq.ops import _qpegate, X + + +def test_qpe_str(): + unitary = X + gate = _qpegate.QPE(unitary) + assert str(gate) == "QPE(X)" diff --git a/projectq/setups/decompositions/__init__.py b/projectq/setups/decompositions/__init__.py index aab71b28c..db29dc7a8 100755 --- a/projectq/setups/decompositions/__init__.py +++ b/projectq/setups/decompositions/__init__.py @@ -31,7 +31,8 @@ swap2cnot, toffoli2cnotandtgate, time_evolution, - uniformlycontrolledr2cnot) + uniformlycontrolledr2cnot, + phaseestimation) all_defined_decomposition_rules = [ rule @@ -54,6 +55,7 @@ swap2cnot, toffoli2cnotandtgate, time_evolution, - uniformlycontrolledr2cnot] + uniformlycontrolledr2cnot, + phaseestimation] for rule in module.all_defined_decomposition_rules ] diff --git a/projectq/setups/decompositions/phaseestimation.py b/projectq/setups/decompositions/phaseestimation.py new file mode 100644 index 000000000..faf7523cf --- /dev/null +++ b/projectq/setups/decompositions/phaseestimation.py @@ -0,0 +1,129 @@ +# Copyright 2018 ProjectQ-Framework (www.projectq.ch) +# +# 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. + +""" +Registers a decomposition for phase estimation. + +(reference https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm) + +The Quantum Phase Estimation (QPE) executes the algorithm up to the inverse +QFT included. The following steps measuring the ancillas and computing the +phase should be executed outside of the QPE. + +The decomposition uses as ancillas (qpe_ancillas) the first qubit/qureg in +the Command and as system qubits teh second qubit/qureg in the Command. + +The unitary operator for which the phase estimation is estimated (unitary) +is the gate in Command + +Example: + .. code-block:: python + + # Example using a ProjectQ gate + + n_qpe_ancillas = 3 + qpe_ancillas = eng.allocate_qureg(n_qpe_ancillas) + system_qubits = eng.allocate_qureg(1) + angle = cmath.pi*2.*0.125 + U = Ph(angle) # unitary_specfic_to_the_problem() + + # Apply Quantum Phase Estimation + QPE(U) | (qpe_ancillas, system_qubits) + + All(Measure) | qpe_ancillas + # Compute the phase from the ancilla measurement + #(https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm) + phasebinlist = [int(q) for q in qpe_ancillas] + phase_in_bin = ''.join(str(j) for j in phasebinlist) + phase_int = int(phase_in_bin,2) + phase = phase_int / (2 ** n_qpe_ancillas) + print (phase) + + # Example using a function (two_qubit_gate). + # Instead of applying QPE on a gate U one could provide a function + + def two_qubit_gate(system_q, time): + CNOT | (system_q[0], system_q[1]) + Ph(2.0*cmath.pi*(time * 0.125)) | system_q[1] + CNOT | (system_q[0], system_q[1]) + + n_qpe_ancillas = 3 + qpe_ancillas = eng.allocate_qureg(n_qpe_ancillas) + system_qubits = eng.allocate_qureg(2) + X | system_qubits[0] + + # Apply Quantum Phase Estimation + QPE(two_qubit_gate) | (qpe_ancillas, system_qubits) + + All(Measure) | qpe_ancillas + # Compute the phase from the ancilla measurement + #(https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm) + phasebinlist = [int(q) for q in qpe_ancillas] + phase_in_bin = ''.join(str(j) for j in phasebinlist) + phase_int = int(phase_in_bin,2) + phase = phase_int / (2 ** n_qpe_ancillas) + print (phase) + +Attributes: + unitary (BasicGate): Unitary Operation either a ProjectQ gate or a function f. + Calling the function with the parameters system_qubits(Qureg) and time (integer), + i.e. f(system_qubits, time), applies to the system qubits a unitary defined in f + with parameter time. + + +""" + +import numpy as np + +from projectq.cengines import DecompositionRule +from projectq.meta import Control, Loop, get_control_count +from projectq.ops import H, Tensor, get_inverse, QFT + +from projectq.ops import QPE + + +def _decompose_QPE(cmd): + """ Decompose the Quantum Phase Estimation gate. """ + eng = cmd.engine + + # Ancillas is the first qubit/qureg. System-qubit is the second qubit/qureg + qpe_ancillas = cmd.qubits[0] + system_qubits = cmd.qubits[1] + + # Hadamard on the ancillas + Tensor(H) | qpe_ancillas + + # The Unitary Operator + U = cmd.gate.unitary + + # Control U on the system_qubits + if (callable(U)): + # If U is a function + for i in range(len(qpe_ancillas)): + with Control(eng, qpe_ancillas[i]): + U(system_qubits, time=2**i) + else: + for i in range(len(qpe_ancillas)): + ipower = int(2**i) + with Loop(eng, ipower): + with Control(eng, qpe_ancillas[i]): + U | system_qubits + + # Inverse QFT on the ancillas + get_inverse(QFT) | qpe_ancillas + +#: Decomposition rules +all_defined_decomposition_rules = [ + DecompositionRule(QPE, _decompose_QPE) +] diff --git a/projectq/setups/decompositions/phaseestimation_test.py b/projectq/setups/decompositions/phaseestimation_test.py new file mode 100644 index 000000000..828e522a4 --- /dev/null +++ b/projectq/setups/decompositions/phaseestimation_test.py @@ -0,0 +1,162 @@ +# Copyright 2017 ProjectQ-Framework (www.projectq.ch) +# +# 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. + +"Tests for projectq.setups.decompositions.phaseestimation.py." + +import cmath +import numpy as np +import pytest + +from projectq import MainEngine +from projectq.backends import Simulator +from projectq.cengines import (AutoReplacer, DecompositionRuleSet, + DummyEngine, InstructionFilter, MainEngine) + +from projectq.ops import X, H, All, Measure, Tensor, Ph, CNOT, StatePreparation + +from projectq.ops import (BasicGate) + +from projectq.ops import QPE +from projectq.setups.decompositions import phaseestimation as pe +from projectq.setups.decompositions import qft2crandhadamard as dqft +import projectq.setups.decompositions.stateprep2cnot as stateprep2cnot +import projectq.setups.decompositions.uniformlycontrolledr2cnot as ucr2cnot + + +def test_simple_test_X_eigenvectors(): + rule_set = DecompositionRuleSet(modules=[pe, dqft]) + eng = MainEngine(backend=Simulator(), + engine_list=[AutoReplacer(rule_set), + ]) + results = np.array([]) + for i in range(100): + autovector = eng.allocate_qureg(1) + X | autovector + H | autovector + unit = X + ancillas = eng.allocate_qureg(1) + QPE(unit) | (ancillas, autovector) + All(Measure) | ancillas + fasebinlist = [int(q) for q in ancillas] + fasebin = ''.join(str(j) for j in fasebinlist) + faseint = int(fasebin, 2) + phase = faseint / (2. ** (len(ancillas))) + results = np.append(results, phase) + All(Measure) | autovector + eng.flush() + + num_phase = (results == 0.5).sum() + assert num_phase/100. >= 0.35, "Statistics phase calculation are not correct (%f vs. %f)" % (num_phase/100., 0.35) + + +def test_Ph_eigenvectors(): + rule_set = DecompositionRuleSet(modules=[pe, dqft]) + eng = MainEngine(backend=Simulator(), + engine_list=[AutoReplacer(rule_set), + ]) + results = np.array([]) + for i in range(100): + autovector = eng.allocate_qureg(1) + theta = cmath.pi*2.*0.125 + unit = Ph(theta) + ancillas = eng.allocate_qureg(3) + QPE(unit) | (ancillas, autovector) + All(Measure) | ancillas + fasebinlist = [int(q) for q in ancillas] + fasebin = ''.join(str(j) for j in fasebinlist) + faseint = int(fasebin, 2) + phase = faseint / (2. ** (len(ancillas))) + results = np.append(results, phase) + All(Measure) | autovector + eng.flush() + + num_phase = (results == 0.125).sum() + assert num_phase/100. >= 0.35, "Statistics phase calculation are not correct (%f vs. %f)" % (num_phase/100., 0.35) + + +def two_qubit_gate(system_q, time): + CNOT | (system_q[0], system_q[1]) + Ph(2.0*cmath.pi*(time * 0.125)) | system_q[1] + CNOT | (system_q[0], system_q[1]) + + +def test_2qubitsPh_andfunction_eigenvectors(): + rule_set = DecompositionRuleSet(modules=[pe, dqft]) + eng = MainEngine(backend=Simulator(), + engine_list=[AutoReplacer(rule_set), + ]) + results = np.array([]) + for i in range(100): + autovector = eng.allocate_qureg(2) + X | autovector[0] + ancillas = eng.allocate_qureg(3) + QPE(two_qubit_gate) | (ancillas, autovector) + All(Measure) | ancillas + fasebinlist = [int(q) for q in ancillas] + fasebin = ''.join(str(j) for j in fasebinlist) + faseint = int(fasebin, 2) + phase = faseint / (2. ** (len(ancillas))) + results = np.append(results, phase) + All(Measure) | autovector + eng.flush() + + num_phase = (results == 0.125).sum() + assert num_phase/100. >= 0.35, "Statistics phase calculation are not correct (%f vs. %f)" % (num_phase/100., 0.35) + + +def test_X_no_eigenvectors(): + rule_set = DecompositionRuleSet(modules=[pe, dqft, stateprep2cnot, ucr2cnot]) + eng = MainEngine(backend=Simulator(), + engine_list=[AutoReplacer(rule_set), + ]) + results = np.array([]) + results_plus = np.array([]) + results_minus = np.array([]) + for i in range(100): + autovector = eng.allocate_qureg(1) + amplitude0 = (np.sqrt(2) + np.sqrt(6))/4. + amplitude1 = (np.sqrt(2) - np.sqrt(6))/4. + StatePreparation([amplitude0, amplitude1]) | autovector + unit = X + ancillas = eng.allocate_qureg(1) + QPE(unit) | (ancillas, autovector) + All(Measure) | ancillas + fasebinlist = [int(q) for q in ancillas] + fasebin = ''.join(str(j) for j in fasebinlist) + faseint = int(fasebin, 2) + phase = faseint / (2. ** (len(ancillas))) + results = np.append(results, phase) + Tensor(H) | autovector + if np.allclose(phase, .0, rtol=1e-1): + results_plus = np.append(results_plus, phase) + All(Measure) | autovector + autovector_result = int(autovector) + assert autovector_result == 0 + elif np.allclose(phase, .5, rtol=1e-1): + results_minus = np.append(results_minus, phase) + All(Measure) | autovector + autovector_result = int(autovector) + assert autovector_result == 1 + eng.flush() + + total = len(results_plus) + len(results_minus) + plus_probability = len(results_plus)/100. + assert total == pytest.approx(100, abs=5) + assert plus_probability == pytest.approx(1./4., abs = 1e-1), "Statistics on |+> probability are not correct (%f vs. %f)" % (plus_probability, 1./4.) + + +def test_string(): + unit = X + gate = QPE(unit) + assert (str(gate) == "QPE(X)") From 6a46f45bec528ccb342661a07db09cdd2d8ad8a2 Mon Sep 17 00:00:00 2001 From: Fernando Date: Thu, 18 Jul 2019 10:44:46 +0200 Subject: [PATCH 2/7] Correct statistics in qpe test (#328) * Amplitude Amplification algorithm as a Gate in ops * Amplitude Amplification algorithm as a Gate in ops, correct test_string_functions * Amplitude Amplification algorithm as a Gate in ops, correct coverage for _qaagate_test * Amplitude Amplification algorithm as a Gate in ops, correct test estimation statistics in phaseestimation_test * Try to triger Travis test because an apt get failure in the travis test * resend docs/projectq.ops.rst file * resend file versions previous to AA * Try to triger Travis test because it never ran * Try to triger Travis test again to try to get the tests ran --- projectq/setups/decompositions/phaseestimation_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/projectq/setups/decompositions/phaseestimation_test.py b/projectq/setups/decompositions/phaseestimation_test.py index 828e522a4..1b8f8e63c 100644 --- a/projectq/setups/decompositions/phaseestimation_test.py +++ b/projectq/setups/decompositions/phaseestimation_test.py @@ -112,7 +112,7 @@ def test_2qubitsPh_andfunction_eigenvectors(): eng.flush() num_phase = (results == 0.125).sum() - assert num_phase/100. >= 0.35, "Statistics phase calculation are not correct (%f vs. %f)" % (num_phase/100., 0.35) + assert num_phase/100. >= 0.34, "Statistics phase calculation are not correct (%f vs. %f)" % (num_phase/100., 0.34) def test_X_no_eigenvectors(): From f3582b7456d2ef860d810bd2ac7a272146467b18 Mon Sep 17 00:00:00 2001 From: Fernando Date: Thu, 18 Jul 2019 11:51:49 +0200 Subject: [PATCH 3/7] Amplitude Amplification algorithm as a Gate in ops (#327) * Amplitude Amplification algorithm as a Gate in ops * Amplitude Amplification algorithm as a Gate in ops, correct test_string_functions * Amplitude Amplification algorithm as a Gate in ops, correct coverage for _qaagate_test * Amplitude Amplification algorithm as a Gate in ops, correct test estimation statistics in phaseestimation_test * Try to triger Travis test because an apt get failure in the travis test * Address changes proposed by Damien * Address changes proposed by Damien, missed one * Address comments by Damien including eliminate the usage of algorith_inverse and eliminate QPE files form the PR * Address comments by Damien including eliminate the usage of algorith_inverse and eliminate QPE files form the PR, second try * Address comments by Damien forgot _qaagate_test * Update amplitudeamplification_test.py Wrap lines to 80 characters * Update amplitudeamplification.py Wrap lines to 80 characters * Update _qaagate.py Minor style correction * Update amplitudeamplification_test.py Cleanup imports --- docs/projectq.ops.rst | 1 + docs/projectq.setups.decompositions.rst | 8 + projectq/ops/__init__.py | 1 + projectq/ops/_qaagate.py | 81 ++++++++ projectq/ops/_qaagate_test.py | 27 +++ projectq/setups/decompositions/__init__.py | 6 +- .../decompositions/amplitudeamplification.py | 111 +++++++++++ .../amplitudeamplification_test.py | 182 ++++++++++++++++++ 8 files changed, 415 insertions(+), 2 deletions(-) create mode 100755 projectq/ops/_qaagate.py create mode 100755 projectq/ops/_qaagate_test.py create mode 100644 projectq/setups/decompositions/amplitudeamplification.py create mode 100644 projectq/setups/decompositions/amplitudeamplification_test.py diff --git a/docs/projectq.ops.rst b/docs/projectq.ops.rst index 58d66ee16..6b9ec5936 100755 --- a/docs/projectq.ops.rst +++ b/docs/projectq.ops.rst @@ -55,6 +55,7 @@ The operations collection consists of various default gates and is a work-in-pro projectq.ops.StatePreparation projectq.ops.QPE projectq.ops.FlipBits + projectq.ops.QAA Module contents diff --git a/docs/projectq.setups.decompositions.rst b/docs/projectq.setups.decompositions.rst index d900f5bfc..26cd392fc 100755 --- a/docs/projectq.setups.decompositions.rst +++ b/docs/projectq.setups.decompositions.rst @@ -27,6 +27,7 @@ The decomposition package is a collection of gate decomposition / replacement ru projectq.setups.decompositions.toffoli2cnotandtgate projectq.setups.decompositions.uniformlycontrolledr2cnot projectq.setups.decompositions.phaseestimation + projectq.setups.decompositions.amplitudeamplification Submodules @@ -180,6 +181,13 @@ projectq.setups.decompositions.phaseestimation module :members: :undoc-members: +projectq.setups.decompositions.amplitudeamplification module +--------------------------------------------------------------- + +.. automodule:: projectq.setups.decompositions.amplitudeamplification + :members: + :undoc-members: + Module contents --------------- diff --git a/projectq/ops/__init__.py b/projectq/ops/__init__.py index 388350259..cac384d9e 100755 --- a/projectq/ops/__init__.py +++ b/projectq/ops/__init__.py @@ -38,3 +38,4 @@ UniformlyControlledRz) from ._state_prep import StatePreparation from ._qpegate import QPE +from ._qaagate import QAA diff --git a/projectq/ops/_qaagate.py b/projectq/ops/_qaagate.py new file mode 100755 index 000000000..751b9dc60 --- /dev/null +++ b/projectq/ops/_qaagate.py @@ -0,0 +1,81 @@ +# Copyright 2019 ProjectQ-Framework (www.projectq.ch) +# +# 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. + +from ._basics import BasicGate + + +class QAA(BasicGate): + """ + Quantum Aplitude Aplification gate. + + (Quick reference https://en.wikipedia.org/wiki/Amplitude_amplification. + Complete reference G. Brassard, P. Hoyer, M. Mosca, A. Tapp (2000) + Quantum Amplitude Amplification and Estimation + https://arxiv.org/abs/quant-ph/0005055) + + Quantum Amplitude Amplification (QAA) executes the algorithm, but not + the final measurement required to obtain the marked state(s) with high + probability. The starting state on wich the QAA algorithm is executed + is the one resulting of aplying the Algorithm on the |0> state. + + Example: + .. code-block:: python + + def func_algorithm(eng,system_qubits): + All(H) | system_qubits + + def func_oracle(eng,system_qubits,qaa_ancilla): + # This oracle selects the state |010> as the one marked + with Compute(eng): + All(X) | system_qubits[0::2] + with Control(eng, system_qubits): + X | qaa_ancilla + Uncompute(eng) + + system_qubits = eng.allocate_qureg(3) + # Prepare the qaa_ancilla qubit in the |-> state + qaa_ancilla = eng.allocate_qubit() + X | qaa_ancilla + H | qaa_ancilla + + # Creates the initial state form the Algorithm + func_algorithm(eng, system_qubits) + # Apply Quantum Amplitude Amplification the correct number of times + num_it = int(math.pi/4.*math.sqrt(1 << 3)) + with Loop(eng, num_it): + QAA(func_algorithm, func_oracle) | (system_qubits, qaa_ancilla) + + All(Measure) | system_qubits + + Warning: + No qubit allocation/deallocation may take place during the call + to the defined Algorithm :code:`func_algorithm` + + Attributes: + func_algorithm: Algorithm that initialite the state and to be used + in the QAA algorithm + func_oracle: The Oracle that marks the state(s) as "good" + system_qubits: the system we are interested on + qaa_ancilla: auxiliary qubit that helps to invert the amplitude of the + "good" states + + """ + def __init__(self, algorithm, oracle): + BasicGate.__init__(self) + self.algorithm = algorithm + self.oracle = oracle + + def __str__(self): + return 'QAA(Algorithm = {0}, Oracle = {1})'.format( + str(self.algorithm.__name__), str(self.oracle.__name__)) diff --git a/projectq/ops/_qaagate_test.py b/projectq/ops/_qaagate_test.py new file mode 100755 index 000000000..3e15e6801 --- /dev/null +++ b/projectq/ops/_qaagate_test.py @@ -0,0 +1,27 @@ +# Copyright 2019 ProjectQ-Framework (www.projectq.ch) +# +# 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. + +"""Tests for projectq.ops._qaagate.""" + +from projectq.ops import _qaagate, All, H, X + + +def test_qaa_str(): + + def func_algorithm(): All(H) + + def func_oracle(): All(X) + + gate = _qaagate.QAA(func_algorithm, func_oracle) + assert str(gate) == "QAA(Algorithm = func_algorithm, Oracle = func_oracle)" diff --git a/projectq/setups/decompositions/__init__.py b/projectq/setups/decompositions/__init__.py index db29dc7a8..883f04581 100755 --- a/projectq/setups/decompositions/__init__.py +++ b/projectq/setups/decompositions/__init__.py @@ -32,7 +32,8 @@ toffoli2cnotandtgate, time_evolution, uniformlycontrolledr2cnot, - phaseestimation) + phaseestimation, + amplitudeamplification) all_defined_decomposition_rules = [ rule @@ -56,6 +57,7 @@ toffoli2cnotandtgate, time_evolution, uniformlycontrolledr2cnot, - phaseestimation] + phaseestimation, + amplitudeamplification] for rule in module.all_defined_decomposition_rules ] diff --git a/projectq/setups/decompositions/amplitudeamplification.py b/projectq/setups/decompositions/amplitudeamplification.py new file mode 100644 index 000000000..517aadeed --- /dev/null +++ b/projectq/setups/decompositions/amplitudeamplification.py @@ -0,0 +1,111 @@ +# Copyright 2019 ProjectQ-Framework (www.projectq.ch) +# +# 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. +""" +Registers a decomposition for quantum amplitude amplification. + +(Quick reference https://en.wikipedia.org/wiki/Amplitude_amplification. +Complete reference G. Brassard, P. Hoyer, M. Mosca, A. Tapp (2000) +Quantum Amplitude Amplification and Estimation +https://arxiv.org/abs/quant-ph/0005055) + +Quantum Amplitude Amplification (QAA) executes the algorithm, but not +the final measurement required to obtain the marked state(s) with high +probability. The starting state on wich the QAA algorithm is executed +is the one resulting of aplying the Algorithm on the |0> state. + +Example: + .. code-block:: python + + def func_algorithm(eng,system_qubits): + All(H) | system_qubits + + def func_oracle(eng,system_qubits,qaa_ancilla): + # This oracle selects the state |010> as the one marked + with Compute(eng): + All(X) | system_qubits[0::2] + with Control(eng, system_qubits): + X | qaa_ancilla + Uncompute(eng) + + system_qubits = eng.allocate_qureg(3) + # Prepare the qaa_ancilla qubit in the |-> state + qaa_ancilla = eng.allocate_qubit() + X | qaa_ancilla + H | qaa_ancilla + + # Creates the initial state form the Algorithm + func_algorithm(eng, system_qubits) + # Apply Quantum Amplitude Amplification the correct number of times + num_it = int(math.pi/4.*math.sqrt(1 << 3)) + with Loop(eng, num_it): + QAA(func_algorithm, func_oracle) | (system_qubits, qaa_ancilla) + + All(Measure) | system_qubits + +Warning: + No qubit allocation/deallocation may take place during the call + to the defined Algorithm :code:`func_algorithm` + +Attributes: + func_algorithm: Algorithm that initialite the state and to be used + in the QAA algorithm + func_oracle: The Oracle that marks the state(s) as "good" + system_qubits: the system we are interested on + qaa_ancilla: auxiliary qubit that helps to invert the amplitude of the + "good" states + +""" + +import math +import numpy as np + +from projectq.cengines import DecompositionRule +from projectq.meta import Control, Compute, Uncompute, CustomUncompute, Dagger +from projectq.ops import X, Z, Ph, All + +from projectq.ops import QAA + + +def _decompose_QAA(cmd): + """ Decompose the Quantum Amplitude Apmplification algorithm as a gate. """ + eng = cmd.engine + + # System-qubit is the first qubit/qureg. Ancilla qubit is the second qubit + system_qubits = cmd.qubits[0] + qaa_ancilla = cmd.qubits[1] + + # The Oracle and the Algorithm + Oracle = cmd.gate.oracle + A = cmd.gate.algorithm + + # Apply the oracle to invert the amplitude of the good states, S_Chi + Oracle(eng, system_qubits, qaa_ancilla) + + # Apply the inversion of the Algorithm, + # the inversion of the aplitude of |0> and the Algorithm + + with Compute(eng): + with Dagger(eng): + A(eng, system_qubits) + All(X) | system_qubits + with Control(eng, system_qubits[0:-1]): + Z | system_qubits[-1] + with CustomUncompute(eng): + All(X) | system_qubits + A(eng, system_qubits) + Ph(math.pi) | system_qubits[0] + + +#: Decomposition rules +all_defined_decomposition_rules = [DecompositionRule(QAA, _decompose_QAA)] diff --git a/projectq/setups/decompositions/amplitudeamplification_test.py b/projectq/setups/decompositions/amplitudeamplification_test.py new file mode 100644 index 000000000..f99681713 --- /dev/null +++ b/projectq/setups/decompositions/amplitudeamplification_test.py @@ -0,0 +1,182 @@ +# Copyright 2019 ProjectQ-Framework (www.projectq.ch) +# +# 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. + +"Tests for projectq.setups.decompositions.amplitudeamplification.py." + +import math +import pytest + +from projectq import MainEngine +from projectq.backends import Simulator +from projectq.cengines import (AutoReplacer, DecompositionRuleSet, MainEngine) + +from projectq.ops import (X, H, Ry, All, Measure) +from projectq.meta import Loop, Control, Compute, Uncompute + +from projectq.ops import QAA +from projectq.setups.decompositions import amplitudeamplification as aa + + +def hache_algorithm(eng, qreg): + All(H) | qreg + + +def simple_oracle(eng, system_q, control): + # This oracle selects the state |1010101> as the one marked + with Compute(eng): + All(X) | system_q[1::2] + with Control(eng, system_q): + X | control + Uncompute(eng) + + +def test_simple_grover(): + rule_set = DecompositionRuleSet(modules=[aa]) + + eng = MainEngine(backend=Simulator(), + engine_list=[ + AutoReplacer(rule_set), + ]) + + system_qubits = eng.allocate_qureg(7) + + # Prepare the control qubit in the |-> state + control = eng.allocate_qubit() + X | control + H | control + + # Creates the initial state form the Algorithm + hache_algorithm(eng, system_qubits) + + # Get the amplitude of the marked state before the AA + # to calculate the number of iterations + eng.flush() + prob1010101 = eng.backend.get_probability('1010101', system_qubits) + + total_amp_before = math.sqrt(prob1010101) + theta_before = math.asin(total_amp_before) + + # Apply Quantum Amplitude Amplification the correct number of times + # Theta is calculated previously using get_probability + # We calculate also the theoretical final probability + # of getting the good state + num_it = int(math.pi / (4. * theta_before) + 1) + theoretical_prob = math.sin((2 * num_it + 1.) * theta_before)**2 + with Loop(eng, num_it): + QAA(hache_algorithm, simple_oracle) | (system_qubits, control) + + # Get the probabilty of getting the marked state after the AA + # to compare with the theoretical probability after teh AA + eng.flush() + prob1010101 = eng.backend.get_probability('1010101', system_qubits) + total_prob_after = prob1010101 + + All(Measure) | system_qubits + H | control + Measure | control + result = [int(q) for q in system_qubits] + control_result = int(control) + + eng.flush() + + assert total_prob_after == pytest.approx(theoretical_prob, abs=1e-6), ( + "The obtained probability is less than expected %f vs. %f" % + (total_prob_after, theoretical_prob)) + + +def complex_algorithm(eng, qreg): + All(H) | qreg + with Control(eng, qreg[0]): + All(X) | qreg[1:] + All(Ry(math.pi / 4)) | qreg[1:] + with Control(eng, qreg[-1]): + All(X) | qreg[1:-1] + + +def complex_oracle(eng, system_q, control): + # This oracle selects the subspace |000000>+|111111> as the good one + with Compute(eng): + with Control(eng, system_q[0]): + All(X) | system_q[1:] + H | system_q[0] + All(X) | system_q + + with Control(eng, system_q): + X | control + + Uncompute(eng) + + +def test_complex_aa(): + rule_set = DecompositionRuleSet(modules=[aa]) + + eng = MainEngine(backend=Simulator(), + engine_list=[ + AutoReplacer(rule_set), + ]) + + system_qubits = eng.allocate_qureg(6) + + # Prepare the control qubit in the |-> state + control = eng.allocate_qubit() + X | control + H | control + + # Creates the initial state form the Algorithm + complex_algorithm(eng, system_qubits) + + # Get the probabilty of getting the marked state before the AA + # to calculate the number of iterations + eng.flush() + prob000000 = eng.backend.get_probability('000000', system_qubits) + prob111111 = eng.backend.get_probability('111111', system_qubits) + + total_amp_before = math.sqrt(prob000000 + prob111111) + theta_before = math.asin(total_amp_before) + + # Apply Quantum Amplitude Amplification the correct number of times + # Theta is calculated previously using get_probability + # We calculate also the theoretical final probability + # of getting the good state + num_it = int(math.pi / (4. * theta_before) + 1) + theoretical_prob = math.sin((2 * num_it + 1.) * theta_before)**2 + with Loop(eng, num_it): + QAA(complex_algorithm, complex_oracle) | (system_qubits, control) + + # Get the probabilty of getting the marked state after the AA + # to compare with the theoretical probability after the AA + eng.flush() + prob000000 = eng.backend.get_probability('000000', system_qubits) + prob111111 = eng.backend.get_probability('111111', system_qubits) + total_prob_after = prob000000 + prob111111 + + All(Measure) | system_qubits + H | control + Measure | control + result = [int(q) for q in system_qubits] + control_result = int(control) + + eng.flush() + + assert total_prob_after == pytest.approx(theoretical_prob, abs=1e-2), ( + "The obtained probability is less than expected %f vs. %f" % + (total_prob_after, theoretical_prob)) + + +def test_string_functions(): + algorithm = hache_algorithm + oracle = simple_oracle + gate = QAA(algorithm, oracle) + assert (str(gate) == + "QAA(Algorithm = hache_algorithm, Oracle = simple_oracle)") From 9db33835239eeaf1e1eb1a318ab3a471646d726d Mon Sep 17 00:00:00 2001 From: David Wierichs Date: Wed, 14 Aug 2019 10:11:13 +0200 Subject: [PATCH 4/7] 3 additional 2-qubit gates (#330) * Modified _gates.py: Documentation, 2-qubit rotations, 1qubit-rotation string attributes. * Strings of rotation gates fixed. * Added two-qubit rotation gate tests. * Resource Counter import Rzz added. * Added Rzz test and fixed expected outcome. * removed wrongfully pushed dev gates. * Update _gates.py Remove unneeded import * Removed hardcoded "Phase" name for Ph(angle) gate --- projectq/backends/_resource_test.py | 5 +-- projectq/ops/_gates.py | 48 +++++++++++++++++++++++++++-- projectq/ops/_gates_test.py | 36 ++++++++++++++++++++++ 3 files changed, 84 insertions(+), 5 deletions(-) diff --git a/projectq/backends/_resource_test.py b/projectq/backends/_resource_test.py index 664b687ad..cf7122c01 100755 --- a/projectq/backends/_resource_test.py +++ b/projectq/backends/_resource_test.py @@ -20,7 +20,7 @@ from projectq.cengines import DummyEngine, MainEngine, NotYetMeasuredError from projectq.meta import LogicalQubitIDTag -from projectq.ops import All, Allocate, CNOT, Command, H, Measure, QFT, Rz, X +from projectq.ops import All, Allocate, CNOT, Command, H, Measure, QFT, Rz, Rzz, X from projectq.types import WeakQubitRef from projectq.backends import ResourceCounter @@ -74,6 +74,7 @@ def test_resource_counter(): CNOT | (qubit1, qubit3) Rz(0.1) | qubit1 Rz(0.3) | qubit1 + Rzz(0.5) | qubit1 All(Measure) | qubit1 + qubit3 @@ -81,7 +82,7 @@ def test_resource_counter(): int(qubit1) assert resource_counter.max_width == 2 - assert resource_counter.depth_of_dag == 5 + assert resource_counter.depth_of_dag == 6 str_repr = str(resource_counter) assert str_repr.count(" HGate : 1") == 1 diff --git a/projectq/ops/_gates.py b/projectq/ops/_gates.py index d524e4b3b..8967b1da5 100755 --- a/projectq/ops/_gates.py +++ b/projectq/ops/_gates.py @@ -16,17 +16,29 @@ Contains definitions of standard gates such as * Hadamard (H) * Pauli-X (X / NOT) +* Pauli-Y (Y) * Pauli-Z (Z) +* S and its inverse (S / Sdagger) * T and its inverse (T / Tdagger) +* SqrtX gate (SqrtX) * Swap gate (Swap) +* SqrtSwap gate (SqrtSwap) +* Entangle (Entangle) * Phase gate (Ph) +* Rotation-X (Rx) +* Rotation-Y (Ry) * Rotation-Z (Rz) +* Rotation-XX on two qubits (Rxx) +* Rotation-YY on two qubits (Ryy) +* Rotation-ZZ on two qubits (Rzz) * Phase-shift (R) * Measurement (Measure) and meta gates, i.e., * Allocate / Deallocate qubits * Flush gate (end of circuit) +* Barrier +* FlipBits """ import math @@ -110,7 +122,7 @@ def __str__(self): #: Shortcut (instance of) :class:`projectq.ops.SGate` S = SGate() -#: Shortcut (instance of) :class:`projectq.ops.SGate` +#: Inverse (and shortcut) of :class:`projectq.ops.SGate` Sdag = Sdagger = get_inverse(S) @@ -125,7 +137,7 @@ def __str__(self): #: Shortcut (instance of) :class:`projectq.ops.TGate` T = TGate() -#: Shortcut (instance of) :class:`projectq.ops.TGate` +#: Inverse (and shortcut) of :class:`projectq.ops.TGate` Tdag = Tdagger = get_inverse(T) @@ -217,7 +229,7 @@ def matrix(self): class Ry(BasicRotationGate): - """ RotationX gate class """ + """ RotationY gate class """ @property def matrix(self): return np.matrix([[math.cos(0.5 * self.angle), @@ -234,6 +246,36 @@ def matrix(self): [0, cmath.exp(.5 * 1j * self.angle)]]) +class Rxx(BasicRotationGate): + """ RotationXX gate class """ + @property + def matrix(self): + return np.matrix([[cmath.cos(.5 * self.angle), 0, 0, -1j*cmath.sin(.5 * self.angle)], + [0, cmath.cos( .5 * self.angle), -1j*cmath.sin(.5 * self.angle), 0], + [0, -1j*cmath.sin(.5 * self.angle), cmath.cos( .5 * self.angle), 0], + [-1j*cmath.sin(.5 * self.angle), 0, 0, cmath.cos( .5 * self.angle)]]) + + +class Ryy(BasicRotationGate): + """ RotationYY gate class """ + @property + def matrix(self): + return np.matrix([[cmath.cos(.5 * self.angle), 0, 0, 1j*cmath.sin(.5 * self.angle)], + [0, cmath.cos( .5 * self.angle), -1j*cmath.sin(.5 * self.angle), 0], + [0, -1j*cmath.sin(.5 * self.angle), cmath.cos( .5 * self.angle), 0], + [1j*cmath.sin(.5 * self.angle), 0, 0, cmath.cos( .5 * self.angle)]]) + + +class Rzz(BasicRotationGate): + """ RotationZZ gate class """ + @property + def matrix(self): + return np.matrix([[cmath.exp(-.5 * 1j * self.angle), 0, 0, 0], + [0, cmath.exp( .5 * 1j * self.angle), 0, 0], + [0, 0, cmath.exp( .5 * 1j * self.angle), 0], + [0, 0, 0, cmath.exp(-.5 * 1j * self.angle)]]) + + class R(BasicPhaseGate): """ Phase-shift gate (equivalent to Rz up to a global phase) """ @property diff --git a/projectq/ops/_gates_test.py b/projectq/ops/_gates_test.py index efcd63b0a..88efa3a19 100755 --- a/projectq/ops/_gates_test.py +++ b/projectq/ops/_gates_test.py @@ -156,6 +156,42 @@ def test_rz(angle): assert np.allclose(gate.matrix, expected_matrix) +@pytest.mark.parametrize("angle", [0, 0.2, 2.1, 4.1, 2 * math.pi, + 4 * math.pi]) +def test_rxx(angle): + gate = _gates.Rxx(angle) + expected_matrix = np.matrix([[cmath.cos(.5 * angle), 0, 0, -1j * cmath.sin(.5 * angle)], + [0, cmath.cos(.5 * angle), -1j * cmath.sin(.5 * angle), 0], + [0, -1j * cmath.sin(.5 * angle), cmath.cos(.5 * angle), 0], + [-1j * cmath.sin(.5 * angle), 0, 0, cmath.cos(.5 * angle)]]) + assert gate.matrix.shape == expected_matrix.shape + assert np.allclose(gate.matrix, expected_matrix) + + +@pytest.mark.parametrize("angle", [0, 0.2, 2.1, 4.1, 2 * math.pi, + 4 * math.pi]) +def test_ryy(angle): + gate = _gates.Ryy(angle) + expected_matrix = np.matrix([[cmath.cos(.5 * angle), 0, 0, 1j * cmath.sin(.5 * angle)], + [0, cmath.cos(.5 * angle), -1j * cmath.sin(.5 * angle), 0], + [0, -1j * cmath.sin(.5 * angle), cmath.cos(.5 * angle), 0], + [ 1j * cmath.sin(.5 * angle), 0, 0, cmath.cos(.5 * angle)]]) + assert gate.matrix.shape == expected_matrix.shape + assert np.allclose(gate.matrix, expected_matrix) + + +@pytest.mark.parametrize("angle", [0, 0.2, 2.1, 4.1, 2 * math.pi, + 4 * math.pi]) +def test_rzz(angle): + gate = _gates.Rzz(angle) + expected_matrix = np.matrix([[cmath.exp(-.5 * 1j * angle), 0, 0, 0], + [0, cmath.exp( .5 * 1j * angle), 0, 0], + [0, 0, cmath.exp( .5 * 1j * angle), 0], + [0, 0, 0, cmath.exp(-.5 * 1j * angle)]]) + assert gate.matrix.shape == expected_matrix.shape + assert np.allclose(gate.matrix, expected_matrix) + + @pytest.mark.parametrize("angle", [0, 0.2, 2.1, 4.1, 2 * math.pi]) def test_ph(angle): gate = _gates.Ph(angle) From f27286f520ff8b1ca2dec9142b01b59a2490fcd5 Mon Sep 17 00:00:00 2001 From: Melven Roehrig-Zoellner Date: Tue, 20 Aug 2019 10:12:09 +0200 Subject: [PATCH 5/7] C++ simulator performance improvements (#329) * C++ simulator performance: make the swap-gate run in native C++ It was defined as a BasicMathGate before which made it run as python code through the emulate_math_wrapper. The new variant just uses its matrix representation to run it in native code. * C++ simulator performance: add dedicated C++ code for common math gates The BasicMathGate uses a C++ python wrapper (emulate_math_wrapper) to allow generic calculations which makes it very slow. This detects some math gates and provides a native C++ implementation for it. * C++ simulator performance: use larger memory alignment * C++ simulator performance: recycle large StateVector memory buffers This avoids costly std::vector copying/reallocations by using some static std::vector to reuse the allocated buffer (just by std::swap'ing a vector into a buffer for later use when it would be deallocated otherwise). * C++ simulator performance: improve compiler flags * Add test coverage for constant math emulation * Revert "Add test coverage for constant math emulation" This reverts commit 3bb8a2cc7fd595db48b0f4d260124ccfe60d7fcf. * Add test coverage for constant math emulation --- .../backends/_sim/_cppkernels/simulator.hpp | 154 +++++++++++++----- projectq/backends/_sim/_cppsim.cpp | 3 + projectq/backends/_sim/_simulator.py | 28 +++- projectq/backends/_sim/_simulator_test.py | 52 ++++++ projectq/ops/_gates.py | 3 +- setup.py | 1 + 6 files changed, 196 insertions(+), 45 deletions(-) diff --git a/projectq/backends/_sim/_cppkernels/simulator.hpp b/projectq/backends/_sim/_cppkernels/simulator.hpp index c4e611eba..d248ed038 100755 --- a/projectq/backends/_sim/_cppkernels/simulator.hpp +++ b/projectq/backends/_sim/_cppkernels/simulator.hpp @@ -38,7 +38,7 @@ class Simulator{ public: using calc_type = double; using complex_type = std::complex; - using StateVector = std::vector>; + using StateVector = std::vector>; using Map = std::map; using RndEngine = std::mt19937; using Term = std::vector>; @@ -55,11 +55,18 @@ class Simulator{ void allocate_qubit(unsigned id){ if (map_.count(id) == 0){ map_[id] = N_++; - auto newvec = StateVector(1UL << N_); - #pragma omp parallel for schedule(static) + StateVector newvec; // avoid large memory allocations + if( tmpBuff1_.capacity() >= (1UL << N_) ) + std::swap(newvec, tmpBuff1_); + newvec.resize(1UL << N_); +#pragma omp parallel for schedule(static) for (std::size_t i = 0; i < newvec.size(); ++i) newvec[i] = (i < vec_.size())?vec_[i]:0.; - vec_ = std::move(newvec); + std::swap(vec_, newvec); + // recycle large memory + std::swap(tmpBuff1_, newvec); + if( tmpBuff1_.capacity() < tmpBuff2_.capacity() ) + std::swap(tmpBuff1_, tmpBuff2_); } else throw(std::runtime_error( @@ -113,12 +120,18 @@ class Simulator{ } } else{ - StateVector newvec((1UL << (N_-1))); - #pragma omp parallel for schedule(static) + StateVector newvec; // avoid costly memory reallocations + if( tmpBuff1_.capacity() >= (1UL << (N_-1)) ) + std::swap(tmpBuff1_, newvec); + newvec.resize((1UL << (N_-1))); + #pragma omp parallel for schedule(static) if(0) for (std::size_t i = 0; i < vec_.size(); i += 2*delta) std::copy_n(&vec_[i + static_cast(value)*delta], delta, &newvec[i/2]); - vec_ = std::move(newvec); + std::swap(vec_, newvec); + std::swap(tmpBuff1_, newvec); + if( tmpBuff1_.capacity() < tmpBuff2_.capacity() ) + std::swap(tmpBuff1_, tmpBuff2_); for (auto& p : map_){ if (p.second > pos) @@ -189,8 +202,8 @@ class Simulator{ } template - void apply_controlled_gate(M const& m, std::vector ids, - std::vector ctrl){ + void apply_controlled_gate(M const& m, const std::vector& ids, + const std::vector& ctrl){ auto fused_gates = fused_gates_; fused_gates.insert(m, ids, ctrl); @@ -209,8 +222,8 @@ class Simulator{ } template - void emulate_math(F const& f, QuReg quregs, std::vector ctrl, - unsigned num_threads=1){ + void emulate_math(F const& f, QuReg quregs, const std::vector& ctrl, + bool parallelize = false){ run(); auto ctrlmask = get_control_mask(ctrl); @@ -218,37 +231,76 @@ class Simulator{ for (unsigned j = 0; j < quregs[i].size(); ++j) quregs[i][j] = map_[quregs[i][j]]; - StateVector newvec(vec_.size(), 0.); - std::vector res(quregs.size()); - - #pragma omp parallel for schedule(static) firstprivate(res) num_threads(num_threads) - for (std::size_t i = 0; i < vec_.size(); ++i){ - if ((ctrlmask&i) == ctrlmask){ - for (unsigned qr_i = 0; qr_i < quregs.size(); ++qr_i){ - res[qr_i] = 0; - for (unsigned qb_i = 0; qb_i < quregs[qr_i].size(); ++qb_i) - res[qr_i] |= ((i >> quregs[qr_i][qb_i])&1) << qb_i; - } - f(res); - auto new_i = i; - for (unsigned qr_i = 0; qr_i < quregs.size(); ++qr_i){ - for (unsigned qb_i = 0; qb_i < quregs[qr_i].size(); ++qb_i){ - if (!(((new_i >> quregs[qr_i][qb_i])&1) == ((res[qr_i] >> qb_i)&1))) - new_i ^= (1UL << quregs[qr_i][qb_i]); - } - } - newvec[new_i] += vec_[i]; - } - else - newvec[i] += vec_[i]; + StateVector newvec; // avoid costly memory reallocations + if( tmpBuff1_.capacity() >= vec_.size() ) + std::swap(newvec, tmpBuff1_); + newvec.resize(vec_.size()); +#pragma omp parallel for schedule(static) + for (std::size_t i = 0; i < vec_.size(); i++) + newvec[i] = 0; + +//#pragma omp parallel reduction(+:newvec[:newvec.size()]) if(parallelize) // requires OpenMP 4.5 + { + std::vector res(quregs.size()); + //#pragma omp for schedule(static) + for (std::size_t i = 0; i < vec_.size(); ++i){ + if ((ctrlmask&i) == ctrlmask){ + for (unsigned qr_i = 0; qr_i < quregs.size(); ++qr_i){ + res[qr_i] = 0; + for (unsigned qb_i = 0; qb_i < quregs[qr_i].size(); ++qb_i) + res[qr_i] |= ((i >> quregs[qr_i][qb_i])&1) << qb_i; + } + f(res); + auto new_i = i; + for (unsigned qr_i = 0; qr_i < quregs.size(); ++qr_i){ + for (unsigned qb_i = 0; qb_i < quregs[qr_i].size(); ++qb_i){ + if (!(((new_i >> quregs[qr_i][qb_i])&1) == ((res[qr_i] >> qb_i)&1))) + new_i ^= (1UL << quregs[qr_i][qb_i]); + } + } + newvec[new_i] += vec_[i]; + } + else + newvec[i] += vec_[i]; + } } - vec_ = std::move(newvec); + std::swap(vec_, newvec); + std::swap(tmpBuff1_, newvec); + } + + // faster version without calling python + template + inline void emulate_math_addConstant(int a, const QuReg& quregs, const std::vector& ctrl) + { + emulate_math([a](std::vector &res){for(auto& x: res) x = x + a;}, quregs, ctrl, true); + } + + // faster version without calling python + template + inline void emulate_math_addConstantModN(int a, int N, const QuReg& quregs, const std::vector& ctrl) + { + emulate_math([a,N](std::vector &res){for(auto& x: res) x = (x + a) % N;}, quregs, ctrl, true); + } + + // faster version without calling python + template + inline void emulate_math_multiplyByConstantModN(int a, int N, const QuReg& quregs, const std::vector& ctrl) + { + emulate_math([a,N](std::vector &res){for(auto& x: res) x = (x * a) % N;}, quregs, ctrl, true); } calc_type get_expectation_value(TermsDict const& td, std::vector const& ids){ run(); calc_type expectation = 0.; - auto current_state = vec_; + + StateVector current_state; // avoid costly memory reallocations + if( tmpBuff1_.capacity() >= vec_.size() ) + std::swap(tmpBuff1_, current_state); + current_state.resize(vec_.size()); +#pragma omp parallel for schedule(static) + for (std::size_t i = 0; i < vec_.size(); ++i) + current_state[i] = vec_[i]; + for (auto const& term : td){ auto const& coefficient = term.second; apply_term(term.first, ids, {}); @@ -260,17 +312,29 @@ class Simulator{ auto const a2 = std::real(vec_[i]); auto const b2 = std::imag(vec_[i]); delta += a1 * a2 - b1 * b2; + // reset vec_ + vec_[i] = current_state[i]; } expectation += coefficient * delta; - vec_ = current_state; } + std::swap(current_state, tmpBuff1_); return expectation; } void apply_qubit_operator(ComplexTermsDict const& td, std::vector const& ids){ run(); - auto new_state = StateVector(vec_.size(), 0.); - auto current_state = vec_; + StateVector new_state, current_state; // avoid costly memory reallocations + if( tmpBuff1_.capacity() >= vec_.size() ) + std::swap(tmpBuff1_, new_state); + if( tmpBuff2_.capacity() >= vec_.size() ) + std::swap(tmpBuff2_, current_state); + new_state.resize(vec_.size()); + current_state.resize(vec_.size()); +#pragma omp parallel for schedule(static) + for (std::size_t i = 0; i < vec_.size(); ++i){ + new_state[i] = 0; + current_state[i] = vec_[i]; + } for (auto const& term : td){ auto const& coefficient = term.second; apply_term(term.first, ids, {}); @@ -280,7 +344,9 @@ class Simulator{ vec_[i] = current_state[i]; } } - vec_ = std::move(new_state); + std::swap(vec_, new_state); + std::swap(tmpBuff1_, new_state); + std::swap(tmpBuff2_, current_state); } calc_type get_probability(std::vector const& bit_string, @@ -452,6 +518,8 @@ class Simulator{ #pragma omp parallel kernel(vec_, ids[4], ids[3], ids[2], ids[1], ids[0], m, ctrlmask); break; + default: + throw std::invalid_argument("Gates with more than 5 qubits are not supported!"); } fused_gates_ = Fusion(); @@ -500,6 +568,12 @@ class Simulator{ unsigned fusion_qubits_min_, fusion_qubits_max_; RndEngine rnd_eng_; std::function rng_; + + // large array buffers to avoid costly reallocations + static StateVector tmpBuff1_, tmpBuff2_; }; +Simulator::StateVector Simulator::tmpBuff1_; +Simulator::StateVector Simulator::tmpBuff2_; + #endif diff --git a/projectq/backends/_sim/_cppsim.cpp b/projectq/backends/_sim/_cppsim.cpp index 74498d4e2..cab68d0ee 100755 --- a/projectq/backends/_sim/_cppsim.cpp +++ b/projectq/backends/_sim/_cppsim.cpp @@ -50,6 +50,9 @@ PYBIND11_PLUGIN(_cppsim) { .def("measure_qubits", &Simulator::measure_qubits_return) .def("apply_controlled_gate", &Simulator::apply_controlled_gate) .def("emulate_math", &emulate_math_wrapper) + .def("emulate_math_addConstant", &Simulator::emulate_math_addConstant) + .def("emulate_math_addConstantModN", &Simulator::emulate_math_addConstantModN) + .def("emulate_math_multiplyByConstantModN", &Simulator::emulate_math_multiplyByConstantModN) .def("get_expectation_value", &Simulator::get_expectation_value) .def("apply_qubit_operator", &Simulator::apply_qubit_operator) .def("emulate_time_evolution", &Simulator::emulate_time_evolution) diff --git a/projectq/backends/_sim/_simulator.py b/projectq/backends/_sim/_simulator.py index 2218c3471..19e884d6d 100755 --- a/projectq/backends/_sim/_simulator.py +++ b/projectq/backends/_sim/_simulator.py @@ -33,10 +33,12 @@ TimeEvolution) from projectq.types import WeakQubitRef +FALLBACK_TO_PYSIM = False try: from ._cppsim import Simulator as SimulatorBackend except ImportError: from ._pysim import Simulator as SimulatorBackend + FALLBACK_TO_PYSIM = True class Simulator(BasicEngine): @@ -384,14 +386,34 @@ def _handle(self, cmd): ID = cmd.qubits[0][0].id self._simulator.deallocate_qubit(ID) elif isinstance(cmd.gate, BasicMathGate): + # improve performance by using C++ code for some commomn gates + from projectq.libs.math import (AddConstant, + AddConstantModN, + MultiplyByConstantModN) qubitids = [] for qr in cmd.qubits: qubitids.append([]) for qb in qr: qubitids[-1].append(qb.id) - math_fun = cmd.gate.get_math_function(cmd.qubits) - self._simulator.emulate_math(math_fun, qubitids, - [qb.id for qb in cmd.control_qubits]) + if FALLBACK_TO_PYSIM: + math_fun = cmd.gate.get_math_function(cmd.qubits) + self._simulator.emulate_math(math_fun, qubitids, + [qb.id for qb in cmd.control_qubits]) + else: + # individual code for different standard gates to make it faster! + if isinstance(cmd.gate, AddConstant): + self._simulator.emulate_math_addConstant(cmd.gate.a, qubitids, + [qb.id for qb in cmd.control_qubits]) + elif isinstance(cmd.gate, AddConstantModN): + self._simulator.emulate_math_addConstantModN(cmd.gate.a, cmd.gate.N, qubitids, + [qb.id for qb in cmd.control_qubits]) + elif isinstance(cmd.gate, MultiplyByConstantModN): + self._simulator.emulate_math_multiplyByConstantModN(cmd.gate.a, cmd.gate.N, qubitids, + [qb.id for qb in cmd.control_qubits]) + else: + math_fun = cmd.gate.get_math_function(cmd.qubits) + self._simulator.emulate_math(math_fun, qubitids, + [qb.id for qb in cmd.control_qubits]) elif isinstance(cmd.gate, TimeEvolution): op = [(list(term), coeff) for (term, coeff) in cmd.gate.hamiltonian.terms.items()] diff --git a/projectq/backends/_sim/_simulator_test.py b/projectq/backends/_sim/_simulator_test.py index 9e301be0d..9f7d298cb 100755 --- a/projectq/backends/_sim/_simulator_test.py +++ b/projectq/backends/_sim/_simulator_test.py @@ -683,3 +683,55 @@ def receive(command_list): qubit1[0].id: qubit0[0].id} assert (sim._convert_logical_to_mapped_qureg(qubit0 + qubit1) == qubit1 + qubit0) + + +def test_simulator_constant_math_emulation(): + if "cpp_simulator" not in get_available_simulators(): + pytest.skip("No C++ simulator") + return + + results = [[[1, 1, 0, 0, 0]], [[0, 1, 0, 0, 0]], [[0, 1, 1, 1, 0]]] + + import projectq.backends._sim._simulator as _sim + from projectq.backends._sim._pysim import Simulator as PySim + from projectq.backends._sim._cppsim import Simulator as CppSim + from projectq.libs.math import (AddConstant, AddConstantModN, + MultiplyByConstantModN) + + def gate_filter(eng, cmd): + g = cmd.gate + if isinstance(g, BasicMathGate): + return False + return eng.next_engine.is_available(cmd) + + def run_simulation(sim): + eng = MainEngine(sim, []) + quint = eng.allocate_qureg(5) + AddConstant(3) | quint + All(Measure) | quint + eng.flush() + results[0].append([int(qb) for qb in quint]) + + AddConstantModN(4, 5) | quint + All(Measure) | quint + eng.flush() + results[1].append([int(qb) for qb in quint]) + + MultiplyByConstantModN(15, 16) | quint + All(Measure) | quint + eng.flush() + results[2].append([int(qb) for qb in quint]) + + cppsim = Simulator(gate_fusion=False) + cppsim._simulator = CppSim(1) + run_simulation(cppsim) + + _sim.FALLBACK_TO_PYSIM = True + pysim = Simulator() + pysim._simulator = PySim(1) + # run_simulation(pysim) + + for result in results: + ref = result[0] + for res in result[1:]: + assert ref == res diff --git a/projectq/ops/_gates.py b/projectq/ops/_gates.py index 8967b1da5..be2240d00 100755 --- a/projectq/ops/_gates.py +++ b/projectq/ops/_gates.py @@ -157,10 +157,9 @@ def __str__(self): SqrtX = SqrtXGate() -class SwapGate(SelfInverseGate, BasicMathGate): +class SwapGate(SelfInverseGate): """ Swap gate class (swaps 2 qubits) """ def __init__(self): - BasicMathGate.__init__(self, lambda x, y: (y, x)) SelfInverseGate.__init__(self) self.interchangeable_qubit_indices = [[0, 1]] diff --git a/setup.py b/setup.py index 604f006af..a0254b0e1 100755 --- a/setup.py +++ b/setup.py @@ -124,6 +124,7 @@ def build_extensions(self): opts.append('/arch:AVX') else: opts.append('-march=native') + opts.append('-ffast-math') opts.append(openmp) if ct == 'unix': From 44ba35baf84b08741c51c5d96460df4d0a802d31 Mon Sep 17 00:00:00 2001 From: Nguyen Damien Date: Fri, 23 Aug 2019 21:12:36 +0200 Subject: [PATCH 6/7] Update constant math documentation to include list of pre-conditions (#331) --- projectq/libs/math/_gates.py | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/projectq/libs/math/_gates.py b/projectq/libs/math/_gates.py index fe1df6784..ef5cade99 100755 --- a/projectq/libs/math/_gates.py +++ b/projectq/libs/math/_gates.py @@ -90,6 +90,14 @@ class AddConstantModN(BasicMathGate): qunum = eng.allocate_qureg(5) # 5-qubit number X | qunum[1] # qunum is now equal to 2 AddConstantModN(3, 4) | qunum # qunum is now equal to 1 + + .. note:: + + Pre-conditions: + + * c < N + * c >= 0 + * The value stored in the quantum register must be lower than N """ def __init__(self, a, N): """ @@ -145,6 +153,14 @@ def SubConstantModN(a, N): qunum = eng.allocate_qureg(3) # 3-qubit number X | qunum[1] # qunum is now equal to 2 SubConstantModN(4,5) | qunum # qunum is now -2 = 6 = 1 (mod 5) + + .. note:: + + Pre-conditions: + + * c < N + * c >= 0 + * The value stored in the quantum register must be lower than N """ return AddConstantModN(N - a, N) @@ -162,6 +178,15 @@ class MultiplyByConstantModN(BasicMathGate): qunum = eng.allocate_qureg(5) # 5-qubit number X | qunum[2] # qunum is now equal to 4 MultiplyByConstantModN(3,5) | qunum # qunum is now 2. + + .. note:: + + Pre-conditions: + + * c < N + * c >= 0 + * gcd(c, N) == 1 + * The value stored in the quantum register must be lower than N """ def __init__(self, a, N): """ From d34941c39e4968f55621ea4800bd15a41cc96b8e Mon Sep 17 00:00:00 2001 From: cm Date: Fri, 25 Oct 2019 22:59:11 +0800 Subject: [PATCH 7/7] Add files via upload --- examples/bernstein_vazirani_tutorial.ipynb | 302 +++++++++++++++++++++ 1 file changed, 302 insertions(+) create mode 100644 examples/bernstein_vazirani_tutorial.ipynb diff --git a/examples/bernstein_vazirani_tutorial.ipynb b/examples/bernstein_vazirani_tutorial.ipynb new file mode 100644 index 000000000..8d93c0339 --- /dev/null +++ b/examples/bernstein_vazirani_tutorial.ipynb @@ -0,0 +1,302 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Bernstein-Varzirani Algorithm\n", + "The Bernstein-Varzirani algorithm solves the following problem:\n", + "1. Given an n-bit secret number x\n", + "2. There is an oracle holding the secret key that one can query with a number y, and the oracle will return $x \\bullet y$, the bit-wise dot product modulo 2 of the two numbers. For example, given a 4-bit secret number $x=13$ with binary representation 1101 and $y=15$ with binary representation 1111, then $x \\bullet y =$ 3 mod 2, which is 1\n", + "3. The question is: how many oracle queries does one need to figure out x?\n", + "\n", + "With the classical approach, if we query the oracle with 1, 2(10), 4(100) and 8(1000), then we can figure out each bit of the secret key and thus we need four queries. For a general n-bit number, we need n queries. Surprisingly, in the quantum computing paradigm, the Bernstein-Varzirani algorithm needs only one measurement, regardless of the size of the secret number" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import projectq\n", + "from projectq.backends import CircuitDrawer\n", + "from projectq.ops import H, Z, All, Measure, Barrier\n", + "\n", + "from pathlib import Path" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "n = 3 #No. of qubits" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The classical approach" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def query_oracle(y):\n", + " \"\"\"\n", + " The oracle returns the bit-wise dot product modulo 2 \n", + " of the secret key and y\n", + " \"\"\"\n", + " x = 6 #secret key\n", + " s = bin(x & y).count('1') % 2\n", + " \n", + " return s" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Secret key is 6 after 3 tries\n" + ] + } + ], + "source": [ + "b = ''.join([str(query_oracle(2**(k))) for k in reversed(range(n))])\n", + "print(\"Secret key is {} after {} tries\".format(int(b, 2), n))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The quantum computing approach\n", + "There are five steps in the Bernstein-Varzirani circuit construction:\n", + "1. Initialize all qubits to 0 (This is the default value anyway)\n", + "2. Apply Hadamard transform to all qubits\n", + "3. Implement the oracle \n", + "4. Apply Hadamard transform to all qubits\n", + "5. Measure the circuit\n", + "\n", + "The workings of the Bernstein-Varzirani algorithm can only be understood by working through the mathematics. After step 2, the application of the Hadamard transform to all qubits results in\n", + "$$\\begin{align} H^{\\otimes n} |0\\rangle = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} | y \\rangle \\end{align}$$\n", + "where $N = 2^n$\n", + "\n", + "The oracle seeks to implement the following operation in step 3:\n", + "$$\\begin{align} \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} | y \\rangle \\rightarrow \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle \\end{align}$$\n", + "The oracle effectively multiplies each superposition state $|y \\rangle$ with -1 raised to the power of $x \\bullet y$, the bit-wise inner product modulo 2. But why does the oracle implement this operation?\n", + "\n", + "Using the identity (proof given below)\n", + "$$\\begin{align} H^{\\otimes n} |x\\rangle = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle \\end{align}$$\n", + "If we further apply another Hadamard transform to all qubits in step 4, then $$ H^{\\otimes n} H^{\\otimes n} |x\\rangle = H^{\\otimes n} \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle $$\n", + "i.e. $$ |x\\rangle = H^{\\otimes n} \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle $$\n", + "since the Hadamard transform is symmetric and unitary (meaning that the transform is its own inverse). Thus measuring the circuit in step 5 will give the secret key.\n", + "\n", + "Summarizing,\n", + "$$ H^{\\otimes n} |0\\rangle ^{\\otimes n} \\rightarrow \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} | y \\rangle \\rightarrow \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle \\rightarrow H^{\\otimes n} \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle = |x\\rangle$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The oracle implementation\n", + "First note that in implementing\n", + "$$ \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{(x_{n-1} y_{n-1} \\oplus \\ldots \\oplus x_0 y_0)} | y_{n-1}\\ldots y_0\\rangle $$\n", + "$$ = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x_{n-1} y_{n-1}} (-1)^{x_{n-2} y_{n-2}} \\ldots (-1)^{x_0 y_0}) | y_{n-1}\\ldots y_0\\rangle $$\n", + "the individual product ${x_{n-1} y_{n-1}, \\ldots, x_0 y_0}$ only needs to be evaluated for bits in x that are 1. \n", + "\n", + "For the bits in x that are 1, we only introduce a minus sign for the bits in y that are 1 i.e.\n", + "$ |0 \\rangle \\rightarrow |0 \\rangle$ and $ |1 \\rangle \\rightarrow -|1 \\rangle$. This is precisely what the Z gate does. Thus the oracle can be implemented as: if a bit in x is 1, apply the Z gate to the corresponding bit in y." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "def make_bv_circuit(engine, n):\n", + " circuit = engine.allocate_qureg(n)\n", + " All(H) | circuit\n", + " Barrier | circuit\n", + " \n", + " # Oracle\n", + " x = 6 # secret key\n", + " for i in range(n):\n", + " if x & 1: # Only apply Z if the current bit of the secret key x is 1\n", + " Z | circuit[i]\n", + " x >>= 1 # Move the next bit to the 1 position\n", + " Barrier | circuit\n", + " \n", + " All(H) | circuit\n", + " All(Measure) | circuit\n", + " engine.flush()\n", + "\n", + " return circuit" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Draw circuit diagram" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "drawing_engine = CircuitDrawer()\n", + "main_engine = projectq.MainEngine(drawing_engine)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "circuit = make_bv_circuit(main_engine, n)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR0AAACBCAQAAABXsKZ8AAAAAmJLR0QA/4ePzL8AAAAJcEhZcwAAAEgAAABIAEbJaz4AAAAHdElNRQfjChIULw91KGppAAAKh0lEQVR42u2dX6gcVx3HP3vN9ZKKKXNTYyAS60QF9aVhEi1iHwJzQSSIPuxFKLUUYbcIYlOqOyhUFJFdERRE5a4++OLL3Uf/vNxRqgVBegdpDa1B3BjaUDXJHdqYFk2T9WFmZ2dnZ2Z3zszendn7+zwkOztz5vzOme+cOX/uzrc2oIYiK3yMD2ZO9Tp/YY+9DLk8xAOsZs7lIle4qlq0AljjQaW4X+QlbsTuW+WTrGeO42Uuc52bwBE+o36xJzmUK/X7Bz/PmqT2JP/mtUxJTg5+kDmXr+erlkJQibuVsnOFdYX6vsBt/0ZdZX3w/aIKV3typahTzcxqkcoXppKvcUjmyP5LR1gOVjzpmIkHmLOeSVhmmjHfTZNO2rNXOCDYGDHfRp+Epi+jHo6fqk5vtgz69P1PGkZkSyie/avvHlsx3473dUxaWFhYbPn5d6nPnoWLxYavOXBphraE4tmf+k6SwHir0wpaGJsWm4AL6IGgU9HRcXD8Z5y31ZEn3tzYn/p2cWjE7gm3OhpmIBI7kFovS7sjLBvJYgxLJ/qI1AHoySjr4NJnKINJ0iaMhg8qB0O6LMtPH5s+OjouDhomBt2UB+Asc41dWrED+wQ2QsEI86eI+u7ioGMGvZo6fWx+yKWc0umjEzf3U+MjcYfvBJ8sOpO7T3CGd2dYxUrIZQonOMtJrqtUY0GsKsZ9hmO8nnDG++O+nlLfH+AO69wCjnA0Lr1Nj/pEV1ingcMvsDASuslh6biRfU5wlj5gT6Rd4bhC5Vxll0tcm/l41VyeW/jK+X1Kce8mrpyvcUzhjH9jlxd4DTjKyehOFwsjdtbGmwjU2cKmSSNmrigsHQcXzf9s0A+k1Jh1UlCoFi4W7eCSRxlOBJqYWEyOosanBLvBfj0kFyOmxREqT7pwxicC23QnHkrj0ulgYOKtXFn+d3URznLSSRGOixPp3LaxIuJZiaTYwKWBHRpRmbM+rmwsukCT7sSWUDz56tuiniicuIlAjXbQmnhMjrCcsTkcDW3WMZ+JSTthSyiePPVtpy6SeqtPUTQaWKE8pv2pV0MeV8uHi5261pW0+GBAqF2ZJh1TnjfLR/oiaR8SH2XtkBw86SS3LCKcJcRI6eUQrMXH0wj6L15fJ1k6MqOzhNRz7B31kOTP2gVFRDqCIvl+pXOx9kTm3zcucl2p6rxcu5D5iv019Plm7anCfgV3M590rvFnBplTHeJuQeEfLAb0ua2Q7u3+NbrDZaX0CeSRzl32uKmU8s3iCnCAuM2NiYWkWRjwBgD/Ybc8vzl/o8CKEaYxULxRh7yV8HdAikg3WVBEpCMoItIRFBHpCIqIdARFRDqCIiIdQRGRjqCISEdQRKQjKCLSERQR6QiKiHQERUQ6giIiHUERkY6giEhHUESkIygi0hEUEekIiuT7s/ZDvE0p3f8UfoIzr/JXqQSrSrf6gNsMgFqxXmR5pFPj7D4YOM6TapUgr4HjO8tj4FjbJwPH+VGtEoiBY+UNHKtVAjFwFEqGGDgKsyAGjoIS5TZwjB7kKP0wvwD6kRe2Tol+gXGLgeMQjRY7/kvudHbYSn3hXSId1nM54nTZoOefwWUj8q7gecWtFrUYOHrl7mGEXu/dY5ftkF/PvuGwHVTVJlqC9UYZ4hYDx6SIe5iJBnAptNhTSRagB8W16dHIei7FuPNGPT+qaOCo4y6ivzOsBJcmhsr75hcU93yIf2+7R1kNHA3qNP1L8I6O/6GFhk0Xh21semxhRLY1LPqjC+590Oii0cBgE5M6zbTSDBvnDn1288Stcd4KoiYSpxsbtX+PPs2bftQOZijuL823yvs4OGgYuDjoGBipD8CyGTia/r8mZ4ILvOnSBjpssoOJzilsWnRwMCLbDbZGyzwdXvGNnLY5RZ+/Y9KhNcuN4PgH5oj729wzipqJOBOj/gr/CEVth+I+za3C69un64tleNPVcXD4Kc8XYuA42XbX+HDc4coGjvcHn2yu8Aw/49dAjcNem2P4ox0d77EyXH8c39aGuXyWh/m8f0kP8y3+AAww2WQ0j5VoDtlET660aAni4z7Ce8JRR+OMjRpO82m+GIl6GPcqH1Wo70/xANf4L3A4/m2z8QaOBgZNflyIgePkKHWFhzlPVpINHE1G7423ucMWf6THCsfbD3YZv6fMmKSRXD4EvOpv/ZJXsIHvTpQg1hyyg8O2fzm7kxUXLUF83M/yWPcn0ZbATNkC4Djwr0jUw7jVDBx/w58CA8dHojunGTie9g0c45yzymvg2GfoxXXuDFuY2HH3VDp25P8ZM+5QD8ZZmR8Cw7i/wDd0tah/73t/7oMnUIdWYjd4ZOBoxPr1ldfAUWPYwX+sHbo/Z76UvyM8OMgw8m1C6D7MLB0vbp0fsa0Q9QvAe1WiVsGinphFeCJQi/HgK6+Bo0YLlw7wvmFPwAt95mAu8WzwrGnMPsPbw/bHRF4FZhxnD+PWwevUZoz6Krs8mj1qFSzqiQOBqIFjnHgORVJsYNDADo1BzFlbXBvbtxD0ThHemopOgzqjzriJwwYu8Ft+1TnvXQKTJvewCWwGdu5eLsNtGxvowime4XG+TBsbgz4a28C2f3gKTTRcLLw1ImfapFZS3Dbf47wdRG1MxBkbNcDTfC6I2vEPiI07V32nCifZwHGsVzQAUma+vCqPZ4VHBpnhq3yCd81QtlAue4OdwV62XL7GxzmBN0OupZ17biVY4/GsUfvnPpoQ9RqPKkT7FA9xLwBHeWL47c6gnZJmb9CK/X439D3frISBY/r1T8XFXtjc7q2yRu3Sm2LgGN9uiYHjgaeTusDSJ7l/LgaOBxwttTVMb5HqQbsjBo4HkFaOvaNhg/xZu6CISEdQ5KAbOFarBEtk4Hi98gaOVSrBQNGAcWTg2Oet4sLJZ+B4Q9HXrSwGjtUqwW2uK/3SfWTg6JTlN+fLYOBYpRKIgaOwHIh0BEVEOoIiIh1BEZGOoIhIR1BEpCMoItIRFBHpCIqIdARFRDqCIiIdQRGRjqCISEdQRKQjKCLSERQR6QiKiHQERUQ6giIiHUERMXCsUgnEwJEbJXkpdbVKIAaOYuCoiBg4Vsr+sPolEANHoWSIgaMwC3GuYGLgKEzFjn1BXFkMHEtDRgPHUkR60A0cS0JmA8cFIgaOpSKjgeMCEQPHkpHTwHHpqKKB44LIbeC4VFTRwHFh5DRwrDAudsjA0cBAT3xvO5TPwLEkZDZwXCBF1HcPGw0zZOBo0+Myz6e0u2UycIwn0WQxlROc5STXVc+dycAxnlXFuM9wLOHFbashq8gQRRg4mhPDAROTJt9Jcbopk4FjPAkmi1NzeY4rU19TW4yBYzxr3KcU9y4v+S58k2cs3MARrMRRpI3BOc5V0MBxYeQ0cKwUaa5Y1TVwXBg5DRwrRJpwqmvguDByGjhWiHQfvuoYOJaGjAaOC2SeBo7WxDkmDRwnR1jO2ByOhjZrm22GBnfRreqg8jr0xZCnvu3URVIXYvZq1Mfmlith4CgUi6qBo4krBo4Hm/RFUjFwFBLRUhd1xcBRSEQMHIWFItIRFKmCgePF2gXu5XCmNP8sQQku1i6gsZbhzHd5NXV/qQwc/w9+31fX9CH7ZQAAACV0RVh0ZGF0ZTpjcmVhdGUAMjAxOS0xMC0xOFQyMDo0NzoxNSswODowMK2N2cUAAAAldEVYdGRhdGU6bW9kaWZ5ADIwMTktMTAtMThUMjA6NDc6MTUrMDg6MDDc0GF5AAAAFHRFWHRwZGY6VmVyc2lvbgBQREYtMS41IAVcCzkAAABKdEVYdHNpZ25hdHVyZQA2YTczNWVhMDZhYjBhZWQxNGFkNzdlYzI5MDM3MDg0OWE1ZGM5MDFkYmIwZjc0ODE2MThmYzIyYzJlOTE5OGM1OseklwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "p = Path('diagram')\n", + "if not p.exists(): #if the diagram directory doesn't exist, create it\n", + " p.mkdir()\n", + "with open('diagram/bv.tex', 'w') as f:\n", + " latex = drawing_engine.get_latex() #get circuit diagram as latex\n", + " f.write(latex) \n", + "#Change the pdf scale to 1.8 from 0.8 to have better visual effect\n", + "!sed -i 's@tikzpicture\\}\\[scale=0.8@tikzpicture\\}\\[scale=1.8@g' diagram/bv.tex\n", + "!cd diagram; pdflatex bv.tex > /dev/null #convert tex to latex, piping to /dev/null to silent output \n", + "\n", + "#Wand package needed to convert pdf to image\n", + "from wand.image import Image as WImage\n", + "img = WImage(filename='diagram/bv.pdf')\n", + "img" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run circuit" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Measurement outcome: 110\n", + "Secret key: 6\n" + ] + } + ], + "source": [ + "main_engine = projectq.MainEngine()\n", + "circuit = make_bv_circuit(main_engine, n)\n", + "measurement = ''.join(str(int(c)) for c in reversed(circuit))\n", + "print('Measurement outcome:', measurement)\n", + "print('Secret key:', int(measurement, 2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Proof of the key identity" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\\begin{align} H^{\\otimes n} |x\\rangle = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} (-1)^{x \\bullet y} | y \\rangle \\end{align}$$ where $ N = 2^n $ and the $ \\bullet$ operator denotes bitwise dot product modulo 2 i.e. $ x \\bullet y = x_0y_0 \\oplus x_1y_1 \\oplus \\ldots x_{N-1}y_{N-1}$ where $\\oplus$ denotes the sum modulo 2 operation or the XOR operation. We prove this identity by induction.\n", + "\n", + "Base case, n = 1:\n", + "$$ H |x_0\\rangle = \\frac{1}{\\sqrt{2}} (|0\\rangle + (-1)^{x_0}|1\\rangle) = \\frac{1}{\\sqrt{2}} \\sum_{y=0}^{1} (-1)^{x_0y} | y \\rangle$$\n", + "\n", + "Assuming n-1 (n > 1) is true i.e. $$ H^{\\otimes (n-1)} |x_{n-2} \\ldots x_0\\rangle = \\frac{1}{\\sqrt{2^{n-1}}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{(x_{n-2} \\ldots x_0) \\bullet (y_{n-2} \\ldots y_0)} | y \\rangle$$ \n", + "\n", + "$$ = \\frac{1}{\\sqrt{2^{n-1}}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{(x_{n-2} y_{n-2} \\oplus \\ldots \\oplus x_0 y_0)} | y \\rangle$$ \n", + "\n", + "where $y = y_{n-2} \\ldots y_0$ is the binary representation\n", + "\n", + "\n", + "Then $$ H^{\\otimes n} |x_{n-1} \\ldots x_0 \\rangle = \\frac{1}{\\sqrt{2}} (|0\\rangle + (-1)^{x_{n-1}}|1\\rangle) \\times \\frac{1}{\\sqrt{2^{n-1}}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{(x_{n-2} y_{n-2} \\oplus \\ldots \\oplus x_0 y_0)} | y \\rangle $$\n", + "\n", + "$$ \\begin{align} = \\frac{1}{\\sqrt{2^n}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{(x_{n-2} y_{n-2} \\oplus \\ldots \\oplus x_0 y_0)} |0\\rangle| y \\rangle + \\frac{1}{\\sqrt{2^n}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{x_{n-1}} (-1)^{(x_{n-2} y_{n-2} \\oplus \\ldots \\oplus x_0 y_0)} |1\\rangle| y \\rangle \\end{align} $$\n", + "\n", + "$$ \\begin{align} = \\frac{1}{\\sqrt{2^n}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{(x_{n-1}0 \\oplus x_{n-2}y_{n-2} \\oplus \\ldots \\oplus x_0y_0)} | y \\rangle + \\frac{1}{\\sqrt{2^n}} \\sum_{y=2^{n-1}}^{2^{n} - 1} (-1)^{(x_{n-1}1 \\oplus x_{n-2}y_{n-2}\\oplus \\ldots \\oplus x_0y_0)} | y \\rangle \\end{align}$$\n", + "\n", + "$$ = \\frac{1}{\\sqrt{2^n}} \\sum_{y=0}^{2^{n-1} - 1} (-1)^{(x_{n-1}x_{n-2} \\ldots x_0) \\bullet (0y_{n-2} \\ldots y_0)} | y \\rangle + \\frac{1}{\\sqrt{2^n}} \\sum_{y=2^{n-1}}^{2^{n} - 1} (-1)^{(x_{n-1} \\ldots x_0) \\bullet (1y_{{n}-2} \\ldots y_0)} | y \\rangle $$\n", + "\n", + "$$ = \\frac{1}{\\sqrt{2^{n}}} \\sum_{y=0}^{2^{n} - 1} (-1)^{(x_{n-1} \\ldots x_0) \\bullet (y_{n-1} \\ldots y_0)} | y \\rangle $$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "projectq", + "language": "python", + "name": "projectq" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.8" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}