Enable arbitrary state initialization for Qulacs and Qiskit backends (#…

…371)

src/tequila/simulators/simulator_api.py

@@ -9,6 +9,7 @@
 from tequila.utils.exceptions import TequilaException, TequilaWarning
 from tequila.simulators.simulator_base import BackendCircuit, BackendExpectationValue
 from tequila.circuit.noise import NoiseModel
+from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
 
 SUPPORTED_BACKENDS = ["qulacs", "qulacs_gpu", "qibo", "qiskit", "qiskit_gpu", "cirq", "pyquil", "symbolic", "qlm"]
 SUPPORTED_NOISE_BACKENDS = ["qiskit", "qiskit_gpu", "cirq", "pyquil"]  # qulacs removed in v.1.9
@@ -22,7 +23,6 @@
     from tequila.objective import Objective, Variable
     from tequila.circuit.gates import QCircuit
     import numbers.Real as RealNumber
-    from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
 
 """
 Check which simulators are installed
@@ -369,8 +369,9 @@ def simulate(objective: typing.Union['Objective', 'QCircuit', 'QTensor'],
              backend: str = None,
              noise: NoiseModel = None,
              device: str = None,
+             initial_state: Union[int, QubitWaveFunction] = 0,
              *args,
-             **kwargs) -> Union[RealNumber, 'QubitWaveFunction']:
+             **kwargs) -> Union[RealNumber, QubitWaveFunction]:
     """Simulate a tequila objective or circuit
 
     Parameters
@@ -388,6 +389,8 @@ def simulate(objective: typing.Union['Objective', 'QCircuit', 'QTensor'],
         specify a noise model to apply to simulation/sampling
     device:
         a device upon which (or in emulation of which) to sample
+    initial_state: int or QubitWaveFunction:
+        the initial state of the circuit
     *args :
 
     **kwargs :
@@ -409,7 +412,7 @@ def simulate(objective: typing.Union['Objective', 'QCircuit', 'QTensor'],
     compiled_objective = compile(objective=objective, samples=samples, variables=variables, backend=backend,
                                  noise=noise, device=device, *args, **kwargs)
 
-    return compiled_objective(variables=variables, samples=samples, *args, **kwargs)
+    return compiled_objective(variables=variables, samples=samples, initial_state=initial_state, *args, **kwargs)
 
 
 def draw(objective, variables=None, backend: str = None, name=None, *args, **kwargs):

src/tequila/simulators/simulator_base.py

@@ -7,6 +7,7 @@
 from tequila import BitString
 from tequila.objective.objective import Variable, format_variable_dictionary
 from tequila.circuit import compiler
+from typing import Union
 
 import numbers, typing, numpy, copy, warnings
 
@@ -107,6 +108,11 @@ class BackendCircuit():
         "cc_max": True
     }
 
+    # Can be overwritten by backends that allow basis state initialization when sampling
+    supports_sampling_initialization: bool = False
+    # Can be overwritten by backends that allow initializing arbitrary states
+    supports_generic_initialization: bool = False
+
     @property
     def n_qubits(self) -> numbers.Integral:
         return len(self.qubit_map)
@@ -328,7 +334,7 @@ def update_variables(self, variables):
         """
         self.circuit = self.create_circuit(abstract_circuit=self.abstract_circuit, variables=variables)
 
-    def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunction:
+    def simulate(self, variables, initial_state: Union[int, QubitWaveFunction] = 0, *args, **kwargs) -> QubitWaveFunction:
         """
         simulate the circuit via the backend.
 
@@ -348,25 +354,34 @@ def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunc
             the wavefunction of the system produced by the action of the circuit on the initial state.
 
         """
+        if isinstance(initial_state, QubitWaveFunction) and not self.supports_generic_initialization:
+            raise TequilaException("Backend does not support arbitrary initial states")
+
         self.update_variables(variables)
         if isinstance(initial_state, BitString):
             initial_state = initial_state.integer
-        if isinstance(initial_state, QubitWaveFunction):
-            if initial_state.length() != 1:
-                raise TequilaException("only product states as initial states accepted")
-            initial_state = list(initial_state.keys())[0].integer
 
         all_qubits = list(range(self.abstract_circuit.n_qubits))
         active_qubits = self.qubit_map.keys()
 
         # Keymap is only necessary if not all qubits are active
         keymap_required = sorted(active_qubits) != all_qubits
 
+        # Combining keymap and general initial states is awkward, because it's not clear what should happen with
+        # different states on non-active qubits. For now, this is simply not allowed.
+        # A better solution might be to check if all components of the initial state differ only on the active qubits.
+        if keymap_required and isinstance(initial_state, QubitWaveFunction):
+            raise TequilaException("Can only set non-basis initial state if all qubits are used")
+
         if keymap_required:
             # maps from reduced register to full register
             keymap = KeyMapSubregisterToRegister(subregister=active_qubits, register=all_qubits)
 
-        mapped_initial_state = keymap.inverted(initial_state).integer if keymap_required else int(initial_state)
+        if not isinstance(initial_state, QubitWaveFunction):
+            mapped_initial_state = keymap.inverted(initial_state).integer if keymap_required else int(initial_state)
+        else:
+            mapped_initial_state = initial_state
+
         result = self.do_simulate(variables=variables, initial_state=mapped_initial_state, *args,
                                   **kwargs)
 
@@ -375,7 +390,7 @@ def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunc
 
         return result
 
-    def sample(self, variables, samples, read_out_qubits=None, circuit=None, *args, **kwargs):
+    def sample(self, variables, samples, read_out_qubits=None, circuit=None, initial_state=0, *args, **kwargs):
         """
         Sample the circuit. If circuit natively equips paulistrings, sample therefrom.
         Parameters
@@ -395,6 +410,12 @@ def sample(self, variables, samples, read_out_qubits=None, circuit=None, *args,
             The result of sampling, a recreated QubitWaveFunction in the sampled basis.
 
         """
+        if initial_state != 0 and not self.supports_sampling_initialization:
+            raise TequilaException("Backend does not support initial states for sampling")
+
+        if isinstance(initial_state, QubitWaveFunction) and not self.supports_generic_initialization:
+            raise TequilaException("Backend does not support arbitrary initial states")
+
         self.update_variables(variables)
         if read_out_qubits is None:
             read_out_qubits = self.abstract_qubits
@@ -406,7 +427,9 @@ def sample(self, variables, samples, read_out_qubits=None, circuit=None, *args,
             circuit = self.add_measurement(circuit=self.circuit, target_qubits=read_out_qubits)
         else:
             circuit = self.add_measurement(circuit=circuit, target_qubits=read_out_qubits)
-        return self.do_sample(samples=samples, circuit=circuit, read_out_qubits=read_out_qubits, *args, **kwargs)
+
+        return self.do_sample(samples=samples, circuit=circuit, read_out_qubits=read_out_qubits,
+                              initial_state=initial_state, *args, **kwargs)
 
     def sample_all_z_hamiltonian(self, samples: int, hamiltonian, variables, *args, **kwargs):
         """
@@ -511,7 +534,7 @@ def sample_paulistring(self, samples: int, paulistring, variables, *args,
         E = E / samples * paulistring.coeff
         return E
 
-    def do_sample(self, samples, circuit, noise, abstract_qubits=None, *args, **kwargs) -> QubitWaveFunction:
+    def do_sample(self, samples, circuit, noise, abstract_qubits=None, initial_state=0, *args, **kwargs) -> QubitWaveFunction:
         """
         helper function for sampling. MUST be overwritten by inheritors.
 
@@ -777,7 +800,7 @@ def __call__(self, variables, samples: int = None, *args, **kwargs):
                 raise TequilaException(
                     "BackendExpectationValue received not all variables. Circuit depends on variables {}, you gave {}".format(
                         self._variables, variables))
-        
+
         if samples is None:
             data = self.simulate(variables=variables, *args, **kwargs)
         else:
@@ -856,7 +879,7 @@ def sample(self, variables, samples, *args, **kwargs) -> numpy.array:
             samples = max(1, int(self.abstract_expectationvalue.samples * total_samples))
             suggested = samples
             # samples are not necessarily set (either the user has to set it or some functions like optimize_measurements)
- 
+
         if suggested is not None and suggested != samples:
             warnings.warn("simulating with samples={}, but expectationvalue carries suggested samples={}\nTry calling with samples='auto-total#ofsamples'".format(samples, suggested), TequilaWarning)
 

src/tequila/simulators/simulator_qiskit.py

@@ -5,9 +5,11 @@
 from tequila import BitString, BitNumbering, BitStringLSB
 from tequila.utils.keymap import KeyMapRegisterToSubregister
 from tequila.utils import to_float
+from typing import Union
 import warnings
 import numpy as np
 import qiskit, qiskit_aer, qiskit.providers.fake_provider
+from qiskit import QuantumCircuit
 
 HAS_NOISE = True
 try:
@@ -160,6 +162,9 @@ class BackendCircuitQiskit(BackendCircuit):
 
     numbering = BitNumbering.LSB
 
+    supports_sampling_initialization = True
+    supports_generic_initialization = True
+
     def __init__(self, abstract_circuit: QCircuit, variables, qubit_map=None, noise=None,
                  device=None, *args, **kwargs):
         """
@@ -262,6 +267,20 @@ def make_classical_map(self, qubit_map: dict):
 
         return classical_map
 
+    def add_state_init(self, circuit: QuantumCircuit, initial_state: Union[int, QubitWaveFunction]) -> QuantumCircuit:
+        if initial_state == 0:
+            return circuit
+
+        if isinstance(initial_state, QubitWaveFunction):
+            statevector = initial_state.to_array(self.numbering)
+        else:
+            statevector = np.zeros(2 ** self.n_qubits)
+            statevector[reverse_int_bits(initial_state, self.n_qubits)] = 1.0
+
+        init_circuit = qiskit.QuantumCircuit(self.q, self.c)
+        init_circuit.set_statevector(statevector)
+        return init_circuit.compose(circuit)
+
     def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunction:
         """
         Helper function for performing simulation.
@@ -298,20 +317,15 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
 
         circuit = self.circuit.assign_parameters(self.resolver)
 
-        if initial_state != 0:
-            state = np.zeros(2 ** self.n_qubits)
-            state[reverse_int_bits(initial_state, self.n_qubits)] = 1.0
-            init_circuit = qiskit.QuantumCircuit(self.q, self.c)
-            init_circuit.set_statevector(state)
-            circuit = init_circuit.compose(circuit)
+        circuit = self.add_state_init(circuit, initial_state)
 
         circuit.save_statevector()
 
         backend_result = qiskit_backend.run(circuit, optimization_level=optimization_level).result()
 
         return QubitWaveFunction.from_array(array=backend_result.get_statevector(circuit).data, numbering=self.numbering)
 
-    def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubits, *args,
+    def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubits, initial_state=0, *args,
                   **kwargs) -> QubitWaveFunction:
         """
         Helper function for performing sampling.
@@ -321,6 +335,8 @@ def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubit
             the circuit from which to sample.
         samples:
             the number of samples to take.
+        initial_state:
+            initial state of the circuit
         args
         kwargs
 
@@ -371,14 +387,14 @@ def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubit
                 else:
                     use_basis = qiskit_backend.configuration().basis_gates
                 circuit = circuit.assign_parameters(self.resolver)  # this is necessary -- see qiskit-aer issue 1346
+                circuit = self.add_state_init(circuit, initial_state)
                 circuit = qiskit.transpile(circuit, backend=qiskit_backend,
                                            basis_gates=use_basis,
                                            optimization_level=optimization_level
                                            )
 
                 job = qiskit_backend.run(circuit, shots=samples)
-                return self.convert_measurements(job,
-                                                 target_qubits=read_out_qubits)
+                return self.convert_measurements(job, target_qubits=read_out_qubits)
 
     def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveFunction:
         """

src/tequila/simulators/simulator_qulacs.py

@@ -1,3 +1,5 @@
+from typing import Union
+
 import qulacs
 import numbers, numpy
 import warnings
@@ -63,6 +65,11 @@ class BackendCircuitQulacs(BackendCircuit):
 
     numbering = BitNumbering.LSB
 
+    quantum_state_class = qulacs.QuantumState
+
+    supports_sampling_initialization = True
+    supports_generic_initialization = True
+
     def __init__(self, abstract_circuit, noise=None, *args, **kwargs):
         """
 
@@ -110,10 +117,17 @@ def __init__(self, abstract_circuit, noise=None, *args, **kwargs):
 
             self.circuit=self.add_noise_to_circuit(noise)
 
-    def initialize_state(self, n_qubits:int=None) -> qulacs.QuantumState:
+    def initialize_state(self, n_qubits: int = None, initial_state: Union[int, QubitWaveFunction] = None) -> qulacs.QuantumState:
         if n_qubits is None:
             n_qubits = self.n_qubits
-        return qulacs.QuantumState(n_qubits)
+
+        state = self.quantum_state_class(n_qubits)
+        if isinstance(initial_state, int):
+            state.set_computational_basis(reverse_int_bits(initial_state, self.n_qubits))
+        elif isinstance(initial_state, QubitWaveFunction):
+            state.load(initial_state.to_array(self.numbering))
+
+        return state
 
     def update_variables(self, variables):
         """
@@ -130,7 +144,7 @@ def update_variables(self, variables):
         for k, angle in enumerate(self.variables):
             self.circuit.set_parameter(k, angle(variables))
 
-    def do_simulate(self, variables, initial_state, *args, **kwargs):
+    def do_simulate(self, variables, initial_state: Union[int, QubitWaveFunction], *args, **kwargs):
         """
         Helper function to perform simulation.
 
@@ -148,8 +162,7 @@ def do_simulate(self, variables, initial_state, *args, **kwargs):
         QubitWaveFunction:
             QubitWaveFunction representing result of the simulation.
         """
-        state = self.initialize_state(self.n_qubits)
-        state.set_computational_basis(reverse_int_bits(initial_state, self.n_qubits))
+        state = self.initialize_state(self.n_qubits, initial_state)
         self.circuit.update_quantum_state(state)
 
         wfn = QubitWaveFunction.from_array(array=state.get_vector(), numbering=self.numbering)
@@ -205,8 +218,7 @@ def do_sample(self, samples, circuit, noise_model=None, initial_state=0, *args,
         QubitWaveFunction:
             the results of sampling, as a Qubit Wave Function.
         """
-        state = self.initialize_state(self.n_qubits)
-        state.set_computational_basis(reverse_int_bits(initial_state, self.n_qubits))
+        state = self.initialize_state(self.n_qubits, initial_state)
         circuit.update_quantum_state(state)
         sampled = state.sampling(samples)
         return self.convert_measurements(backend_result=sampled, target_qubits=self.measurements)

src/tequila/simulators/simulator_qulacs_gpu.py

@@ -1,16 +1,18 @@
 import qulacs
+from qulacs_core import QuantumStateGpu
+
 from tequila import TequilaException
 from tequila.simulators.simulator_qulacs import BackendCircuitQulacs, BackendExpectationValueQulacs
 
+
 class TequilaQulacsGpuException(TequilaException):
     def __str__(self):
         return "Error in qulacs gpu backend:" + self.message
 
+
 class BackendCircuitQulacsGpu(BackendCircuitQulacs):
-    def initialize_state(self, n_qubits:int=None) -> qulacs.QuantumState:
-        if n_qubits is None:
-            n_qubits = self.n_qubits
-        return qulacs.QuantumStateGpu(n_qubits)
+    quantum_state_class = QuantumStateGpu
+
 
 class BackendExpectationValueQulacsGpu(BackendExpectationValueQulacs):
     BackendCircuitType = BackendCircuitQulacsGpu