V1.9.9 merge (#383)

* Add initial state support for expectation values (#373)

* Improve circuit compilation performance (#375)

* Wave function compact visualization (#376)

* wfn compact visualization

* minor edit

* Update src/tequila/wavefunction/qubit_wavefunction.py

Co-authored-by: Oliver Hüttenhofer <[email protected]>

* Update src/tequila/wavefunction/qubit_wavefunction.py

Co-authored-by: Oliver Hüttenhofer <[email protected]>

* Update test_chemistry.py

for h2 hcb is a bit slower (overhead due to changes in orbopt function, not dramatic)

---------

Co-authored-by: Oliver Hüttenhofer <[email protected]>

* spex backend (#379)

* adding spex as a supported backend
* simulator_spex.py
* simple test case
* using ExpPauliTerm in circuit
* cashing using a circuit hash
* using SPEX and OMP_NUM_THREADS for parallelisation if available
* using a threshold to eliminate elements with negligible amplitudes
* compress_qubit_indicies and moving thresholding to c++
* deleting test_spex_simulator
* Update ci_backends.yml
* num_qubits param for spex backend (from LSB to MSB)

* minor fixup spex backend (#381)

* import error in spex backend
* spex-backend: error in hadamard compiling
* test_gradients: samplers/simulators missmatch

* Bug in orbital_optimizer.py to vqe_solver (#377)

* Bug in orbital_optimizer.py to vqe_solver

* Remove Qiskit noise feature and transpile circuit (#378)

* geometry_to_xyz (#382)

* Update version.py

* Update version.py

---------

Co-authored-by: Oliver Hüttenhofer <[email protected]>
Co-authored-by: davibinco <[email protected]>
Co-authored-by: Mikel-Ma <[email protected]>
Co-authored-by: JdelArco98 <[email protected]>

Author: 5 people

Diff for: .github/workflows/ci_backends.yml

@@ -16,25 +16,27 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ['3.9']
-
+        python-version: ['3.10']
+        include:
+          - os: ubuntu-latest
+            cxx: /usr/bin/g++-14
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python ${{ matrix.python-version }}
       uses: actions/setup-python@v1
       with:
         python-version: ${{ matrix.python-version }}
+    - name: Install spex
+      run: |
+        python -m pip install --upgrade pip
+        pip install git+https://github.com/Mikel-Ma/spex@devel
     - name: Install basic dependencies
       run: |
         python -m pip install --upgrade pip
         pip install -r requirements.txt
-        if [ $ver -eq 7 ]; then
-            # myqlm does not work in python3.7
-            pip install "cirq" "qiskit>=0.30" "qulacs" "qibo==0.1.1"
-        elif [ $ver -eq 8 ]; then
-            pip install "cirq" "qiskit" "qulacs" "qibo==0.1.1" "myqlm" "cirq-google"
-        fi
-        pip install -e .
+        pip install "qiskit" "qulacs"
+        # pip install cirq # installs too much stuff currently
+        pip install -e . 
     - name: Lint with flake8
       run: |
         pip install flake8

Diff for: src/tequila/circuit/circuit.py

@@ -9,6 +9,7 @@
 
 from .qpic import export_to
 
+
 class QCircuit():
     """
     Fundamental class representing an abstract circuit.
@@ -254,8 +255,8 @@ def replace_gates(self, positions: list, circuits: list, replace: list = None):
             the gates to add at the corresponding positions
         replace: list of bool: (Default value: None)
             Default is None which corresponds to all true
-            decide if gates shall be replaces or if the new parts shall be inserted without replacement
-            if replace[i] = true: gate at position [i] will be replaces by gates[i]
+            decide if gates shall be replaced or if the new parts shall be inserted without replacement
+            if replace[i] = true: gate at position [i] will be replaced by gates[i]
             if replace[i] = false: gates[i] circuit will be inserted at position [i] (beaming before gate previously at position [i])
         Returns
         -------
@@ -271,8 +272,9 @@ def replace_gates(self, positions: list, circuits: list, replace: list = None):
         dataset = zip(positions, circuits, replace)
         dataset = sorted(dataset, key=lambda x: x[0])
 
-        offset = 0
-        new_gatelist = self.gates
+        new_gatelist = []
+        last_idx = -1
+
         for idx, circuit, do_replace in dataset:
 
             # failsafe
@@ -283,13 +285,14 @@ def replace_gates(self, positions: list, circuits: list, replace: list = None):
             else:
                 gatelist = [circuit]
 
-            pos = idx + offset
-            if do_replace:
-                new_gatelist = new_gatelist[:pos] + gatelist + new_gatelist[pos + 1:]
-                offset += len(gatelist) - 1
-            else:
-                new_gatelist = new_gatelist[:pos] + gatelist + new_gatelist[pos:]
-                offset += len(gatelist)
+            new_gatelist += self.gates[last_idx + 1:idx]
+            new_gatelist += gatelist
+            if not do_replace:
+                new_gatelist.append(self.gates[idx])
+
+            last_idx = idx
+
+        new_gatelist += self.gates[last_idx + 1:]
 
         result = QCircuit(gates=new_gatelist)
         result.n_qubits = max(result.n_qubits, self.n_qubits)
@@ -386,7 +389,7 @@ def __iadd__(self, other):
         for k, v in other._parameter_map.items():
             self._parameter_map[k] += [(x[0] + offset, x[1]) for x in v]
 
-        self._gates += copy.deepcopy(other.gates)
+        self._gates += other.gates
         self._min_n_qubits = max(self._min_n_qubits, other._min_n_qubits)
 
         return self
@@ -593,7 +596,7 @@ def _inpl_control_circ(self, control) -> None:
 
         This is an in-place method, so it mutates self and doesn't return any value.
 
-        Raise TequilaWarning if there any qubits in common between self and control.
+        Raise TequilaWarning if there are any qubits in common between self and control.
         """
         gates = self.gates
         control = list_assignment(control)
@@ -610,8 +613,10 @@ def _inpl_control_circ(self, control) -> None:
                 if len(control_lst) < len(gate.control) + len(control):
                     # warnings.warn("Some of the controls {} were already included in the control "
                     #               "of a gate {}.".format(control, gate), TequilaWarning)
-                    raise TequilaWarning(f'Some of the controls {control} were already included '
-                                         f'in the control of a gate {gate}.')
+                    raise TequilaWarning(f"Some of the controls {control} were already included "
+                                         f"in the control of a gate {gate}. "
+                                         f"This might be because the same instance of a gate is used in multiple places, "
+                                         f"e.g. because the same circuit was appended twice.")
             else:
                 control_lst, not_control = list(control), list()
 

Diff for: src/tequila/objective/objective.py

@@ -1,4 +1,5 @@
 import typing, copy, numbers
+from typing import Union
 from tequila.grouping.compile_groups import compile_commuting_parts
 from tequila import TequilaException
 from tequila.utils import JoinedTransformation
@@ -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 = 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]

Diff for: src/tequila/quantumchemistry/chemistry_tools.py

@@ -442,7 +442,14 @@ def read_xyz_from_file(filename):
             for i in range(natoms):
                 coord += content[2 + i]
             return coord, comment
-
+    def get_xyz(self)->str:
+        geom = self.parameters.get_geometry()
+        f = ''
+        f += f'{len(geom)}\n'
+        f += f'{self.parameters.name}\n'
+        for at in geom:
+            f += f'{at[0]} {at[1][0]} {at[1][1]} {at[1][2]}\n'
+        return f
 
 @dataclass
 class ClosedShellAmplitudes:

Diff for: src/tequila/quantumchemistry/orbital_optimizer.py

@@ -217,7 +217,8 @@ def kernel(self, h1, h2, *args, **kwargs):
                 vqe_solver_arguments = self.vqe_solver_arguments
             result = self.vqe_solver(H=H, circuit=self.circuit, molecule=molecule, **vqe_solver_arguments)
             if hasattr(self.vqe_solver, "compute_rdms"):
-                rdm1, rdm2 = self.vqe_solver.compute_rdms(U=self.circuit, variables=result.variables, molecule=molecule, use_hcb=restrict_to_hcb)
+                rdm1,rdm2 = self.vqe_solver.compute_rdms(U=self.circuit, variables=result.variables, molecule=molecule, use_hcb=restrict_to_hcb)
+                rdm2 = self.reorder(rdm2, 'dirac', 'mulliken')
         elif self.circuit is None:
             raise Exception("Orbital Optimizer: Either provide a callable vqe_solver or a circuit")
         else:

Diff for: src/tequila/simulators/simulator_api.py

@@ -11,8 +11,9 @@
 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
+SUPPORTED_BACKENDS = ["qulacs", "qulacs_gpu", "qibo", "qiskit", "qiskit_gpu", "cirq", "pyquil", "symbolic", "qlm", "spex"]
+# TODO: Reenable noise for Qiskit
+SUPPORTED_NOISE_BACKENDS = ["cirq", "pyquil"]  # qulacs removed in v.1.9
 BackendTypes = namedtuple('BackendTypes', 'CircType ExpValueType')
 INSTALLED_SIMULATORS = {}
 INSTALLED_SAMPLERS = {}
@@ -30,6 +31,15 @@
 """
 
 
+HAS_SPEX = True
+try:
+    from tequila.simulators.simulator_spex import BackendCircuitSpex, BackendExpectationValueSpex
+
+    INSTALLED_SIMULATORS["spex"] = BackendTypes(BackendCircuitSpex, BackendExpectationValueSpex)
+except ImportError:
+    HAS_SPEX = False
+
+
 HAS_QISKIT = True
 try:
     from tequila.simulators.simulator_qiskit import BackendCircuitQiskit, BackendExpectationValueQiskit

Diff for: 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)

Diff for: src/tequila/simulators/simulator_qiskit.py

@@ -9,7 +9,7 @@
 import warnings
 import numpy as np
 import qiskit, qiskit_aer, qiskit.providers.fake_provider
-from qiskit import QuantumCircuit
+from qiskit import QuantumCircuit, transpile
 
 HAS_NOISE = True
 try:
@@ -316,11 +316,9 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
             optimization_level = kwargs['optimization_level']
 
         circuit = self.circuit.assign_parameters(self.resolver)
-
         circuit = self.add_state_init(circuit, initial_state)
-
         circuit.save_statevector()
-
+        circuit = transpile(circuit, qiskit_backend)
         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)
@@ -349,11 +347,11 @@ def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubit
         if 'optimization_level' in kwargs:
             optimization_level = kwargs['optimization_level']
         if self.device is None:
-            qiskit_backend = self.retrieve_device('aer_simulator')
+            qiskit_backend = self.retrieve_device(self.STATEVECTOR_DEVICE_NAME)
         else:
             qiskit_backend = self.retrieve_device(self.device)
 
-        if isinstance(qiskit_backend, IBMBackend):
+        if HAS_IBMQ and isinstance(qiskit_backend, IBMBackend):
             if self.noise_model is not None:
                 raise TequilaException(
                     'Cannot combine backend {} with custom noise models.'.format(str(qiskit_backend)))
@@ -370,7 +368,7 @@ def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubit
                 if self.noise_model is not None:
                     from_back = self.noise_model
                 basis = from_back.basis_gates
-                use_backend = self.retrieve_device('aer_simulator')
+                use_backend = self.retrieve_device(self.STATEVECTOR_DEVICE_NAME)
                 use_backend.set_options(noise_model=from_back)
                 circuit = qiskit.transpile(circuit, backend=use_backend,
                                            basis_gates=basis,