Merge pull request #364 from tequilahub/update-to-1.9.7

Update to 1.9.7

requirements.txt

@@ -1,6 +1,6 @@
 # necessary
 numpy
-scipy < 1.11 # for now, until issues with scipy/pyscf are fixed # can in principle be smaller, we recommend >= 1.5 for consistency with our tutorials (numerical gradients mostly)
+scipy
 sympy
 #jax
 #jaxlib

src/tequila/circuit/circuit.py

@@ -46,10 +46,7 @@ def export_to(self, *args, **kwargs):
         Convenience: see src/tequila/circuit/qpic.py - export_to for more
         Parameters
         """
-        # this way we allow calling U.export_to("asd.png") instead of having to specify U.export_to(filename="asd.png")
-        if "circuit" not in kwargs:
-            kwargs["circuit"]=self
-        return export_to(*args, **kwargs)
+        return export_to(self, *args, **kwargs)
 
     @property
     def moments(self):
@@ -389,7 +386,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 += other.gates
+        self._gates += copy.deepcopy(other.gates)
         self._min_n_qubits = max(self._min_n_qubits, other._min_n_qubits)
 
         return self

src/tequila/quantumchemistry/qc_base.py

@@ -1737,7 +1737,7 @@ def rdm2(self):
 
     def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_free: bool = True,
                      get_rdm1: bool = True, get_rdm2: bool = True, ordering="dirac", use_hcb: bool = False,
-                     rdm_trafo: QubitHamiltonian = None, decompose=None):
+                     rdm_trafo: QubitHamiltonian = None, evaluate=True):
         """
         Computes the one- and two-particle reduced density matrices (rdm1 and rdm2) given
         a unitary U. This method uses the standard ordering in physics as denoted below.
@@ -1768,7 +1768,10 @@ def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_fre
         rdm_trafo :
             The rdm operators can be transformed, e.g., a^dagger_i a_j -> U^dagger a^dagger_i a_j U,
             where U represents the transformation. The default is set to None, implying that U equas the identity.
-
+        evaluate :
+            if true, the tequila expectation values are evaluated directly via the tq.simulate command.
+            the protocol is optimized to avoid repetation of wavefunction simulation
+            if false, the rdms are returned as tq.QTensors
         Returns
         -------
         """
@@ -1891,13 +1894,14 @@ def _build_2bdy_operators_spinfree() -> list:
                     ops += [op]
             return ops
 
-        def _assemble_rdm1(evals) -> numpy.ndarray:
+        def _assemble_rdm1(evals, rdm1=None) -> numpy.ndarray:
             """
             Returns spin-ful or spin-free one-particle RDM built by symmetry conditions
             Same symmetry with or without spin, so we can use the same function
             """
             N = n_MOs if spin_free else n_SOs
-            rdm1 = numpy.zeros([N, N])
+            if rdm1 is None:
+                rdm1 = numpy.zeros([N, N])
             ctr: int = 0
             for p in range(N):
                 for q in range(p + 1):
@@ -1908,10 +1912,11 @@ def _assemble_rdm1(evals) -> numpy.ndarray:
 
             return rdm1
 
-        def _assemble_rdm2_spinful(evals) -> numpy.ndarray:
+        def _assemble_rdm2_spinful(evals, rdm2=None) -> numpy.ndarray:
             """ Returns spin-ful two-particle RDM built by symmetry conditions """
             ctr: int = 0
-            rdm2 = numpy.zeros([n_SOs, n_SOs, n_SOs, n_SOs])
+            if rdm2 is None:
+                rdm2 = numpy.zeros([n_SOs, n_SOs, n_SOs, n_SOs])
             for p in range(n_SOs):
                 for q in range(p):
                     for r in range(n_SOs):
@@ -1933,10 +1938,11 @@ def _assemble_rdm2_spinful(evals) -> numpy.ndarray:
 
             return rdm2
 
-        def _assemble_rdm2_spinfree(evals) -> numpy.ndarray:
+        def _assemble_rdm2_spinfree(evals, rdm2=None) -> numpy.ndarray:
             """ Returns spin-free two-particle RDM built by symmetry conditions """
             ctr: int = 0
-            rdm2 = numpy.zeros([n_MOs, n_MOs, n_MOs, n_MOs])
+            if rdm2 is None:
+                rdm2 = numpy.zeros([n_MOs, n_MOs, n_MOs, n_MOs])
             for p, q, r, s in product(range(n_MOs), repeat=4):
                 if p * n_MOs + q >= r * n_MOs + s and (p >= q or r >= s):
                     rdm2[p, q, r, s] = evals[ctr]
@@ -2012,18 +2018,25 @@ def _build_2bdy_operators_hcb() -> list:
         # Transform operator lists to QubitHamiltonians
         if (not use_hcb):
             qops = [_get_qop_hermitian(op) for op in qops]
-        
+
         # Compute expected values
-        if rdm_trafo is None:
-            if decompose is not None:
-                print("MANIPULATED")
-                X = decompose(H=qops, U=U)
-                evals = simulate(X, variables=variables)
+        rdm1 = None
+        rdm2 = None
+        from tequila import QTensor
+        if evaluate:
+            if rdm_trafo is None:
+                evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables)
             else:
+                qops = [rdm_trafo.dagger()*qops[i]*rdm_trafo for i in range(len(qops))]
                 evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables)
         else:
-            qops = [rdm_trafo.dagger()*qops[i]*rdm_trafo for i in range(len(qops))]
-            evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables)
+            if rdm_trafo is None:
+                evals = [ExpectationValue(H=x, U=U) for x in qops]
+                N = n_MOs if spin_free else n_SOs
+                rdm1 = QTensor(shape=[N,N])
+                rdm2 = QTensor(shape=[N, N, N, N])
+            else:
+                raise TequilaException("compute_rdms: rdm_trafo was set but evaluate flag is False (not supported)")
 
         # Assemble density matrices
         # If self._rdm1, self._rdm2 exist, reset them if they are of the other spin-type
@@ -2044,11 +2057,11 @@ def _reset_rdm(rdm):
             len_1 = 0
         evals_1, evals_2 = evals[:len_1], evals[len_1:]
         # Build matrices using the expectation values
-        self._rdm1 = _assemble_rdm1(evals_1) if get_rdm1 else self._rdm1
+        self._rdm1 = _assemble_rdm1(evals_1, rdm1=rdm1) if get_rdm1 else self._rdm1
         if spin_free or use_hcb:
-            self._rdm2 = _assemble_rdm2_spinfree(evals_2) if get_rdm2 else self._rdm2
+            self._rdm2 = _assemble_rdm2_spinfree(evals_2, rdm2=rdm2) if get_rdm2 else self._rdm2
         else:
-            self._rdm2 = _assemble_rdm2_spinful(evals_2) if get_rdm2 else self._rdm2
+            self._rdm2 = _assemble_rdm2_spinful(evals_2, rdm2=rdm2) if get_rdm2 else self._rdm2
 
         if get_rdm2:
             rdm2 = NBodyTensor(elems=self.rdm2, ordering="dirac", verify=False)

src/tequila/simulators/simulator_base.py

@@ -356,15 +356,23 @@ def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunc
                 raise TequilaException("only product states as initial states accepted")
             initial_state = list(initial_state.keys())[0].integer
 
-        all_qubits = [i for i in range(self.abstract_circuit.n_qubits)]
+        all_qubits = list(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)
+        # Keymap is only necessary if not all qubits are active
+        keymap_required = sorted(active_qubits) != all_qubits
 
-        result = self.do_simulate(variables=variables, initial_state=keymap.inverted(initial_state).integer, *args,
+        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)
+        result = self.do_simulate(variables=variables, initial_state=mapped_initial_state, *args,
                                   **kwargs)
-        result.apply_keymap(keymap=keymap, initial_state=initial_state)
+
+        if keymap_required:
+            result.apply_keymap(keymap=keymap, initial_state=initial_state)
+
         return result
 
     def sample(self, variables, samples, read_out_qubits=None, circuit=None, *args, **kwargs):

src/tequila/utils/bitstrings.py

@@ -1,6 +1,7 @@
 from enum import Enum
 from typing import List
 from functools import total_ordering
+from math import ceil, log2
 
 
 class BitNumbering(Enum):
@@ -35,7 +36,7 @@ def nbits(self, value):
 
     def update_nbits(self):
         current = self.nbits
-        min_needed = len(format(self._value, 'b'))
+        min_needed = ceil(log2(self._value + 1))
         self._nbits = max(current, min_needed)
         return self
 
@@ -177,15 +178,23 @@ def numbering(self) -> BitNumbering:
         return BitNumbering.LSB
 
 
+def _reverse_int_bits(x: int, nbits: int) -> int:
+    if nbits is None:
+        nbits = x.bit_length()
+    assert nbits <= 32
+
+    x = ((x & 0x55555555) << 1) | ((x & 0xAAAAAAAA) >> 1)
+    x = ((x & 0x33333333) << 2) | ((x & 0xCCCCCCCC) >> 2)
+    x = ((x & 0x0F0F0F0F) << 4) | ((x & 0xF0F0F0F0) >> 4)
+    x = ((x & 0x00FF00FF) << 8) | ((x & 0xFF00FF00) >> 8)
+    x = ((x & 0x0000FFFF) << 16) | ((x & 0xFFFF0000) >> 16)
+    return x >> (32 - nbits)
+
+
 def initialize_bitstring(integer: int, nbits: int = None, numbering_in: BitNumbering = BitNumbering.MSB,
                          numbering_out: BitNumbering = BitNumbering.MSB):
-    if numbering_in == BitNumbering.MSB:
-        if numbering_out == BitNumbering.MSB:
-            return BitString.from_int(integer=integer, nbits=nbits)
-        else:
-            return BitString.from_binary(binary=BitStringLSB.from_int(integer=integer, nbits=nbits).binary, nbits=nbits)
+    integer = _reverse_int_bits(integer, nbits) if numbering_in != numbering_out else integer
+    if numbering_out == BitNumbering.MSB:
+        return BitString.from_int(integer=integer, nbits=nbits)
     else:
-        if numbering_out == BitNumbering.LSB:
-            return BitStringLSB.from_int(integer=integer, nbits=nbits)
-        else:
-            return BitStringLSB.from_binary(binary=BitString.from_int(integer=integer, nbits=nbits).binary, nbits=nbits)
+        return BitStringLSB.from_int(integer=integer, nbits=nbits)

src/tequila/version.py

@@ -1,2 +1,2 @@
-__version__ = "1.9.6"
+__version__ = "1.9.7"
 __author__ = "Tequila Developers "

src/tequila/wavefunction/qubit_wavefunction.py

@@ -145,18 +145,12 @@ def from_array(cls, arr: numpy.ndarray, keymap=None, threshold: float = 1.e-6,
         maxkey = len(arr) - 1
         maxbit = initialize_bitstring(integer=maxkey, numbering_in=numbering, numbering_out=cls.numbering).nbits
         for ii, v in enumerate(arr):
-            i = initialize_bitstring(integer=ii, nbits=maxbit, numbering_in=numbering, numbering_out=cls.numbering)
-            if not numpy.isclose(abs(v), 0.0, atol=threshold):
+            if abs(v) > threshold:
+                i = initialize_bitstring(integer=ii, nbits=maxbit, numbering_in=numbering, numbering_out=cls.numbering)
                 key = i if keymap is None else keymap(i)
                 state[key] = v
         result = QubitWaveFunction(state, n_qubits=n_qubits)
 
-        if cls.numbering != numbering:
-            if cls.numbering == BitNumbering.MSB:
-                result.apply_keymap(keymap=KeyMapLSB2MSB())
-            else:
-                result.apply_keymap(keymap=KeyMapMSB2LSB())
-
         return result
 
     @classmethod