Add initial state support for expectation values

src/tequila/objective/objective.py

@@ -1,6 +1,7 @@
 import typing, copy, numbers
+from typing import Union
 from tequila.grouping.compile_groups import compile_commuting_parts
-from tequila import TequilaException
+from tequila import TequilaException, QubitWaveFunction
 from tequila.utils import JoinedTransformation
 from tequila.hamiltonian import paulis
 import numpy as onp
@@ -545,14 +546,16 @@ def __str__(self):
                "variables          = {}\n" \
                "types              = {}".format(unique, measurements, variables, types)
 
-    def __call__(self, variables=None, *args, **kwargs):
+    def __call__(self, variables=None, initial_state: Union[int, QubitWaveFunction] = 0, *args, **kwargs):
         """
         Return the output of the calculation the objective represents.
 
         Parameters
         ----------
         variables: dict:
             dictionary instantiating all variables that may appear within the objective.
+        initial_state: int or QubitWaveFunction:
+            the initial state of the circuit
         args
         kwargs
 
@@ -579,7 +582,7 @@ def __call__(self, variables=None, *args, **kwargs):
         ev_array = []
         for E in self.args:
             if E not in evaluated:  #
-                expval_result = E(variables=variables, *args, **kwargs)
+                expval_result = E(variables=variables, initial_state=initial_state, *args, **kwargs)
                 evaluated[E] = expval_result
             else:
                 expval_result = evaluated[E]

src/tequila/simulators/simulator_base.py

@@ -431,7 +431,8 @@ def sample(self, variables, samples, read_out_qubits=None, circuit=None, initial
         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):
+    def sample_all_z_hamiltonian(self, samples: int, hamiltonian, variables,
+                                 initial_state: Union[int, QubitWaveFunction] = 0, *args, **kwargs):
         """
         Sample from a Hamiltonian which only consists of Pauli-Z and unit operators
         Parameters
@@ -440,6 +441,8 @@ def sample_all_z_hamiltonian(self, samples: int, hamiltonian, variables, *args,
             number of samples to take
         hamiltonian
             the tequila hamiltonian
+        initial_state
+            the initial state of the circuit
         args
             arguments for do_sample
         kwargs
@@ -458,7 +461,7 @@ def sample_all_z_hamiltonian(self, samples: int, hamiltonian, variables, *args,
                                                                                                  self.n_qubits))
 
         # run simulators
-        counts = self.sample(samples=samples, read_out_qubits=abstract_qubits_H, variables=variables, *args, **kwargs)
+        counts = self.sample(samples=samples, read_out_qubits=abstract_qubits_H, variables=variables, initial_state=initial_state, *args, **kwargs)
         read_out_map = {q: i for i, q in enumerate(abstract_qubits_H)}
 
         # compute energy
@@ -481,8 +484,8 @@ def sample_all_z_hamiltonian(self, samples: int, hamiltonian, variables, *args,
             assert n_samples == samples
         return E
 
-    def sample_paulistring(self, samples: int, paulistring, variables, *args,
-                           **kwargs) -> numbers.Real:
+    def sample_paulistring(self, samples: int, paulistring, variables, initial_state: Union[int, QubitWaveFunction] = 0,
+                           *args, **kwargs) -> numbers.Real:
         """
         Sample an individual pauli word (pauli string) and return the average result thereof.
         Parameters
@@ -520,8 +523,8 @@ def sample_paulistring(self, samples: int, paulistring, variables, *args,
         # on construction: tq.ExpectationValue(H=H, U=U, optimize_measurements=True)
         circuit = self.create_circuit(circuit=copy.deepcopy(self.circuit), abstract_circuit=basis_change)
         # run simulators
-        counts = self.sample(samples=samples, circuit=circuit, read_out_qubits=qubits, variables=variables, *args,
-                             **kwargs)
+        counts = self.sample(samples=samples, circuit=circuit, read_out_qubits=qubits, variables=variables,
+                             initial_state=initial_state, *args, **kwargs)
         # compute energy
         E = 0.0
         n_samples = 0
@@ -792,7 +795,7 @@ def __copy__(self):
     def __deepcopy__(self, memodict={}):
         return type(self)(self.abstract_expectationvalue, **self._input_args)
 
-    def __call__(self, variables, samples: int = None, *args, **kwargs):
+    def __call__(self, variables, samples: int = None, initial_state: Union[int, QubitWaveFunction] = 0, *args, **kwargs):
 
         variables = format_variable_dictionary(variables=variables)
         if self._variables is not None and len(self._variables) > 0:
@@ -802,9 +805,9 @@ def __call__(self, variables, samples: int = None, *args, **kwargs):
                         self._variables, variables))
 
         if samples is None:
-            data = self.simulate(variables=variables, *args, **kwargs)
+            data = self.simulate(variables=variables, initial_state=initial_state, *args, **kwargs)
         else:
-            data = self.sample(variables=variables, samples=samples, *args, **kwargs)
+            data = self.sample(variables=variables, samples=samples, initial_state=initial_state, *args, **kwargs)
 
         if self._shape is None and self._contraction is None:
             # this is the default
@@ -852,7 +855,7 @@ def update_variables(self, variables):
         """wrapper over circuit update_variables"""
         self._U.update_variables(variables=variables)
 
-    def sample(self, variables, samples, *args, **kwargs) -> numpy.array:
+    def sample(self, variables, samples, initial_state: Union[int, QubitWaveFunction] = 0, *args, **kwargs) -> numpy.array:
         """
         sample the expectationvalue.
 
@@ -862,6 +865,8 @@ def sample(self, variables, samples, *args, **kwargs) -> numpy.array:
             variables to supply to the unitary.
         samples: int:
             number of samples to perform.
+        initial_state: int or QubitWaveFunction:
+            the initial state of the circuit
         args
         kwargs
 
@@ -891,23 +896,25 @@ def sample(self, variables, samples, *args, **kwargs) -> numpy.array:
             if len(H.qubits) == 0:
                 E = sum([ps.coeff for ps in H.paulistrings])
             elif H.is_all_z():
-                E = self.U.sample_all_z_hamiltonian(samples=samples, hamiltonian=H, variables=variables, *args,
-                                                    **kwargs)
+                E = self.U.sample_all_z_hamiltonian(samples=samples, hamiltonian=H, variables=variables, initial_state=initial_state,
+                                                *args, **kwargs)
             else:
                 for ps in H.paulistrings:
-                    E += self.U.sample_paulistring(samples=samples, paulistring=ps, variables=variables, *args,
-                                                   **kwargs)
+                    E += self.U.sample_paulistring(samples=samples, paulistring=ps, variables=variables, initial_state=initial_state,
+                                                   *args, **kwargs)
             result.append(to_float(E))
         return numpy.asarray(result)
 
-    def simulate(self, variables, *args, **kwargs):
+    def simulate(self, variables, initial_state: Union[int, QubitWaveFunction], *args, **kwargs):
         """
         Simulate the expectationvalue.
 
         Parameters
         ----------
         variables:
             variables to supply to the unitary.
+        initial_state: int or QubitWaveFunction:
+            the initial state of the circuit
         args
         kwargs
 
@@ -922,7 +929,7 @@ def simulate(self, variables, *args, **kwargs):
             final_E = 0.0
             # TODO inefficient,
             # Always better to overwrite this function
-            wfn = self.U.simulate(variables=variables, *args, **kwargs)
+            wfn = self.U.simulate(variables=variables, initial_state=initial_state, *args, **kwargs)
             final_E += wfn.compute_expectationvalue(operator=H)
             result.append(to_float(final_E))
         return numpy.asarray(result)

src/tequila/simulators/simulator_qulacs.py

@@ -429,13 +429,15 @@ class BackendExpectationValueQulacs(BackendExpectationValue):
     use_mapping = True
     BackendCircuitType = BackendCircuitQulacs
 
-    def simulate(self, variables, *args, **kwargs) -> numpy.array:
+    def simulate(self, variables, initial_state: Union[int, QubitWaveFunction] = 0, *args, **kwargs) -> numpy.array:
         """
         Perform simulation of this expectationvalue.
         Parameters
         ----------
         variables:
             variables, to be supplied to the underlying circuit.
+        initial_state: int or QubitWaveFunction:
+            the initial state of the circuit
         args
         kwargs
 
@@ -453,7 +455,7 @@ def simulate(self, variables, *args, **kwargs) -> numpy.array:
             return numpy.asarray[self.H]
 
         self.U.update_variables(variables)
-        state = self.U.initialize_state(self.n_qubits)
+        state = self.U.initialize_state(self.n_qubits, initial_state)
         self.U.circuit.update_quantum_state(state)
         result = []
         for H in self.H:
@@ -495,7 +497,7 @@ def initialize_hamiltonian(self, hamiltonians):
             result.append(qulacs_H)
         return result
 
-    def sample(self, variables, samples, *args, **kwargs) -> numpy.array:
+    def sample(self, variables, samples, initial_state: Union[int, QubitWaveFunction] = 0, *args, **kwargs) -> numpy.array:
         """
         Sample this Expectation Value.
         Parameters
@@ -504,6 +506,8 @@ def sample(self, variables, samples, *args, **kwargs) -> numpy.array:
             variables, to supply to the underlying circuit.
         samples: int:
             the number of samples to take.
+        initial_state: int or QubitWaveFunction:
+            the initial state of the circuit
         args
         kwargs
 
@@ -513,13 +517,13 @@ def sample(self, variables, samples, *args, **kwargs) -> numpy.array:
             the result of sampling as a number.
         """
         self.update_variables(variables)
-        state = self.U.initialize_state(self.n_qubits)
+        state = self.U.initialize_state(self.n_qubits, initial_state)
         self.U.circuit.update_quantum_state(state)
         result = []
-        for H in self._reduced_hamiltonians: # those are the hamiltonians which where non-used qubits are already traced out
+        for H in self._reduced_hamiltonians:  # those are the hamiltonians which where non-used qubits are already traced out
             E = 0.0
             if H.is_all_z() and not self.U.has_noise:
-                E = super().sample(samples=samples, variables=variables, *args, **kwargs)
+                E = super().sample(samples=samples, variables=variables, initial_state=initial_state, *args, **kwargs)
             else:
                 for ps in H.paulistrings:
                     # change basis, measurement is destructive so the state will be copied
@@ -530,8 +534,8 @@ def sample(self, variables, samples, *args, **kwargs) -> numpy.array:
                     qbc = self.U.create_circuit(abstract_circuit=bc, variables=None)
                     Esamples = []
                     for sample in range(samples):
-                        if self.U.has_noise and sample>0:
-                            state = self.U.initialize_state(self.n_qubits)
+                        if self.U.has_noise and sample > 0:
+                            state = self.U.initialize_state(self.n_qubits, initial_state)
                             self.U.circuit.update_quantum_state(state)
                             state_tmp = state
                         else:
@@ -540,7 +544,7 @@ def sample(self, variables, samples, *args, **kwargs) -> numpy.array:
                             qbc.update_quantum_state(state_tmp)
                         ps_measure = 1.0
                         for idx in ps.keys():
-                            assert idx in self.U.abstract_qubits # assert that the hamiltonian was really reduced
+                            assert idx in self.U.abstract_qubits  # assert that the hamiltonian was really reduced
                             M = qulacs.gate.Measurement(self.U.qubit(idx), self.U.qubit(idx))
                             M.update_quantum_state(state_tmp)
                             measured = state_tmp.get_classical_value(self.U.qubit(idx))

tests/test_simulator_backends.py

@@ -371,9 +371,20 @@ def test_initial_state_from_wavefunction(simulator):
     assert result.isclose(QubitWaveFunction.from_array(np.array([100.0, 0.0])))
 
     state = QubitWaveFunction.from_array(np.array([1.0, -1.0])).normalize()
+    result = tq.simulate(U, initial_state=state, backend=simulator)
+    assert result.isclose(QubitWaveFunction.from_basis_state(n_qubits=1, basis_state=1))
     result = tq.simulate(U, initial_state=state, backend=simulator, samples=100)
     assert result.isclose(QubitWaveFunction.from_array(np.array([0.0, 100.0])))
 
+    U = tq.gates.X(target=0)
+    H = tq.paulis.Z(qubit=0)
+    E = tq.ExpectationValue(U, H)
+    state = QubitWaveFunction.from_array(np.array([0.0, 1.0])).normalize()
+    result = tq.simulate(E, initial_state=state, backend=simulator)
+    assert numpy.isclose(result, 1.0)
+    result = tq.simulate(E, initial_state=state, backend=simulator, samples=100)
+    assert numpy.isclose(result, 1.0)
+
 
 
 @pytest.mark.parametrize("backend", tequila.simulators.simulator_api.INSTALLED_SIMULATORS.keys())