diff --git a/examples/linq/1.the_basics.ipynb b/examples/linq/1.the_basics.ipynb index b6c0457e6..ce9958816 100755 --- a/examples/linq/1.the_basics.ipynb +++ b/examples/linq/1.the_basics.ipynb @@ -329,9 +329,10 @@ } ], "source": [ - "from tangelo.linq import SUPPORTED_GATES\n", + "from tangelo.linq import get_supported_gates\n", "\n", - "for backend, gates in SUPPORTED_GATES.items():\n", + "supported_gates = get_supported_gates()\n", + "for backend, gates in supported_gates.items():\n", " print(f'{backend} : {gates}')" ] }, @@ -844,7 +845,7 @@ "# Option2: Assume quantum circuit simulation was performed separately (by this package or different services)\n", "freqs, _ = sim.simulate(c) # circuit c must be consistent with operator, regarding measurement along non-z axes\n", "term, coef = tuple(op.terms.items())[0] # This yields ((0, 'Z'),), 1.0\n", - "expval2 = coef * Simulator.get_expectation_value_from_frequencies_oneterm(term, freqs)\n", + "expval2 = coef * sim.get_expectation_value_from_frequencies_oneterm(term, freqs)\n", "print(expval2)" ] }, diff --git a/examples/linq/2.qpu_connection.ipynb b/examples/linq/2.qpu_connection.ipynb index cf0b9bc2a..3adabee1a 100755 --- a/examples/linq/2.qpu_connection.ipynb +++ b/examples/linq/2.qpu_connection.ipynb @@ -767,7 +767,7 @@ } ], "source": [ - "from tangelo.linq import Simulator\n", + "from tangelo.linq import get_expectation_value_from_frequencies_oneterm\n", "from tangelo.linq.helpers.circuits.measurement_basis import pauli_string_to_of\n", "\n", "# Rescale counts to obtain frequency histogram\n", @@ -775,16 +775,16 @@ "print(f\"Frequencies: {freqs}\")\n", "\n", "# Compute expectation value for ZZ operator\n", - "exp_ZZ = Simulator.get_expectation_value_from_frequencies_oneterm(pauli_string_to_of(\"ZZ\"), freqs)\n", + "exp_ZZ = get_expectation_value_from_frequencies_oneterm(pauli_string_to_of(\"ZZ\"), freqs)\n", "print(f\"Expectation value for ZZ operator {exp_ZZ}\")" ] } ], "metadata": { "kernelspec": { - "display_name": "main_release0.3.1", + "display_name": "Python 3.9.7 64-bit ('qsdkmain': venv)", "language": "python", - "name": "main_release0.3.1" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -796,7 +796,12 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.10" + "version": "3.9.10" + }, + "vscode": { + "interpreter": { + "hash": "95050af2697fca56ed7491a4fb0b04c1282c0de0a7e0a7cacd318a8297b0b1d8" + } } }, "nbformat": 4, diff --git a/tangelo/algorithms/variational/tests/test_vqe_solver.py b/tangelo/algorithms/variational/tests/test_vqe_solver.py index ef50a5667..65fa0e847 100644 --- a/tangelo/algorithms/variational/tests/test_vqe_solver.py +++ b/tangelo/algorithms/variational/tests/test_vqe_solver.py @@ -16,6 +16,8 @@ import numpy as np from tangelo.linq import Simulator +from tangelo.helpers.utils import installed_backends +from tangelo.linq.target import QiskitSimulator from tangelo.algorithms import BuiltInAnsatze, VQESolver from tangelo.molecule_library import mol_H2_sto3g, mol_H4_sto3g, mol_H4_cation_sto3g, mol_NaH_sto3g, mol_H4_sto3g_symm from tangelo.toolboxes.ansatz_generator.uccsd import UCCSD @@ -197,9 +199,10 @@ def test_simulate_vsqs_h2(self): energy = vqe_solver.simulate() self.assertAlmostEqual(energy, -1.137270, delta=1e-4) + @unittest.skipIf("qiskit" not in installed_backends, "Test Skipped: Backend not available \n") def test_simulate_h2_qiskit(self): """Run VQE on H2 molecule, with UCCSD ansatz, JW qubit mapping, initial - parameters, exact qiskit simulator. + parameters, exact qiskit simulator. Both string and class input. """ backend_options = {"target": "qiskit", "n_shots": None, "noise_model": None} @@ -212,6 +215,16 @@ def test_simulate_h2_qiskit(self): self.assertAlmostEqual(energy, -1.13727042117, delta=1e-6) + backend_options = {"target": QiskitSimulator, "n_shots": None, "noise_model": None} + vqe_options = {"molecule": mol_H2_sto3g, "ansatz": BuiltInAnsatze.UCCSD, "qubit_mapping": "jw", + "initial_var_params": [6.28531447e-06, 5.65431626e-02], "verbose": True, + "backend_options": backend_options} + vqe_solver = VQESolver(vqe_options) + vqe_solver.build() + energy = vqe_solver.simulate() + + self.assertAlmostEqual(energy, -1.13727042117, delta=1e-6) + def test_simulate_h4(self): """Run VQE on H4 molecule, with UCCSD ansatz, JW qubit mapping, initial parameters, exact simulator. diff --git a/tangelo/helpers/utils.py b/tangelo/helpers/utils.py index f2191134b..65cc27103 100644 --- a/tangelo/helpers/utils.py +++ b/tangelo/helpers/utils.py @@ -19,6 +19,7 @@ import os import sys +from importlib import util class HiddenPrints: @@ -34,13 +35,9 @@ def __exit__(self, exc_type, exc_val, exc_tb): def is_package_installed(package_name): - try: - exec(f"import {package_name}") - # DEBUG print(f'{package_name:16s} :: found') - return True - except ModuleNotFoundError: - # DEBUG print(f'{package_name:16s} :: not found') - return False + """Check if module is installed without importing.""" + spam_spec = util.find_spec(package_name) + return spam_spec is not None # List all backends and statevector ones supported by Simulator class diff --git a/tangelo/linq/__init__.py b/tangelo/linq/__init__.py index 78cfd5380..e609cf0b5 100644 --- a/tangelo/linq/__init__.py +++ b/tangelo/linq/__init__.py @@ -15,4 +15,5 @@ from .gate import * from .circuit import Circuit, stack, remove_small_rotations, remove_redundant_gates from .translator import * +from .simulator_base import get_expectation_value_from_frequencies_oneterm from .simulator import Simulator, backend_info diff --git a/tangelo/linq/simulator.py b/tangelo/linq/simulator.py index 69171d5fd..466ebad6d 100644 --- a/tangelo/linq/simulator.py +++ b/tangelo/linq/simulator.py @@ -12,529 +12,43 @@ # See the License for the specific language governing permissions and # limitations under the License. -"""Simulator class, wrapping around the various simulators and abstracting their -differences from the user. Able to run noiseless and noisy simulations, -leveraging the capabilities of different backends, quantum or classical. - -If the user provides a noise model, then a noisy simulation is run with n_shots -shots. If the user only provides n_shots, a noiseless simulation is run, drawing -the desired amount of shots. If the target backend has access to the statevector -representing the quantum state, we leverage any kind of emulation available to -reduce runtime (directly generating shot values from final statevector etc) If -the quantum circuit contains a MEASURE instruction, it is assumed to simulate a -mixed-state and the simulation will be carried by simulating individual shots -(e.g a number of shots is required). - -Some backends may only support a subset of the above. This information is -contained in a separate data-structure. +"""Factory function acting as a unified interface able to instantiate any simulator object +associated to a specific target (qulacs, qiskit, cirq, user-defined...). Some +built-in target simulators support features such as noisy simulation, the ability +to run noiseless simulation with shots, or particular emulation techniques. Target +backends can also be implemented using APIs to quantum devices. """ -import os -import math -from collections import Counter - -import numpy as np -from scipy import stats -from bitarray import bitarray -from openfermion.ops import QubitOperator +from typing import Union, Type +from tangelo.linq.simulator_base import SimulatorBase +from tangelo.linq.target import target_dict, backend_info from tangelo.helpers.utils import default_simulator -from tangelo.linq import Gate, Circuit -from tangelo.linq.helpers.circuits.measurement_basis import measurement_basis_gates -import tangelo.linq.translator as translator - - -# Data-structure showing what functionalities are supported by the backend, in this package -backend_info = dict() -backend_info["qiskit"] = {"statevector_available": True, "statevector_order": "msq_first", "noisy_simulation": True} -backend_info["qulacs"] = {"statevector_available": True, "statevector_order": "msq_first", "noisy_simulation": True} -backend_info["cirq"] = {"statevector_available": True, "statevector_order": "lsq_first", "noisy_simulation": True} -backend_info["qdk"] = {"statevector_available": False, "statevector_order": None, "noisy_simulation": False} - - -class Simulator: - - def __init__(self, target=default_simulator, n_shots=None, noise_model=None): - """Instantiate Simulator object. - - Args: - target (str): One of the available target backends (quantum or - classical). The default simulator is qulacs if found in user's - environment, otherwise cirq. - n_shots (int): Number of shots if using a shot-based simulator. - noise_model: A noise model object assumed to be in the format - expected from the target backend. - """ - self._source = "abstract" - self._target = target if target else default_simulator - self._current_state = None - self._noise_model = noise_model - - # Can be modified later by user as long as long as it retains the same type (ex: cannot change to/from None) - self.n_shots = n_shots - self.freq_threshold = 1e-10 - - # Set additional attributes related to the target backend chosen by the user - for k, v in backend_info[self._target].items(): - setattr(self, k, v) - - # Raise error if user attempts to pass a noise model to a backend not supporting noisy simulation - if self._noise_model and not self.noisy_simulation: - raise ValueError("Target backend does not support noise models.") - - # Raise error if the number of shots has not been passed for a noisy simulation or if statevector unavailable - if not self.n_shots and (not self.statevector_available or self._noise_model): - raise ValueError("A number of shots needs to be specified.") - - def simulate(self, source_circuit, return_statevector=False, initial_statevector=None): - """Perform state preparation corresponding to the input circuit on the - target backend, return the frequencies of the different observables, and - either the statevector or None depending on the availability of the - statevector and if return_statevector is set to True. For the - statevector backends supporting it, an initial statevector can be - provided to initialize the quantum state without simulating all the - equivalent gates. - - Args: - source_circuit: a circuit in the abstract format to be translated - for the target backend. - return_statevector(bool): option to return the statevector as well, - if available. - initial_statevector(list/array) : A valid statevector in the format - supported by the target backend. - - Returns: - dict: A dictionary mapping multi-qubit states to their corresponding - frequency. - numpy.array: The statevector, if available for the target backend - and requested by the user (if not, set to None). - """ - - if source_circuit.is_mixed_state and not self.n_shots: - raise ValueError("Circuit contains MEASURE instruction, and is assumed to prepare a mixed state." - "Please set the Simulator.n_shots attribute to an appropriate value.") - - if source_circuit.width == 0: - raise ValueError("Cannot simulate an empty circuit (e.g identity unitary) with unknown number of qubits.") - - # If the unitary is the identity (no gates) and no noise model, no need for simulation: - # return all-zero state or sample from statevector - if source_circuit.size == 0 and not self._noise_model: - if initial_statevector is not None: - statevector = initial_statevector - frequencies = self._statevector_to_frequencies(initial_statevector) - else: - frequencies = {'0'*source_circuit.width: 1.0} - statevector = np.zeros(2**source_circuit.width) - statevector[0] = 1.0 - return (frequencies, statevector) if return_statevector else (frequencies, None) - - if self._target == "qulacs": - import qulacs - - translated_circuit = translator.translate_qulacs(source_circuit, self._noise_model) - - # Initialize state on GPU if available and desired. Default to CPU otherwise. - if ('QuantumStateGpu' in dir(qulacs)) and (int(os.getenv("QULACS_USE_GPU", 0)) != 0): - state = qulacs.QuantumStateGpu(source_circuit.width) - else: - state = qulacs.QuantumState(source_circuit.width) - if initial_statevector is not None: - state.load(initial_statevector) - - if (source_circuit.is_mixed_state or self._noise_model): - samples = list() - for i in range(self.n_shots): - translated_circuit.update_quantum_state(state) - samples.append(state.sampling(1)[0]) - if initial_statevector is not None: - state.load(initial_statevector) - else: - state.set_zero_state() - python_statevector = None - elif self.n_shots is not None: - translated_circuit.update_quantum_state(state) - python_statevector = np.array(state.get_vector()) if return_statevector else None - samples = state.sampling(self.n_shots) - else: - translated_circuit.update_quantum_state(state) - self._current_state = state - python_statevector = state.get_vector() - frequencies = self._statevector_to_frequencies(python_statevector) - return (frequencies, np.array(python_statevector)) if return_statevector else (frequencies, None) - - frequencies = {self.__int_to_binstr(k, source_circuit.width): v / self.n_shots - for k, v in Counter(samples).items()} - return (frequencies, python_statevector) - - elif self._target == "qiskit": - import qiskit - - translated_circuit = translator.translate_qiskit(source_circuit) - - # If requested, set initial state - if initial_statevector is not None: - if self._noise_model: - raise ValueError("Cannot load an initial state if using a noise model, with Qiskit") - else: - n_qubits = int(math.log2(len(initial_statevector))) - initial_state_circuit = qiskit.QuantumCircuit(n_qubits, n_qubits) - initial_state_circuit.initialize(initial_statevector, list(range(n_qubits))) - translated_circuit = initial_state_circuit.compose(translated_circuit) - - # Drawing individual shots with the qasm simulator, for noisy simulation or simulating mixed states - if self._noise_model or source_circuit.is_mixed_state: - from tangelo.linq.noisy_simulation.noise_models import get_qiskit_noise_model - - meas_range = range(source_circuit.width) - translated_circuit.measure(meas_range, meas_range) - return_statevector = False - backend = qiskit.Aer.get_backend("aer_simulator") - - qiskit_noise_model = get_qiskit_noise_model(self._noise_model) if self._noise_model else None - opt_level = 0 if self._noise_model else None - - job_sim = qiskit.execute(translated_circuit, backend, noise_model=qiskit_noise_model, - shots=self.n_shots, basis_gates=None, optimization_level=opt_level) - sim_results = job_sim.result() - frequencies = {state[::-1]: count/self.n_shots for state, count in sim_results.get_counts(0).items()} - - # Noiseless simulation using the statevector simulator otherwise - else: - backend = qiskit.Aer.get_backend("aer_simulator", method='statevector') - translated_circuit = qiskit.transpile(translated_circuit, backend) - translated_circuit.save_statevector() - sim_results = backend.run(translated_circuit).result() - self._current_state = sim_results.get_statevector(translated_circuit) - frequencies = self._statevector_to_frequencies(self._current_state) - - return (frequencies, np.array(sim_results.get_statevector())) if return_statevector else (frequencies, None) - - elif self._target == "qdk": - - translated_circuit = translator.translate_qsharp(source_circuit) - with open('tmp_circuit.qs', 'w+') as f_out: - f_out.write(translated_circuit) - - # Compile, import and call Q# operation to compute frequencies. Only import qsharp module if qdk is running - # TODO: A try block to catch an exception at compile time, for Q#? Probably as an ImportError. - import qsharp - qsharp.reload() - from MyNamespace import EstimateFrequencies - frequencies_list = EstimateFrequencies.simulate(nQubits=source_circuit.width, nShots=self.n_shots) - print("Q# frequency estimation with {0} shots: \n {1}".format(self.n_shots, frequencies_list)) - - # Convert Q# output to frequency dictionary, apply threshold - frequencies = {bin(i).split('b')[-1]: frequencies_list[i] for i, freq in enumerate(frequencies_list)} - frequencies = {("0"*(source_circuit.width-len(k))+k)[::-1]: v for k, v in frequencies.items() - if v > self.freq_threshold} - return (frequencies, None) - - elif self._target == "cirq": - import cirq - - translated_circuit = translator.translate_cirq(source_circuit, self._noise_model) - - if source_circuit.is_mixed_state or self._noise_model: - # Only DensityMatrixSimulator handles noise well, can use Simulator but it is slower - cirq_simulator = cirq.DensityMatrixSimulator(dtype=np.complex128) - else: - cirq_simulator = cirq.Simulator(dtype=np.complex128) - - # If requested, set initial state - cirq_initial_statevector = initial_statevector if initial_statevector is not None else 0 - - # Calculate final density matrix and sample from that for noisy simulation or simulating mixed states - if self._noise_model or source_circuit.is_mixed_state: - # cirq.dephase_measurements changes measurement gates to Krauss operators so simulators - # can be called once and density matrix sampled repeatedly. - translated_circuit = cirq.dephase_measurements(translated_circuit) - sim = cirq_simulator.simulate(translated_circuit, initial_state=cirq_initial_statevector) - self._current_state = sim.final_density_matrix - indices = list(range(source_circuit.width)) - isamples = cirq.sample_density_matrix(sim.final_density_matrix, indices, repetitions=self.n_shots) - samples = [''.join([str(int(q))for q in isamples[i]]) for i in range(self.n_shots)] - - frequencies = {k: v / self.n_shots for k, v in Counter(samples).items()} - # Noiseless simulation using the statevector simulator otherwise - else: - job_sim = cirq_simulator.simulate(translated_circuit, initial_state=cirq_initial_statevector) - self._current_state = job_sim.final_state_vector - frequencies = self._statevector_to_frequencies(self._current_state) - - return (frequencies, np.array(self._current_state)) if return_statevector else (frequencies, None) - - def get_expectation_value(self, qubit_operator, state_prep_circuit, initial_statevector=None): - r"""Take as input a qubit operator H and a quantum circuit preparing a - state |\psi>. Return the expectation value <\psi | H | \psi>. - - In the case of a noiseless simulation, if the target backend exposes the - statevector then it is used directly to compute expectation values, or - draw samples if required. In the case of a noisy simulator, or if the - statevector is not available on the target backend, individual shots - must be run and the workflow is akin to what we would expect from an - actual QPU. - - Args: - qubit_operator(openfermion-style QubitOperator class): a qubit - operator. - state_prep_circuit: an abstract circuit used for state preparation. - - Returns: - complex: The expectation value of this operator with regards to the - state preparation. - """ - # Check if simulator supports statevector - if initial_statevector is not None and not self.statevector_available: - raise ValueError(f'Statevector not supported in {self._target}') - - # Check that qubit operator does not operate on qubits beyond circuit size. - # Keep track if coefficients are real or not - are_coefficients_real = True - for term, coef in qubit_operator.terms.items(): - if state_prep_circuit.width < len(term): - raise ValueError(f'Term {term} requires more qubits than the circuit contains ({state_prep_circuit.width})') - if type(coef) in {complex, np.complex64, np.complex128}: - are_coefficients_real = False - - # If the underlying operator is hermitian, expectation value is real and can be computed right away - if are_coefficients_real: - if self._noise_model or not self.statevector_available \ - or state_prep_circuit.is_mixed_state or state_prep_circuit.size == 0: - return self._get_expectation_value_from_frequencies(qubit_operator, state_prep_circuit, initial_statevector=initial_statevector) - elif self.statevector_available: - return self._get_expectation_value_from_statevector(qubit_operator, state_prep_circuit, initial_statevector=initial_statevector) - - # Else, separate the operator into 2 hermitian operators, use linearity and call this function twice - else: - qb_op_real, qb_op_imag = QubitOperator(), QubitOperator() - for term, coef in qubit_operator.terms.items(): - qb_op_real.terms[term], qb_op_imag.terms[term] = coef.real, coef.imag - qb_op_real.compress() - qb_op_imag.compress() - exp_real = self.get_expectation_value(qb_op_real, state_prep_circuit, initial_statevector=initial_statevector) - exp_imag = self.get_expectation_value(qb_op_imag, state_prep_circuit, initial_statevector=initial_statevector) - return exp_real if (exp_imag == 0.) else exp_real + 1.0j * exp_imag - - def _get_expectation_value_from_statevector(self, qubit_operator, state_prep_circuit, initial_statevector=None): - r"""Take as input a qubit operator H and a state preparation returning a - ket |\psi>. Return the expectation value <\psi | H | \psi>, computed - without drawing samples (statevector only). Users should not be calling - this function directly, please call "get_expectation_value" instead. - - Args: - qubit_operator(openfermion-style QubitOperator class): a qubit - operator. - state_prep_circuit: an abstract circuit used for state preparation - (only pure states). - - Returns: - complex: The expectation value of this operator with regards to the - state preparation. - """ - n_qubits = state_prep_circuit.width - - expectation_value = 0. - prepared_frequencies, prepared_state = self.simulate(state_prep_circuit, return_statevector=True, initial_statevector=initial_statevector) - - # Use fast built-in qulacs expectation value function if possible - if self._target == "qulacs" and not self.n_shots: - import qulacs - - # Note: This section previously used qulacs.quantum_operator.create_quantum_operator_from_openfermion_text but was changed - # due to a memory leak. We can re-evaluate the implementation if/when Issue #303 (https://github.com/qulacs/qulacs/issues/303) - # is fixed. - operator = qulacs.Observable(n_qubits) - for term, coef in qubit_operator.terms.items(): - pauli_string = "".join(f" {op} {qu}" for qu, op in term) - operator.add_operator(coef, pauli_string) - return operator.get_expectation_value(self._current_state).real - - # Use cirq built-in expectation_from_state_vector/epectation_from_density_matrix - # noise model would require - if self._target == "cirq" and not self.n_shots: - import cirq - - GATE_CIRQ = translator.get_cirq_gates() - qubit_labels = cirq.LineQubit.range(n_qubits) - qubit_map = {q: i for i, q in enumerate(qubit_labels)} - paulisum = 0.*cirq.PauliString(cirq.I(qubit_labels[0])) - for term, coef in qubit_operator.terms.items(): - pauli_list = [GATE_CIRQ[pauli](qubit_labels[index]) for index, pauli in term] - paulisum += cirq.PauliString(pauli_list, coefficient=coef) - if self._noise_model: - exp_value = paulisum.expectation_from_density_matrix(prepared_state, qubit_map) - else: - exp_value = paulisum.expectation_from_state_vector(prepared_state, qubit_map) - return np.real(exp_value) - - # Otherwise, use generic statevector expectation value - for term, coef in qubit_operator.terms.items(): - - if len(term) > n_qubits: # Cannot have a qubit index beyond circuit size - raise ValueError(f"Size of operator {qubit_operator} beyond circuit width ({n_qubits} qubits)") - elif not term: # Empty term: no simulation needed - expectation_value += coef - continue - - if not self.n_shots: - # Directly simulate and compute expectation value using statevector - pauli_circuit = Circuit([Gate(pauli, index) for index, pauli in term], n_qubits=n_qubits) - _, pauli_state = self.simulate(pauli_circuit, return_statevector=True, initial_statevector=prepared_state) - - delta = np.dot(pauli_state.real, prepared_state.real) + np.dot(pauli_state.imag, prepared_state.imag) - expectation_value += coef * delta - - else: - # Run simulation with statevector but compute expectation value with samples directly drawn from it - basis_circuit = Circuit(measurement_basis_gates(term), n_qubits=state_prep_circuit.width) - if basis_circuit.size > 0: - frequencies, _ = self.simulate(basis_circuit, initial_statevector=prepared_state) - else: - frequencies = prepared_frequencies - expectation_term = self.get_expectation_value_from_frequencies_oneterm(term, frequencies) - expectation_value += coef * expectation_term - - return expectation_value - - def _get_expectation_value_from_frequencies(self, qubit_operator, state_prep_circuit, initial_statevector=None): - r"""Take as input a qubit operator H and a state preparation returning a - ket |\psi>. Return the expectation value <\psi | H | \psi> computed - using the frequencies of observable states. - - Args: - qubit_operator(openfermion-style QubitOperator class): a qubit - operator. - state_prep_circuit: an abstract circuit used for state preparation. - - Returns: - complex: The expectation value of this operator with regards to the - state preparation. - """ - n_qubits = state_prep_circuit.width - if not self.statevector_available or state_prep_circuit.is_mixed_state or self._noise_model: - initial_circuit = state_prep_circuit - if initial_statevector is not None and not self.statevector_available: - raise ValueError(f'Backend {self._target} does not support statevectors') - else: - updated_statevector = initial_statevector - else: - initial_circuit = Circuit(n_qubits=n_qubits) - _, updated_statevector = self.simulate(state_prep_circuit, return_statevector=True, initial_statevector=initial_statevector) - - expectation_value = 0. - for term, coef in qubit_operator.terms.items(): - - if len(term) > n_qubits: - raise ValueError(f"Size of operator {qubit_operator} beyond circuit width ({n_qubits} qubits)") - elif not term: # Empty term: no simulation needed - expectation_value += coef - continue - - basis_circuit = Circuit(measurement_basis_gates(term)) - full_circuit = initial_circuit + basis_circuit if (basis_circuit.size > 0) else initial_circuit - frequencies, _ = self.simulate(full_circuit, initial_statevector=updated_statevector) - expectation_term = self.get_expectation_value_from_frequencies_oneterm(term, frequencies) - expectation_value += coef * expectation_term - - return expectation_value - - @staticmethod - def get_expectation_value_from_frequencies_oneterm(term, frequencies): - """Return the expectation value of a single-term qubit-operator, given - the result of a state-preparation. - - Args: - term(openfermion-style QubitOperator object): a qubit operator, with - only a single term. - frequencies(dict): histogram of frequencies of measurements (assumed - to be in lsq-first format). - - Returns: - complex: The expectation value of this operator with regards to the - state preparation. - """ - - if not frequencies.keys(): - return ValueError("Must pass a non-empty dictionary of frequencies.") - n_qubits = len(list(frequencies.keys())[0]) - - # Get term mask - mask = ["0"] * n_qubits - for index, op in term: - mask[index] = "1" - mask = "".join(mask) - - # Compute expectation value of the term - expectation_term = 0. - for basis_state, freq in frequencies.items(): - # Compute sample value using state_binstr and term mask, update term expectation value - sample = (-1) ** ((bitarray(mask) & bitarray(basis_state)).to01().count("1") % 2) - expectation_term += sample * freq - - return expectation_term - - def _statevector_to_frequencies(self, statevector): - """For a given statevector representing the quantum state of a qubit - register, returns a sparse histogram of the probabilities in the - least-significant-qubit (lsq) -first order. e.g the string '100' means - qubit 0 measured in basis state |1>, and qubit 1 & 2 both measured in - state |0>. - - Args: - statevector(list or ndarray(complex)): an iterable 1D data-structure - containing the amplitudes. - - Returns: - dict: A dictionary whose keys are bitstrings representing the - multi-qubit states with the least significant qubit first (e.g - '100' means qubit 0 in state |1>, and qubit 1 and 2 in state - |0>), and the associated value is the corresponding frequency. - Unless threshold=0., this dictionary will be sparse. - """ - - n_qubits = int(math.log2(len(statevector))) - frequencies = dict() - for i, amplitude in enumerate(statevector): - frequency = abs(amplitude)**2 - if (frequency - self.freq_threshold) >= 0.: - frequencies[self.__int_to_binstr(i, n_qubits)] = frequency - - # If n_shots, has been specified, then draw that amount of samples from the distribution - # and return empirical frequencies instead. Otherwise, return the exact frequencies - if not self.n_shots: - return frequencies - else: - xk, pk = [], [] - for k, v in frequencies.items(): - xk.append(int(k[::-1], 2)) - pk.append(frequencies[k]) - distr = stats.rv_discrete(name='distr', values=(np.array(xk), np.array(pk))) - - # Generate samples from distribution. Cut in chunks to ensure samples fit in memory, gradually accumulate - chunk_size = 10**7 - n_chunks = self.n_shots // chunk_size - freqs_shots = Counter() - - for i in range(n_chunks+1): - this_chunk = self.n_shots % chunk_size if i == n_chunks else chunk_size - samples = distr.rvs(size=this_chunk) - freqs_shots += Counter(samples) - freqs_shots = {self.__int_to_binstr_lsq(k, n_qubits): v / self.n_shots for k, v in freqs_shots.items()} - return freqs_shots - def __int_to_binstr(self, i, n_qubits): - """Convert an integer into a bit string of size n_qubits, in the order - specified for the state vector. - """ - bs = bin(i).split('b')[-1] - state_binstr = "0" * (n_qubits - len(bs)) + bs - return state_binstr[::-1] if (self.statevector_order == "msq_first") else state_binstr - def __int_to_binstr_lsq(self, i, n_qubits): - """Convert an integer into a bit string of size n_qubits, in the - least-significant qubit order. - """ - bs = bin(i).split('b')[-1] - state_binstr = "0" * (n_qubits - len(bs)) + bs - return state_binstr[::-1] +def Simulator(target: Union[None, str, Type[SimulatorBase]] = default_simulator, n_shots: Union[None, int] = None, + noise_model=None, **kwargs): + """Return requested target backend object. + + Args: + target (string or Type[SimulatorBase] or None): Supported string identifiers can be found in + target_dict (from tangelo.linq.target). Can also provide a subclass of SimulatorBase. + n_shots (int): Number of shots if using a shot-based simulator. + noise_model: A noise model object assumed to be in the format expected from the target backend. + kwargs: Other arguments that could be passed to a target. Examples are qubits_to_use for a QPU, transpiler + optimization level, error mitigation flags etc. + """ + + if target is None: + target = target_dict[default_simulator] + # If target is a string use target_dict to return built-in Target Simulators + elif isinstance(target, str): + try: + target = target_dict[target] + except KeyError: + raise ValueError(f"Error: backend {target} not supported. Available built-in options: {list(target_dict.keys())}") + # If subclass of SimulatorBase, use target + elif not issubclass(target, SimulatorBase): + raise TypeError(f"target must be a str or a subclass of SimulatorBase but received class {type(target).__name__}") + + return target(n_shots=n_shots, noise_model=noise_model, **kwargs) diff --git a/tangelo/linq/simulator_base.py b/tangelo/linq/simulator_base.py new file mode 100644 index 000000000..1d219a516 --- /dev/null +++ b/tangelo/linq/simulator_base.py @@ -0,0 +1,418 @@ +# Copyright 2021 Good Chemistry Company. +# +# 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. + +"""SimulatorBase class, base class to define various target simulators and abstracting their +differences from the user. Users can define their own target simulator as a child class of SimulatorBase. +Able to run noiseless and noisy simulations, leveraging the capabilities of different backends, +quantum or classical. + +If the user provides a noise model, then a noisy simulation is run with n_shots +shots. If the user only provides n_shots, a noiseless simulation is run, drawing +the desired amount of shots. If the target backend has access to the statevector +representing the quantum state, we leverage any kind of emulation available to +reduce runtime (directly generating shot values from final statevector etc). If +the quantum circuit contains a MEASURE instruction, it is assumed to simulate a +mixed-state and the simulation will be carried by simulating individual shots +(e.g a number of shots is required). + +Some backends may only support a subset of the above. This information is +contained in the static method backend_info defined by each subclass target. +""" + +import abc +import math +from collections import Counter + +import numpy as np +from scipy import stats +from bitarray import bitarray +from openfermion.ops import QubitOperator + +from tangelo.linq import Gate, Circuit +from tangelo.linq.helpers.circuits.measurement_basis import measurement_basis_gates + + +def get_expectation_value_from_frequencies_oneterm(term, frequencies): + """Return the expectation value of a single-term qubit-operator, given + the result of a state-preparation. + + Args: + term(openfermion-style QubitOperator object): a qubit operator, with + only a single term. + frequencies(dict): histogram of frequencies of measurements (assumed + to be in lsq-first format). + + Returns: + complex: The expectation value of this operator with regards to the + state preparation. + """ + + if not frequencies.keys(): + return ValueError("Must pass a non-empty dictionary of frequencies.") + n_qubits = len(list(frequencies.keys())[0]) + + # Get term mask + mask = ["0"] * n_qubits + for index, op in term: + mask[index] = "1" + mask = "".join(mask) + + # Compute expectation value of the term + expectation_term = 0. + for basis_state, freq in frequencies.items(): + # Compute sample value using state_binstr and term mask, update term expectation value + sample = (-1) ** ((bitarray(mask) & bitarray(basis_state)).to01().count("1") % 2) + expectation_term += sample * freq + + return expectation_term + + +class SimulatorBase(abc.ABC): + + def __init__(self, n_shots=None, noise_model=None): + """Instantiate Simulator object. + + Args: + n_shots (int): Number of shots if using a shot-based simulator. + noise_model: A noise model object assumed to be in the format + expected from the target backend. + """ + + self._current_state = None + self._noise_model = noise_model + + # Can be modified later by user as long as long as it retains the same type (ex: cannot change to/from None) + self.n_shots = n_shots + self.freq_threshold = 1e-10 + + # Set additional attributes related to the target backend chosen by the user + for k, v in self.backend_info().items(): + setattr(self, k, v) + + # Raise error if user attempts to pass a noise model to a backend not supporting noisy simulation + if self._noise_model and not self.noisy_simulation: + raise ValueError("Target backend does not support noise models.") + + # Raise error if the number of shots has not been passed for a noisy simulation or if statevector unavailable + if not self.n_shots and (not self.statevector_available or self._noise_model): + raise ValueError("A number of shots needs to be specified.") + + @abc.abstractmethod + def simulate_circuit(self): + """Perform state preparation corresponding to the input circuit on the + target backend, return the frequencies of the different observables, and + either the statevector or None depending on the availability of the + statevector and if return_statevector is set to True. For the + statevector backends supporting it, an initial statevector can be + provided to initialize the quantum state without simulating all the + equivalent gates. + + Args: + source_circuit: a circuit in the abstract format to be translated + for the target backend. + return_statevector(bool): option to return the statevector as well, + if available. + initial_statevector(list/array) : A valid statevector in the format + supported by the target backend. + + Returns: + dict: A dictionary mapping multi-qubit states to their corresponding + frequency. + numpy.array: The statevector, if available for the target backend + and requested by the user (if not, set to None). + """ + pass + + def simulate(self, source_circuit, return_statevector=False, initial_statevector=None): + """Perform state preparation corresponding to the input circuit on the + target backend, return the frequencies of the different observables, and + either the statevector or None depending on the availability of the + statevector and if return_statevector is set to True. For the + statevector backends supporting it, an initial statevector can be + provided to initialize the quantum state without simulating all the + equivalent gates. + + Args: + source_circuit: a circuit in the abstract format to be translated + for the target backend. + return_statevector(bool): option to return the statevector as well, + if available. + initial_statevector(list/array) : A valid statevector in the format + supported by the target backend. + + Returns: + dict: A dictionary mapping multi-qubit states to their corresponding + frequency. + numpy.array: The statevector, if available for the target backend + and requested by the user (if not, set to None). + """ + + if source_circuit.is_mixed_state and not self.n_shots: + raise ValueError("Circuit contains MEASURE instruction, and is assumed to prepare a mixed state." + "Please set the Simulator.n_shots attribute to an appropriate value.") + + if source_circuit.width == 0: + raise ValueError("Cannot simulate an empty circuit (e.g identity unitary) with unknown number of qubits.") + + # If the unitary is the identity (no gates) and no noise model, no need for simulation: + # return all-zero state or sample from statevector + if source_circuit.size == 0 and not self._noise_model: + if initial_statevector is not None: + statevector = initial_statevector + frequencies = self._statevector_to_frequencies(initial_statevector) + else: + frequencies = {'0'*source_circuit.width: 1.0} + statevector = np.zeros(2**source_circuit.width) + statevector[0] = 1.0 + return (frequencies, statevector) if return_statevector else (frequencies, None) + + return self.simulate_circuit(source_circuit, return_statevector=return_statevector, initial_statevector=initial_statevector) + + def get_expectation_value(self, qubit_operator, state_prep_circuit, initial_statevector=None): + r"""Take as input a qubit operator H and a quantum circuit preparing a + state |\psi>. Return the expectation value <\psi | H | \psi>. + + In the case of a noiseless simulation, if the target backend exposes the + statevector then it is used directly to compute expectation values, or + draw samples if required. In the case of a noisy simulator, or if the + statevector is not available on the target backend, individual shots + must be run and the workflow is akin to what we would expect from an + actual QPU. + + Args: + qubit_operator(openfermion-style QubitOperator class): a qubit + operator. + state_prep_circuit: an abstract circuit used for state preparation. + + Returns: + complex: The expectation value of this operator with regards to the + state preparation. + """ + # Check if simulator supports statevector + if initial_statevector is not None and not self.statevector_available: + raise ValueError(f'Statevector not supported in {self.__class__}') + + # Check that qubit operator does not operate on qubits beyond circuit size. + # Keep track if coefficients are real or not + are_coefficients_real = True + for term, coef in qubit_operator.terms.items(): + if state_prep_circuit.width < len(term): + raise ValueError(f'Term {term} requires more qubits than the circuit contains ({state_prep_circuit.width})') + if type(coef) in {complex, np.complex64, np.complex128}: + are_coefficients_real = False + + # If the underlying operator is hermitian, expectation value is real and can be computed right away + if are_coefficients_real: + if self._noise_model or not self.statevector_available \ + or state_prep_circuit.is_mixed_state or state_prep_circuit.size == 0: + return self._get_expectation_value_from_frequencies(qubit_operator, state_prep_circuit, initial_statevector=initial_statevector) + elif self.statevector_available: + return self._get_expectation_value_from_statevector(qubit_operator, state_prep_circuit, initial_statevector=initial_statevector) + + # Else, separate the operator into 2 hermitian operators, use linearity and call this function twice + else: + qb_op_real, qb_op_imag = QubitOperator(), QubitOperator() + for term, coef in qubit_operator.terms.items(): + qb_op_real.terms[term], qb_op_imag.terms[term] = coef.real, coef.imag + qb_op_real.compress() + qb_op_imag.compress() + exp_real = self.get_expectation_value(qb_op_real, state_prep_circuit, initial_statevector=initial_statevector) + exp_imag = self.get_expectation_value(qb_op_imag, state_prep_circuit, initial_statevector=initial_statevector) + return exp_real if (exp_imag == 0.) else exp_real + 1.0j * exp_imag + + def _get_expectation_value_from_statevector(self, qubit_operator, state_prep_circuit, initial_statevector=None): + r"""Take as input a qubit operator H and a state preparation returning a + ket |\psi>. Return the expectation value <\psi | H | \psi>, computed + without drawing samples (statevector only). Users should not be calling + this function directly, please call "get_expectation_value" instead. + + Args: + qubit_operator(openfermion-style QubitOperator class): a qubit + operator. + state_prep_circuit: an abstract circuit used for state preparation + (only pure states). + + Returns: + complex: The expectation value of this operator with regards to the + state preparation. + """ + n_qubits = state_prep_circuit.width + + expectation_value = 0. + prepared_frequencies, prepared_state = self.simulate(state_prep_circuit, return_statevector=True, initial_statevector=initial_statevector) + + if hasattr(self, "expectation_value_from_prepared_state"): + return self.expectation_value_from_prepared_state(qubit_operator, n_qubits, prepared_state) + + # Otherwise, use generic statevector expectation value + for term, coef in qubit_operator.terms.items(): + + if len(term) > n_qubits: # Cannot have a qubit index beyond circuit size + raise ValueError(f"Size of operator {qubit_operator} beyond circuit width ({n_qubits} qubits)") + elif not term: # Empty term: no simulation needed + expectation_value += coef + continue + + if not self.n_shots: + # Directly simulate and compute expectation value using statevector + pauli_circuit = Circuit([Gate(pauli, index) for index, pauli in term], n_qubits=n_qubits) + _, pauli_state = self.simulate(pauli_circuit, return_statevector=True, initial_statevector=prepared_state) + + delta = np.dot(pauli_state.real, prepared_state.real) + np.dot(pauli_state.imag, prepared_state.imag) + expectation_value += coef * delta + + else: + # Run simulation with statevector but compute expectation value with samples directly drawn from it + basis_circuit = Circuit(measurement_basis_gates(term), n_qubits=state_prep_circuit.width) + if basis_circuit.size > 0: + frequencies, _ = self.simulate(basis_circuit, initial_statevector=prepared_state) + else: + frequencies = prepared_frequencies + expectation_term = self.get_expectation_value_from_frequencies_oneterm(term, frequencies) + expectation_value += coef * expectation_term + + return expectation_value + + def _get_expectation_value_from_frequencies(self, qubit_operator, state_prep_circuit, initial_statevector=None): + r"""Take as input a qubit operator H and a state preparation returning a + ket |\psi>. Return the expectation value <\psi | H | \psi> computed + using the frequencies of observable states. + + Args: + qubit_operator(openfermion-style QubitOperator class): a qubit + operator. + state_prep_circuit: an abstract circuit used for state preparation. + + Returns: + complex: The expectation value of this operator with regards to the + state preparation. + """ + n_qubits = state_prep_circuit.width + if not self.statevector_available or state_prep_circuit.is_mixed_state or self._noise_model: + initial_circuit = state_prep_circuit + if initial_statevector is not None and not self.statevector_available: + raise ValueError(f'Backend {self.__class__} does not support statevectors') + else: + updated_statevector = initial_statevector + else: + initial_circuit = Circuit(n_qubits=n_qubits) + _, updated_statevector = self.simulate(state_prep_circuit, return_statevector=True, initial_statevector=initial_statevector) + + expectation_value = 0. + for term, coef in qubit_operator.terms.items(): + + if len(term) > n_qubits: + raise ValueError(f"Size of operator {qubit_operator} beyond circuit width ({n_qubits} qubits)") + elif not term: # Empty term: no simulation needed + expectation_value += coef + continue + + basis_circuit = Circuit(measurement_basis_gates(term)) + full_circuit = initial_circuit + basis_circuit if (basis_circuit.size > 0) else initial_circuit + frequencies, _ = self.simulate(full_circuit, initial_statevector=updated_statevector) + expectation_term = self.get_expectation_value_from_frequencies_oneterm(term, frequencies) + expectation_value += coef * expectation_term + + return expectation_value + + @staticmethod + def get_expectation_value_from_frequencies_oneterm(term, frequencies): + """Return the expectation value of a single-term qubit-operator, given + the result of a state-preparation. + + Args: + term(openfermion-style QubitOperator object): a qubit operator, with + only a single term. + frequencies(dict): histogram of frequencies of measurements (assumed + to be in lsq-first format). + + Returns: + complex: The expectation value of this operator with regards to the + state preparation. + """ + + return get_expectation_value_from_frequencies_oneterm(term, frequencies) + + def _statevector_to_frequencies(self, statevector): + """For a given statevector representing the quantum state of a qubit + register, returns a sparse histogram of the probabilities in the + least-significant-qubit (lsq) -first order. e.g the string '100' means + qubit 0 measured in basis state |1>, and qubit 1 & 2 both measured in + state |0>. + + Args: + statevector(list or ndarray(complex)): an iterable 1D data-structure + containing the amplitudes. + + Returns: + dict: A dictionary whose keys are bitstrings representing the + multi-qubit states with the least significant qubit first (e.g + '100' means qubit 0 in state |1>, and qubit 1 and 2 in state + |0>), and the associated value is the corresponding frequency. + Unless threshold=0., this dictionary will be sparse. + """ + + n_qubits = int(math.log2(len(statevector))) + frequencies = dict() + for i, amplitude in enumerate(statevector): + frequency = abs(amplitude)**2 + if (frequency - self.freq_threshold) >= 0.: + frequencies[self._int_to_binstr(i, n_qubits)] = frequency + + # If n_shots, has been specified, then draw that amount of samples from the distribution + # and return empirical frequencies instead. Otherwise, return the exact frequencies + if not self.n_shots: + return frequencies + else: + xk, pk = [], [] + for k, v in frequencies.items(): + xk.append(int(k[::-1], 2)) + pk.append(frequencies[k]) + distr = stats.rv_discrete(name='distr', values=(np.array(xk), np.array(pk))) + + # Generate samples from distribution. Cut in chunks to ensure samples fit in memory, gradually accumulate + chunk_size = 10**7 + n_chunks = self.n_shots // chunk_size + freqs_shots = Counter() + + for i in range(n_chunks+1): + this_chunk = self.n_shots % chunk_size if i == n_chunks else chunk_size + samples = distr.rvs(size=this_chunk) + freqs_shots += Counter(samples) + freqs_shots = {self._int_to_binstr(k, n_qubits, False): v / self.n_shots for k, v in freqs_shots.items()} + return freqs_shots + + def _int_to_binstr(self, i, n_qubits, use_ordering=True): + """Convert an integer into a bit string of size n_qubits. + + Args: + i (int): integer to convert to bit string. + n_qubits (int): The number of qubits and length of returned bit string. + use_ordering (bool): Flip the order of the returned bit string depending on self.statevector_order being "msq_first" or "lsq_first" + + Returns: + string: The bit string of the integer in lsq-first order. + """ + bs = bin(i).split('b')[-1] + state_binstr = "0" * (n_qubits - len(bs)) + bs + + return state_binstr if use_ordering and (self.statevector_order == "lsq_first") else state_binstr[::-1] + + @staticmethod + @abc.abstractmethod + def backend_info() -> dict: + """A dictionary that includes {'noisy_simulation': True or False, + 'statevector_available': True or False, + 'statevector_order': 'lsq_first' or 'msq_first'""" + pass diff --git a/tangelo/linq/target/__init__.py b/tangelo/linq/target/__init__.py new file mode 100644 index 000000000..277e32b5f --- /dev/null +++ b/tangelo/linq/target/__init__.py @@ -0,0 +1,25 @@ +# Copyright 2021 Good Chemistry Company. +# +# 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 .target_cirq import CirqSimulator +from .target_qiskit import QiskitSimulator +from .target_qulacs import QulacsSimulator +from .target_qdk import QDKSimulator +from tangelo.helpers.utils import all_backends_simulator + + +target_dict = {"qiskit": QiskitSimulator, "cirq": CirqSimulator, "qdk": QDKSimulator, "qulacs": QulacsSimulator} + +# Generate backend info dictionary +backend_info = {sim_id: target_dict[sim_id].backend_info() for sim_id in all_backends_simulator} diff --git a/tangelo/linq/target/target_cirq.py b/tangelo/linq/target/target_cirq.py new file mode 100644 index 000000000..adf38f2de --- /dev/null +++ b/tangelo/linq/target/target_cirq.py @@ -0,0 +1,110 @@ +# Copyright 2021 Good Chemistry Company. +# +# 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 collections import Counter + +import numpy as np + +from tangelo.linq import Circuit +from tangelo.linq.simulator_base import SimulatorBase +import tangelo.linq.translator as translator + + +class CirqSimulator(SimulatorBase): + + def __init__(self, n_shots=None, noise_model=None): + """Instantiate cirq simulator object. + + Args: + n_shots (int): Number of shots if using a shot-based simulator. + noise_model: A noise model object assumed to be in the format + expected from the target backend. + """ + import cirq + super().__init__(n_shots=n_shots, noise_model=noise_model) + self.cirq = cirq + + def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, initial_statevector=None): + """Perform state preparation corresponding to the input circuit on the + target backend, return the frequencies of the different observables, and + either the statevector or None depending on the availability of the + statevector and if return_statevector is set to True. For the + statevector backends supporting it, an initial statevector can be + provided to initialize the quantum state without simulating all the + equivalent gates. + + Args: + source_circuit: a circuit in the abstract format to be translated + for the target backend. + return_statevector(bool): option to return the statevector as well, + if available. + initial_statevector(list/array) : A valid statevector in the format + supported by the target backend. + + Returns: + dict: A dictionary mapping multi-qubit states to their corresponding + frequency. + numpy.array: The statevector, if available for the target backend + and requested by the user (if not, set to None). + """ + + translated_circuit = translator.translate_cirq(source_circuit, self._noise_model) + + if source_circuit.is_mixed_state or self._noise_model: + # Only DensityMatrixSimulator handles noise well, can use Simulator but it is slower + cirq_simulator = self.cirq.DensityMatrixSimulator(dtype=np.complex128) + else: + cirq_simulator = self.cirq.Simulator(dtype=np.complex128) + + # If requested, set initial state + cirq_initial_statevector = initial_statevector if initial_statevector is not None else 0 + + # Calculate final density matrix and sample from that for noisy simulation or simulating mixed states + if self._noise_model or source_circuit.is_mixed_state: + # cirq.dephase_measurements changes measurement gates to Krauss operators so simulators + # can be called once and density matrix sampled repeatedly. + translated_circuit = self.cirq.dephase_measurements(translated_circuit) + sim = cirq_simulator.simulate(translated_circuit, initial_state=cirq_initial_statevector) + self._current_state = sim.final_density_matrix + indices = list(range(source_circuit.width)) + isamples = self.cirq.sample_density_matrix(sim.final_density_matrix, indices, repetitions=self.n_shots) + samples = [''.join([str(int(q))for q in isamples[i]]) for i in range(self.n_shots)] + + frequencies = {k: v / self.n_shots for k, v in Counter(samples).items()} + # Noiseless simulation using the statevector simulator otherwise + else: + job_sim = cirq_simulator.simulate(translated_circuit, initial_state=cirq_initial_statevector) + self._current_state = job_sim.final_state_vector + frequencies = self._statevector_to_frequencies(self._current_state) + + return (frequencies, np.array(self._current_state)) if return_statevector else (frequencies, None) + + def expectation_value_from_prepared_state(self, qubit_operator, n_qubits, prepared_state): + + GATE_CIRQ = translator.get_cirq_gates() + qubit_labels = self.cirq.LineQubit.range(n_qubits) + qubit_map = {q: i for i, q in enumerate(qubit_labels)} + paulisum = 0.*self.cirq.PauliString(self.cirq.I(qubit_labels[0])) + for term, coef in qubit_operator.terms.items(): + pauli_list = [GATE_CIRQ[pauli](qubit_labels[index]) for index, pauli in term] + paulisum += self.cirq.PauliString(pauli_list, coefficient=coef) + if self._noise_model: + exp_value = paulisum.expectation_from_density_matrix(prepared_state, qubit_map) + else: + exp_value = paulisum.expectation_from_state_vector(prepared_state, qubit_map) + return np.real(exp_value) + + @staticmethod + def backend_info(): + return {"statevector_available": True, "statevector_order": "lsq_first", "noisy_simulation": True} diff --git a/tangelo/linq/target/target_qdk.py b/tangelo/linq/target/target_qdk.py new file mode 100644 index 000000000..0ac3c53e1 --- /dev/null +++ b/tangelo/linq/target/target_qdk.py @@ -0,0 +1,69 @@ +# Copyright 2021 Good Chemistry Company. +# +# 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 tangelo.linq import Circuit +from tangelo.linq.simulator_base import SimulatorBase +import tangelo.linq.translator as translator + + +class QDKSimulator(SimulatorBase): + + def __init__(self, n_shots=None, noise_model=None): + import qsharp + super().__init__(n_shots=n_shots, noise_model=noise_model) + self.qsharp = qsharp + + def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, initial_statevector=None): + """Perform state preparation corresponding to the input circuit on the + target backend, return the frequencies of the different observables, and + either the statevector or None depending on the availability of the + statevector and if return_statevector is set to True. For the + statevector backends supporting it, an initial statevector can be + provided to initialize the quantum state without simulating all the + equivalent gates. + + Args: + source_circuit: a circuit in the abstract format to be translated + for the target backend. + return_statevector(bool): option to return the statevector as well, + if available. + initial_statevector(list/array) : A valid statevector in the format + supported by the target backend. + + Returns: + dict: A dictionary mapping multi-qubit states to their corresponding + frequency. + numpy.array: The statevector, if available for the target backend + and requested by the user (if not, set to None). + """ + translated_circuit = translator.translate_qsharp(source_circuit) + with open('tmp_circuit.qs', 'w+') as f_out: + f_out.write(translated_circuit) + + # Compile, import and call Q# operation to compute frequencies. Only import qsharp module if qdk is running + # TODO: A try block to catch an exception at compile time, for Q#? Probably as an ImportError. + self.qsharp.reload() + from MyNamespace import EstimateFrequencies + frequencies_list = EstimateFrequencies.simulate(nQubits=source_circuit.width, nShots=self.n_shots) + print("Q# frequency estimation with {0} shots: \n {1}".format(self.n_shots, frequencies_list)) + + # Convert Q# output to frequency dictionary, apply threshold + frequencies = {bin(i).split('b')[-1]: freq for i, freq in enumerate(frequencies_list)} + frequencies = {("0"*(source_circuit.width-len(k))+k)[::-1]: v for k, v in frequencies.items() + if v > self.freq_threshold} + return (frequencies, None) + + @staticmethod + def backend_info(): + return {"statevector_available": False, "statevector_order": None, "noisy_simulation": False} diff --git a/tangelo/linq/target/target_qiskit.py b/tangelo/linq/target/target_qiskit.py new file mode 100644 index 000000000..c0ce04d93 --- /dev/null +++ b/tangelo/linq/target/target_qiskit.py @@ -0,0 +1,100 @@ +# Copyright 2021 Good Chemistry Company. +# +# 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 math + +import numpy as np + +from tangelo.linq import Circuit +from tangelo.linq.simulator_base import SimulatorBase +import tangelo.linq.translator as translator + + +class QiskitSimulator(SimulatorBase): + """Interface to the qiskit simulator.""" + + def __init__(self, n_shots=None, noise_model=None): + import qiskit + from qiskit_aer import AerSimulator + super().__init__(n_shots=n_shots, noise_model=noise_model) + self.qiskit = qiskit + self.AerSimulator = AerSimulator + + def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, initial_statevector=None): + """Perform state preparation corresponding to the input circuit on the + target backend, return the frequencies of the different observables, and + either the statevector or None depending on the availability of the + statevector and if return_statevector is set to True. For the + statevector backends supporting it, an initial statevector can be + provided to initialize the quantum state without simulating all the + equivalent gates. + + Args: + source_circuit: a circuit in the abstract format to be translated + for the target backend. + return_statevector(bool): option to return the statevector as well, + if available. + initial_statevector(list/array) : A valid statevector in the format + supported by the target backend. + + Returns: + dict: A dictionary mapping multi-qubit states to their corresponding + frequency. + numpy.array: The statevector, if available for the target backend + and requested by the user (if not, set to None). + """ + + translated_circuit = translator.translate_qiskit(source_circuit) + + # If requested, set initial state + if initial_statevector is not None: + if self._noise_model: + raise ValueError("Cannot load an initial state if using a noise model, with Qiskit") + else: + n_qubits = int(math.log2(len(initial_statevector))) + initial_state_circuit = self.qiskit.QuantumCircuit(n_qubits, n_qubits) + initial_state_circuit.initialize(initial_statevector, list(range(n_qubits))) + translated_circuit = initial_state_circuit.compose(translated_circuit) + + # Drawing individual shots with the qasm simulator, for noisy simulation or simulating mixed states + if self._noise_model or source_circuit.is_mixed_state: + from tangelo.linq.noisy_simulation.noise_models import get_qiskit_noise_model + + meas_range = range(source_circuit.width) + translated_circuit.measure(meas_range, meas_range) + return_statevector = False + backend = self.AerSimulator() + + qiskit_noise_model = get_qiskit_noise_model(self._noise_model) if self._noise_model else None + opt_level = 0 if self._noise_model else None + + job_sim = self.qiskit.execute(translated_circuit, backend, noise_model=qiskit_noise_model, + shots=self.n_shots, basis_gates=None, optimization_level=opt_level) + sim_results = job_sim.result() + frequencies = {state[::-1]: count/self.n_shots for state, count in sim_results.get_counts(0).items()} + + # Noiseless simulation using the statevector simulator otherwise + else: + backend = self.AerSimulator(method='statevector') + translated_circuit = self.qiskit.transpile(translated_circuit, backend) + translated_circuit.save_statevector() + sim_results = backend.run(translated_circuit).result() + self._current_state = np.asarray(sim_results.get_statevector(translated_circuit)) + frequencies = self._statevector_to_frequencies(self._current_state) + + return (frequencies, np.array(sim_results.get_statevector())) if return_statevector else (frequencies, None) + + @staticmethod + def backend_info(): + return {"statevector_available": True, "statevector_order": "msq_first", "noisy_simulation": True} diff --git a/tangelo/linq/target/target_qulacs.py b/tangelo/linq/target/target_qulacs.py new file mode 100644 index 000000000..2ee71430d --- /dev/null +++ b/tangelo/linq/target/target_qulacs.py @@ -0,0 +1,106 @@ +# Copyright 2021 Good Chemistry Company. +# +# 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 os +from collections import Counter + +import numpy as np + +from tangelo.linq import Circuit +from tangelo.linq.simulator_base import SimulatorBase +import tangelo.linq.translator as translator + + +class QulacsSimulator(SimulatorBase): + + def __init__(self, n_shots=None, noise_model=None): + import qulacs + super().__init__(n_shots=n_shots, noise_model=noise_model) + self.qulacs = qulacs + + def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, initial_statevector=None): + """Perform state preparation corresponding to the input circuit on the + target backend, return the frequencies of the different observables, and + either the statevector or None depending on the availability of the + statevector and if return_statevector is set to True. For the + statevector backends supporting it, an initial statevector can be + provided to initialize the quantum state without simulating all the + equivalent gates. + + Args: + source_circuit: a circuit in the abstract format to be translated + for the target backend. + return_statevector(bool): option to return the statevector as well, + if available. + initial_statevector(list/array) : A valid statevector in the format + supported by the target backend. + + Returns: + dict: A dictionary mapping multi-qubit states to their corresponding + frequency. + numpy.array: The statevector, if available for the target backend + and requested by the user (if not, set to None). + """ + + translated_circuit = translator.translate_qulacs(source_circuit, self._noise_model) + + # Initialize state on GPU if available and desired. Default to CPU otherwise. + if ('QuantumStateGpu' in dir(self.qulacs)) and (int(os.getenv("QULACS_USE_GPU", 0)) != 0): + state = self.qulacs.QuantumStateGpu(source_circuit.width) + else: + state = self.qulacs.QuantumState(source_circuit.width) + if initial_statevector is not None: + state.load(initial_statevector) + + if (source_circuit.is_mixed_state or self._noise_model): + samples = list() + for i in range(self.n_shots): + translated_circuit.update_quantum_state(state) + samples.append(state.sampling(1)[0]) + if initial_statevector is not None: + state.load(initial_statevector) + else: + state.set_zero_state() + python_statevector = None + elif self.n_shots is not None: + translated_circuit.update_quantum_state(state) + self._current_state = state + python_statevector = np.array(state.get_vector()) if return_statevector else None + samples = state.sampling(self.n_shots) + else: + translated_circuit.update_quantum_state(state) + self._current_state = state + python_statevector = state.get_vector() + frequencies = self._statevector_to_frequencies(python_statevector) + return (frequencies, np.array(python_statevector)) if return_statevector else (frequencies, None) + + frequencies = {self._int_to_binstr(k, source_circuit.width): v / self.n_shots + for k, v in Counter(samples).items()} + return (frequencies, python_statevector) + + def expectation_value_from_prepared_state(self, qubit_operator, n_qubits, prepared_state): + + # TODO: This section previously used qulacs.quantum_operator.create_quantum_operator_from_openfermion_text but was changed + # due to a memory leak. We can re-evaluate the implementation if/when Issue #303 (https://github.com/qulacs/qulacs/issues/303) + # is fixed. + operator = self.qulacs.Observable(n_qubits) + for term, coef in qubit_operator.terms.items(): + pauli_string = "".join(f" {op} {qu}" for qu, op in term) + operator.add_operator(coef, pauli_string) + return operator.get_expectation_value(self._current_state).real + + @staticmethod + def backend_info(): + return {"statevector_available": True, "statevector_order": "msq_first", "noisy_simulation": True} diff --git a/tangelo/linq/tests/test_simulator.py b/tangelo/linq/tests/test_simulator.py index efa50f170..e1739e96b 100644 --- a/tangelo/linq/tests/test_simulator.py +++ b/tangelo/linq/tests/test_simulator.py @@ -27,6 +27,7 @@ from tangelo.linq import Gate, Circuit, translator, Simulator from tangelo.linq.gate import PARAMETERIZED_GATES from tangelo.helpers.utils import installed_simulator, installed_sv_simulator, installed_backends +from tangelo.linq.simulator_base import SimulatorBase, get_expectation_value_from_frequencies_oneterm path_data = os.path.dirname(os.path.abspath(__file__)) + '/data' @@ -179,7 +180,6 @@ def test_simulate_empty_circuit_from_statevector(self): def test_get_exp_value_from_statevector(self): """ Compute the expectation value from the statevector for each statevector backend """ - for b in installed_sv_simulator: simulator = Simulator(target=b) exp_values = np.zeros((len(circuits), len(ops)), dtype=float) @@ -361,6 +361,7 @@ def test_get_exp_value_from_frequencies(self): class TestSimulateMisc(unittest.TestCase): + @unittest.skipIf("qdk" not in installed_backends, "Test Skipped: Backend not available \n") def test_n_shots_needed(self): """ Raise an error if user chooses a target backend that does not provide access to a statevector and @@ -396,9 +397,51 @@ def test_get_exp_value_from_frequencies_oneterm(self): are being provided as input. """ term, coef = ((0, 'Z'),), 1.0 # Data as presented in Openfermion's QubitOperator.terms attribute - exp_value = coef * Simulator.get_expectation_value_from_frequencies_oneterm(term, ref_freqs[2]) + exp_value = coef * get_expectation_value_from_frequencies_oneterm(term, ref_freqs[2]) np.testing.assert_almost_equal(exp_value, -0.41614684, decimal=5) + def test_invalid_target(self): + """ Ensure an error is returned if the target simulator is not supported.""" + self.assertRaises(ValueError, Simulator, 'banana') + + def test_user_provided_simulator(self): + """Test user defined target simulator that disregards the circuit gates and only returns zero state or one state""" + + class TrueFalseSimulator(SimulatorBase): + def __init__(self, n_shots=None, noise_model=None, return_zeros=True): + """Instantiate simulator object that always returns all zeros or all ones ignoring circuit operations.""" + super().__init__(n_shots=n_shots, noise_model=noise_model) + self.return_zeros = return_zeros + + def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, initial_statevector=None): + """Perform state preparation corresponding self.return_zeros.""" + + statevector = np.zeros(2**source_circuit.width, dtype=complex) + if self.return_zeros: + statevector[0] = 1. + else: + statevector[-1] = 1. + + frequencies = self._statevector_to_frequencies(statevector) + + return (frequencies, np.array(statevector)) if return_statevector else (frequencies, None) + + @staticmethod + def backend_info(): + return {"statevector_available": True, "statevector_order": "msq_first", "noisy_simulation": False} + + sim = Simulator(TrueFalseSimulator, n_shots=1, noise_model=None, return_zeros=True) + f, sv = sim.simulate(circuit1, return_statevector=True) + assert_freq_dict_almost_equal(f, {"00": 1}, 1.e-7) + np.testing.assert_almost_equal(np.array([1., 0., 0., 0.]), sv) + self.assertAlmostEqual(sim.get_expectation_value(QubitOperator("Z0", 1.), circuit1), 1.) + + sim = Simulator(TrueFalseSimulator, n_shots=1, noise_model=None, return_zeros=False) + f, sv = sim.simulate(circuit1, return_statevector=True) + assert_freq_dict_almost_equal(f, {"11": 1}, 1.e-7) + np.testing.assert_almost_equal(np.array([0., 0., 0., 1.]), sv) + self.assertAlmostEqual(sim.get_expectation_value(QubitOperator("Z0", 1.), circuit1), -1.) + if __name__ == "__main__": unittest.main() diff --git a/tangelo/linq/translator/__init__.py b/tangelo/linq/translator/__init__.py index 7ee0dd7c7..f3a685362 100644 --- a/tangelo/linq/translator/__init__.py +++ b/tangelo/linq/translator/__init__.py @@ -25,13 +25,15 @@ from .translate_qubitop import translate_operator -# List all supported gates for all backends found -SUPPORTED_GATES = dict() +def get_supported_gates(): + "List all supported gates for all backends found" + SUPPORTED_GATES = dict() -SUPPORTED_GATES["projectq"] = sorted(get_projectq_gates().keys()) -SUPPORTED_GATES["ionq"] = sorted(get_ionq_gates().keys()) -SUPPORTED_GATES["qdk"] = sorted(get_qdk_gates().keys()) -SUPPORTED_GATES["openqasm"] = sorted(get_openqasm_gates().keys()) + SUPPORTED_GATES["projectq"] = sorted(get_projectq_gates().keys()) + SUPPORTED_GATES["ionq"] = sorted(get_ionq_gates().keys()) + SUPPORTED_GATES["qdk"] = sorted(get_qdk_gates().keys()) + SUPPORTED_GATES["openqasm"] = sorted(get_openqasm_gates().keys()) -for v in installed_backends: - SUPPORTED_GATES[v] = sorted(eval(f'get_{v}_gates().keys()')) + for v in installed_backends: + SUPPORTED_GATES[v] = sorted(eval(f'get_{v}_gates().keys()')) + return SUPPORTED_GATES diff --git a/tangelo/toolboxes/circuits/tests/test_lcu.py b/tangelo/toolboxes/circuits/tests/test_lcu.py index 80c86b634..73d5e9135 100644 --- a/tangelo/toolboxes/circuits/tests/test_lcu.py +++ b/tangelo/toolboxes/circuits/tests/test_lcu.py @@ -19,7 +19,7 @@ import numpy as np from scipy.linalg import expm -from tangelo.linq import Simulator, backend_info +from tangelo.linq import Simulator from tangelo.helpers.utils import installed_backends from tangelo.linq.helpers.circuits.statevector import StateVector from tangelo.toolboxes.operators.operators import QubitOperator @@ -53,7 +53,7 @@ def test_get_truncated_taylor_series(self): for backend in backends: sim = Simulator(backend) - statevector_order = backend_info[backend]["statevector_order"] + statevector_order = sim.backend_info()["statevector_order"] sv = StateVector(vec, order=statevector_order) sv_circuit = sv.initializing_circuit() @@ -88,7 +88,7 @@ def test_get_oaa_lcu_circuit(self): for backend in backends: sim = Simulator(backend) - statevector_order = backend_info[backend]["statevector_order"] + statevector_order = sim.backend_info()["statevector_order"] sv = StateVector(vec, order=statevector_order) sv_circuit = sv.initializing_circuit() diff --git a/tangelo/toolboxes/measurements/qubit_terms_grouping.py b/tangelo/toolboxes/measurements/qubit_terms_grouping.py index 14be83bb8..e532e3f56 100644 --- a/tangelo/toolboxes/measurements/qubit_terms_grouping.py +++ b/tangelo/toolboxes/measurements/qubit_terms_grouping.py @@ -21,7 +21,7 @@ import warnings from openfermion.measurements import group_into_tensor_product_basis_sets -from tangelo.linq import Simulator +from tangelo.linq import get_expectation_value_from_frequencies_oneterm def group_qwc(qb_ham, seed=None, n_repeat=1): @@ -133,6 +133,6 @@ def exp_value_from_measurement_bases(sub_ops, histograms): exp_value = 0. for basis, freqs in histograms.items(): for term, coef in sub_ops[basis].terms.items(): - exp_value += Simulator.get_expectation_value_from_frequencies_oneterm(term, freqs) * coef + exp_value += get_expectation_value_from_frequencies_oneterm(term, freqs) * coef return exp_value diff --git a/tangelo/toolboxes/post_processing/histogram.py b/tangelo/toolboxes/post_processing/histogram.py index 8c4c32a5b..9ebd56bce 100644 --- a/tangelo/toolboxes/post_processing/histogram.py +++ b/tangelo/toolboxes/post_processing/histogram.py @@ -18,7 +18,7 @@ from collections import Counter -from tangelo.linq import Simulator +from tangelo.linq import get_expectation_value_from_frequencies_oneterm from tangelo.toolboxes.post_processing.bootstrapping import get_resampled_frequencies @@ -179,7 +179,7 @@ def get_expectation_value(self, term, coeff=1.): Returns: imaginary: Expectation value for this operator. """ - return coeff*Simulator.get_expectation_value_from_frequencies_oneterm(term, self.frequencies) + return coeff*get_expectation_value_from_frequencies_oneterm(term, self.frequencies) def aggregate_histograms(*hists):