Density matrix class and simulation

src/tequila/simulators/simulator_api.py

@@ -176,9 +176,9 @@ def pick_backend(backend: str = None, samples: int = None, noise: NoiseModel = N
                 for f in INSTALLED_SAMPLERS.keys():
                     return f
         else:
-            if samples is None:
-                raise TequilaException(
-                    "Noise requires sampling; please provide a positive, integer value for samples")
+            #if samples is None:
+            #    raise TequilaException(
+            #        "Noise requires sampling; please provide a positive, integer value for samples")
             for f in SUPPORTED_NOISE_BACKENDS:
                 return f
             raise TequilaException(

src/tequila/simulators/simulator_base.py

@@ -3,6 +3,7 @@
 from tequila.utils.keymap import KeyMapSubregisterToRegister
 from tequila.utils.misc import to_float
 from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
+from tequila.wavefunction.density_matrix import DensityMatrix
 from tequila.circuit.compiler import change_basis
 from tequila import BitString
 from tequila.objective.objective import Variable, format_variable_dictionary
@@ -367,6 +368,44 @@ def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunc
         result.apply_keymap(keymap=keymap, initial_state=initial_state)
         return result
 
+    def simulate_density(self, variables, initial_state=0, *args, **kwargs) -> DensityMatrix:
+        """
+        simulate the circuit via the backend.
+
+        Parameters
+        ----------
+        variables:
+            the parameters with which to simulate the circuit.
+        initial_state: Default = 0:
+            one of several types; determines the base state onto which the circuit is applied.
+            Default: the circuit is applied to the all-zero state.
+        args
+        kwargs
+
+        Returns
+        -------
+        DensityMatrix
+            the density of the system produced by the action of the circuit on the initial state and noise, if provided.
+        """
+        self.update_variables(variables)
+        if isinstance(initial_state, BitString):
+            initial_state = initial_state.integer
+        if isinstance(initial_state, QubitWaveFunction):
+            if len(initial_state.keys()) != 1:
+                raise TequilaException("only product states as initial states accepted as of now") # TODO: add initial density state for simulation. Can use qiskit quantum info .DensityMatrix.evolve
+            initial_state = list(initial_state.keys())[0].integer
+
+        all_qubits = [i for i in range(self.abstract_circuit.n_qubits)]
+        active_qubits = self.qubit_map.keys()
+
+        # maps from reduced register to full register
+        keymap = KeyMapSubregisterToRegister(subregister=active_qubits, register=all_qubits)
+
+        result = self.do_simulate_density(variables=variables, initial_state=keymap.inverted(initial_state).integer, *args,
+                                  **kwargs)
+        result.apply_keymap(keymap=keymap, initial_state=initial_state)
+        return result
+
     def sample(self, variables, samples, read_out_qubits=None, circuit=None, *args, **kwargs):
         """
         Sample the circuit. If circuit natively equips paulistrings, sample therefrom.

src/tequila/simulators/simulator_qiskit.py

@@ -1,5 +1,6 @@
 from tequila.simulators.simulator_base import BackendCircuit, QCircuit, BackendExpectationValue
 from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
+from tequila.wavefunction.density_matrix import DensityMatrix
 from tequila import TequilaException, TequilaWarning
 from tequila import BitString, BitNumbering, BitStringLSB
 from tequila.utils.keymap import KeyMapRegisterToSubregister
@@ -256,6 +257,8 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
         QubitWaveFunction:
             the result of simulation.
         """
+        circuit = self.circuit.bind_parameters(self.resolver)
+
         if self.noise_model is None:
             if self.device is None:
                 qiskit_backend = self.retrieve_device('statevector_simulator')
@@ -265,7 +268,12 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
                 else:
                     qiskit_backend = self.retrieve_device(self.device)
         else:
-            raise TequilaQiskitException("wave function simulation with noise cannot be performed presently.")
+            ## density simulation in the presence of noise model
+            #if self.device is None:
+            #    qiskit_backend = self.retrieve_device('aer_simulator_density_matrix') #try statevector simulator too
+            raise TequilaQiskitException("wave function simulation with noise cannot be performed, try simulate_density() or sample() instead.")
+            #qiskit_backend.set_options(noise_model=self.noise_model)
+            #circuit.save_density_matrix()
 
         optimization_level = None
         if "optimization_level" in kwargs:
@@ -278,13 +286,61 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
             print(initial_state, " -> ", i)
             array[i.integer] = 1.0
             opts = {"initial_statevector": array}
+
+
+        qiskit_job = qiskit_backend.run(circuit,optimization_level=optimization_level,**opts)
+        backend_result = qiskit_job.result()
+
+        return QubitWaveFunction.from_array(arr=backend_result.get_statevector(circuit), numbering=self.numbering)
+
+    def do_simulate_density(self, variables, initial_state=0, *args, **kwargs) -> DensityMatrix:
+        """
+        Helper function for performing Density simulation.
+        Parameters
+        ----------
+        variables:
+            variables to pass to the circuit for simulation.
+        initial_state:
+            indicate initial state on which the unitary self.circuit should act.
+        args
+        kwargs
 
+        Returns
+        -------
+        DensityMatrix:
+            the result of simulation.
+        """
         circuit = self.circuit.bind_parameters(self.resolver)
 
-        qiskit_job = qiskit_backend.run(circuit,optimization_level=optimization_level,**opts)
+        if self.device is None:
+            qiskit_backend = self.retrieve_device('aer_simulator_density_matrix')
+        else:
+            if 'aer_simulator_density_matrix' not in str(self.device):
+                raise TequilaException('Density simulation with qiskit requires Qiskit Aer density simulator.')
+            else:
+                qiskit_backend = self.retrieve_device(self.device)
+
+        qiskit_backend.set_options(noise_model=self.noise_model)
+        circuit.save_density_matrix()
 
+        optimization_level = None
+        if "optimization_level" in kwargs:
+            optimization_level = kwargs['optimization_level']
+
+        #inital state stuff to add depending on whether it works
+        opts={}
+        if initial_state != 0:
+            array = numpy.zeros(shape=[2 ** self.n_qubits])
+            i = BitStringLSB.from_binary(BitString.from_int(integer=initial_state, nbits=self.n_qubits).binary)
+            print(initial_state, " -> ", i)
+            array[i.integer] = 1.0
+            opts = {"initial_statevector": array}
+
+
+        qiskit_job = qiskit_backend.run(circuit,optimization_level=optimization_level,**opts)
         backend_result = qiskit_job.result()
-        return QubitWaveFunction.from_array(arr=backend_result.get_statevector(circuit), numbering=self.numbering)
+
+        return DensityMatrix.from_array(density_matrix = backend_result.data()['density_matrix'].data, numbering = BitNumbering.MSB)
 
     def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubits, *args, **kwargs) -> QubitWaveFunction:
         """
@@ -306,6 +362,10 @@ def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubit
         optimization_level = 1
         if 'optimization_level' in kwargs:
             optimization_level = kwargs['optimization_level']
+        if 'save_runs' in kwargs:
+            save_runs = kwargs['save_runs']
+        else:
+            save_runs = False #default
         if self.device is None:
             qiskit_backend = self.retrieve_device('aer_simulator')
         else:
@@ -317,8 +377,8 @@ def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubit
             circuit = circuit.bind_parameters(self.resolver)  # this is necessary in spite of qiskit "fixing" it
             circuit = qiskit.transpile(circuit, qiskit_backend)
             return self.convert_measurements(qiskit_backend.run(circuit,shots=samples,
-                                                            optimization_level=optimization_level),
-                                             target_qubits=read_out_qubits)
+                                                            optimization_level=optimization_level, memory=save_runs),
+                                             target_qubits=read_out_qubits, save_runs=save_runs)
         else:
             if isinstance(qiskit_backend, qiskit.test.mock.FakeBackend):
                 circuit = circuit.bind_parameters(self.resolver)  # this is necessary in spite of qiskit "fixing" it
@@ -335,8 +395,8 @@ def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubit
                                            optimization_level=optimization_level
                                            )
 
-                job=qiskit_backend.run(circuit, shots=samples)
-                return self.convert_measurements(job,target_qubits=read_out_qubits)
+                job=qiskit_backend.run(circuit, shots=samples, memory=save_runs)
+                return self.convert_measurements(job,target_qubits=read_out_qubits, save_runs=save_runs)
             else:
                 if self.noise_model is not None:
                     qiskit_backend.set_options(noise_model=self.noise_model)  # fits better with our methodology.
@@ -349,11 +409,11 @@ def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubit
                                            optimization_level=optimization_level
                                            )
 
-                job = qiskit_backend.run(circuit, shots=samples)
+                job = qiskit_backend.run(circuit, shots=samples, memory=save_runs)
                 return self.convert_measurements(job,
-                                                 target_qubits=read_out_qubits)
+                                                 target_qubits=read_out_qubits, save_runs=save_runs)
 
-    def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveFunction:
+    def convert_measurements(self, backend_result, target_qubits=None, save_runs=False) -> QubitWaveFunction:
         """
         map backend results to QubitWaveFunction
         Parameters
@@ -365,8 +425,18 @@ def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveF
         QubitWaveFunction:
             measurements converted into wave function form.
         """
+        def process_qiskit_memory(memory):
+            """
+            Returns run list from qiskit.result.get_memory() (list of run readouts in lSB)
+            """
+            runs = [BitString.from_bitstring(other=BitStringLSB.from_binary(binary=k)) for k in memory]
+            return runs
         qiskit_counts = backend_result.result().get_counts()
         result = QubitWaveFunction()
+        if save_runs:
+            qiskit_memory = backend_result.result().get_memory()
+            qiskit_runs = process_qiskit_memory(qiskit_memory)
+            result.runs = qiskit_runs #todo bitstring!
         # todo there are faster ways
         for k, v in qiskit_counts.items():
             converted_key = BitString.from_bitstring(other=BitStringLSB.from_binary(binary=k))
@@ -376,6 +446,9 @@ def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveF
             mapped_full = [self.qubit_map[q].number for q in self.abstract_qubits]
             keymap = KeyMapRegisterToSubregister(subregister=mapped_target, register=mapped_full)
             result = result.apply_keymap(keymap=keymap)
+            if save_runs:
+                #map results.runs
+                result.runs = [keymap(input_state=k, initial_state=None) for k in result.runs]
 
         return result