From c57ad618186f1a1c647effba62e489c08ea32bb7 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Thu, 29 Feb 2024 11:08:55 +0100
Subject: [PATCH 1/6] Update to 1.9.4 (#330)

* fixing syntax issue in post_init of dataclass (#327)

* fixing syntax issue in post_init of dataclass
* phoenics dropped due to maintenance resources
* more convenient randomization initialization for OO, avoiding numpy warnings

* Update version.py (#329)
---
 requirements.txt                              |   1 -
 src/tequila/apps/adapt/adapt.py               |   2 +-
 src/tequila/optimizers/__init__.py            |  15 +-
 src/tequila/optimizers/optimizer_gpyopt.py    |   2 -
 src/tequila/optimizers/optimizer_phoenics.py  | 348 ------------------
 .../quantumchemistry/chemistry_tools.py       |   4 +-
 .../quantumchemistry/madness_interface.py     |   2 +-
 .../quantumchemistry/orbital_optimizer.py     |  10 +-
 src/tequila/quantumchemistry/qc_base.py       |   2 +-
 src/tequila/version.py                        |   2 +-
 tests/test_noise_opt.py                       |  14 -
 tests/test_phoenics.py                        |  68 ----
 12 files changed, 13 insertions(+), 457 deletions(-)
 delete mode 100644 src/tequila/optimizers/optimizer_phoenics.py
 delete mode 100644 tests/test_phoenics.py

diff --git a/requirements.txt b/requirements.txt
index 6f9ae975..0dd2e329 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -20,7 +20,6 @@ qulacs # default simulator (best integration), remove if the installation gives
 #qibo <= 0.1.1 # can not be installed in the same environment as gpyopt
 
 #optional optimizers
-#phoenics # version on PyPi  isc urrently broken, we recommend to install from source (AAG github)
 #gpyopt # not in combination with qibo as quantum backend
 
 #optional third party libraries
diff --git a/src/tequila/apps/adapt/adapt.py b/src/tequila/apps/adapt/adapt.py
index 6b8dd9e9..4b9a5549 100644
--- a/src/tequila/apps/adapt/adapt.py
+++ b/src/tequila/apps/adapt/adapt.py
@@ -21,7 +21,7 @@ class AdaptParameters:
     degeneracy_threshold: float = 5.e-4
     silent: bool = False
     
-    def __post__init__(self):
+    def __post_init__(self):
         # avoid stacking of same operator-types in a row
         if "method_options" in self.optimizer_args:
             if "gtol" in self.optimizer_args["method_options"]:
diff --git a/src/tequila/optimizers/__init__.py b/src/tequila/optimizers/__init__.py
index 0022f7de..032606c1 100644
--- a/src/tequila/optimizers/__init__.py
+++ b/src/tequila/optimizers/__init__.py
@@ -16,7 +16,7 @@ class _Optimizers:
     methods: list = None
 
 
-SUPPORTED_OPTIMIZERS = ['scipy', 'phoenics', 'gpyopt', 'gd']
+SUPPORTED_OPTIMIZERS = ['scipy', 'gpyopt', 'gd']
 INSTALLED_OPTIMIZERS = {}
 INSTALLED_OPTIMIZERS['scipy'] = _Optimizers(cls=OptimizerSciPy,
                                             minimize=minimize_scipy,
@@ -37,19 +37,6 @@ class _Optimizers:
 except ImportError:
     has_gpyopt = False
 
-has_phoenics = False
-try:
-    from tequila.optimizers.optimizer_phoenics import OptimizerPhoenics
-    from tequila.optimizers.optimizer_phoenics import minimize as minimize_phoenics
-
-    INSTALLED_OPTIMIZERS['phoenics'] = _Optimizers(cls=OptimizerPhoenics,
-                                                   minimize=minimize_phoenics,
-                                                   methods=OptimizerPhoenics.available_methods())
-    has_phoenics = True
-except ImportError:
-    has_phoenics = False
-
-
 def show_available_optimizers(module=None):
     """
     Returns
diff --git a/src/tequila/optimizers/optimizer_gpyopt.py b/src/tequila/optimizers/optimizer_gpyopt.py
index 50834e07..cfde41b9 100644
--- a/src/tequila/optimizers/optimizer_gpyopt.py
+++ b/src/tequila/optimizers/optimizer_gpyopt.py
@@ -4,8 +4,6 @@
 import numbers
 from tequila.objective.objective import Variable
 import warnings
-
-warnings.simplefilter("ignore")
 import GPyOpt
 from GPyOpt.methods import BayesianOptimization
 import numpy as np
diff --git a/src/tequila/optimizers/optimizer_phoenics.py b/src/tequila/optimizers/optimizer_phoenics.py
deleted file mode 100644
index e2aa77f5..00000000
--- a/src/tequila/optimizers/optimizer_phoenics.py
+++ /dev/null
@@ -1,348 +0,0 @@
-from tequila.objective.objective import Objective
-from tequila.optimizers.optimizer_base import Optimizer, OptimizerResults, dataclass
-import typing
-import numbers
-from tequila.objective.objective import Variable
-import copy
-import warnings
-import pickle
-import time
-from tequila import TequilaException
-warnings.simplefilter("ignore")
-with warnings.catch_warnings():
-    warnings.filterwarnings("ignore", category=DeprecationWarning)
-    warnings.filterwarnings("ignore")
-import phoenics
-
-import numpy as np
-from numpy import pi as pi
-from tequila.simulators.simulator_api import compile_objective
-import os
-
-#numpy, tf, etc can get real, real, real, noisy here. We suppress it.
-warnings.filterwarnings('ignore', category=DeprecationWarning)
-os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'
-warnings.filterwarnings('ignore', category=FutureWarning)
-
-@dataclass
-class PhoenicsResults(OptimizerResults):
-
-    observations: list = None
-    phoenics_instance: phoenics.Phoenics = None
-
-class OptimizerPhoenics(Optimizer):
-    """
-    wrapper to allow optimization of objectives with Phoenics, a bayesian optimizer.
-    See: https://github.com/aspuru-guzik-group/phoenics
-    """
-    @classmethod
-    def available_methods(cls):
-        return ["phoenics"]
-
-    def __init__(self, maxiter, backend=None, save_history=True, minimize=True,
-                 samples=None, silent=None, noise=None, device=None):
-        self._minimize = minimize
-
-        super().__init__(backend=backend, maxiter=maxiter, samples=samples,
-                         noise=noise,device=device,
-                         save_history=save_history, silent=silent)
-
-    def _process_for_sim(self, recommendation, passive_angles):
-        """
-        convert from the phoenics suggestion format to a version recognizable by objectives.
-        Parameters
-        ----------
-        recommendation: dict:
-            the a phoenics suggestion.
-        passive_angles: dict:
-            passive angles not optimized over.
-
-        Returns
-        -------
-        dict:
-            dict of Bariable, float pairs.
-        """
-        rec = copy.deepcopy(recommendation)
-        for part in rec:
-            for k, v in part.items():
-                part[k] = v.item()
-            if passive_angles is not None:
-                for k, v in passive_angles.items():
-                    part[k] = v
-        return rec
-
-    def _process_for_phoenics(self, pset, result, passive_angles=None):
-        """
-        Convert results of a call to an objective into a form interpretable by phoenics.
-        Parameters
-        ----------
-        pset: dict:
-            the parameters evaluated, as a dictionary
-        result:
-            the result of calling some objective, using pset as parameters.
-        passive_angles: dict, optional:
-            passive_angles, not optimized over.
-
-        Returns
-        -------
-        dict:
-            the a dictionary, formatted as phoenics prefers it, for use as an 'observation'.
-        """
-        new = copy.deepcopy(pset)
-        for k, v in new.items():
-            new[k] = np.array([v], dtype=np.float32)
-        if passive_angles is not None:
-            for k in passive_angles.keys():
-                del new[k]
-        new['Energy'] = result
-
-        return new
-
-    def _make_phoenics_object(self, objective, passive_angles=None, conf=None, *args, **kwargs):
-        """
-        instantiate phoenics, to perform optimization.
-
-        Parameters
-        ----------
-        objective: Objective:
-            the objective to optimize over.
-        passive_angles: dict, optional:
-            a dictionary of angles not to optimize over.
-        conf: optional:
-            a user built configuration object or file, from which to initialize a phoenics object.
-            For advanced users only.
-        args
-        kwargs
-
-        Returns
-        -------
-        phoenics.Phoenics
-            a phoenics object configured to optimize an objective.
-        """
-        if conf is not None:
-            if hasattr(conf, 'readlines'):
-                bird = phoenics.Phoenics(config_file=conf)
-            else:
-                bird = phoenics.Phoenics(config_dict=conf)
-
-            return bird
-        op = objective.extract_variables()
-        if passive_angles is not None:
-            for i, thing in enumerate(op):
-                if thing in passive_angles.keys():
-                    op.remove(thing)
-
-        config = {"general": {"auto_desc_gen": "False", "batches": 5, "boosted": "False", "parallel": "False", "scratch_dir":os.getcwd()}}
-        config['parameters'] = [
-            {'name': k, 'periodic': 'True', 'type': 'continuous', 'size': 1, 'low': 0, 'high': 2 * pi} for k in op]
-        if self._minimize is True:
-            config['objectives'] = [{"name": "Energy", "goal": "minimize"}]
-        else:
-            config['objectives'] = [{"name": "Energy", "goal": "maximize"}]
-
-        for k,v in kwargs.items():
-            if hasattr(k, "lower") and k.lower() in config["general"]:
-                config["general"][k.lower()] = v
-
-        bird = phoenics.Phoenics(config_dict=config)
-        return bird
-
-    def __call__(self, objective: Objective,
-                 maxiter=None,
-                 variables: typing.List[Variable] = None,
-                 initial_values: typing.Dict[Variable, numbers.Real] = None,
-                 previous=None,
-                 phoenics_config=None,
-                 file_name=None,
-                 *args,
-                 **kwargs):
-        """
-        Perform optimization with phoenics.
-
-        Parameters
-        ----------
-        objective: Objective
-            the objective to optimize.
-        maxiter: int:
-            (Default value = None)
-            if not None, overwrite the init maxiter with new number.
-        variables: list:
-            (Default value = None)
-            which variables to optimize over. If None: all of the variables in objective are used.
-        initial_values: dict:
-            (Default value = None)
-            an initial point to begin optimization from. Random, if None.
-        previous:
-            previous observations, formatted for phoenics, to use in optimization. For use by advanced users.
-        phoenics_config:
-            a config for a phoenics object.
-        file_name:
-            a file
-        args
-        kwargs
-
-        Returns
-        -------
-        PhoenicsResults:
-            the results of optimization by phoenics.
-
-        """
-
-        objective = objective.contract()
-        active_angles, passive_angles, variables = self.initialize_variables(objective,
-                                                               initial_values=initial_values,
-                                                               variables=variables)
-
-        if maxiter is None:
-            maxiter = 10
-
-        obs = []
-        bird = self._make_phoenics_object(objective, passive_angles, phoenics_config, *args, **kwargs)
-        if previous is not None:
-            if type(previous) is str:
-                try:
-                    obs = pickle.load(open(previous, 'rb'))
-                except:
-                    print(
-                        'failed to load previous observations, which are meant to be a pickle file. Starting fresh.')
-            elif type(previous) is list:
-                if all([type(k) == dict for k in previous]):
-                    obs = previous
-                else:
-                    print('previous observations were not in the correct format (list of dicts). Starting fresh.')
-
-
-
-        if not (type(file_name) == str or file_name == None):
-            raise TequilaException('file_name must be a string, or None. Recieved {}'.format(type(file_name)))
-
-        best = None
-        best_angles = None
-
-        # avoid multiple compilations
-        compiled_objective = compile_objective(objective=objective, backend=self.backend,
-                                               device=self.device,
-                                               samples=self.samples, noise=self.noise)
-
-        if not self.silent:
-            print('phoenics has recieved')
-            print("objective: \n")
-            print(objective)
-            print("noise model : {}".format(self.noise))
-            print("samples     : {}".format(self.samples))
-            print("maxiter     : {}".format(maxiter))
-            print("variables   : {}".format(objective.extract_variables()))
-            print("passive var : {}".format(passive_angles))
-            print('now lets begin')
-        for i in range(0, maxiter):
-            with warnings.catch_warnings():
-                np.testing.suppress_warnings()
-                warnings.simplefilter("ignore")
-                warnings.filterwarnings("ignore", category=FutureWarning)
-
-            precs = bird.recommend(observations=obs)
-
-            runs = []
-            recs = self._process_for_sim(precs, passive_angles=passive_angles)
-
-            start = time.time()
-            for j, rec in enumerate(recs):
-                En = compiled_objective(variables=rec, samples=self.samples, noise=self.noise)
-                runs.append((rec, En))
-                if not self.silent:
-                    if self.print_level > 2:
-                        print("energy = {:+2.8f} , angles=".format(En), rec)
-                    else:
-                        print("energy = {:+2.8f}".format(En))
-            stop = time.time()
-            if not self.silent:
-                print("Quantum Objective evaluations: {}s Wall-Time".format(stop-start))
-
-            for run in runs:
-                angles = run[0]
-                E = run[1]
-                if best is None:
-                    best = E
-                    best_angles = angles
-                else:
-                    if self._minimize:
-                        if E < best:
-                            best = E
-                            best_angles = angles
-                    else:
-                        if E > best:
-                            best = E
-                            best_angles = angles
-
-                if self.save_history:
-                    self.history.energies.append(E)
-                    self.history.angles.append(angles)
-                obs.append(self._process_for_phoenics(angles, E, passive_angles=passive_angles))
-
-        if file_name is not None:
-            with open(file_name, 'wb') as file:
-                pickle.dump(obs, file)
-
-        if not self.silent:
-            print("best energy after {} iterations : {:+2.8f}".format(self.maxiter, best))
-        return PhoenicsResults(energy=best, variables=best_angles, history=self.history, observations=obs,phoenics_instance=bird)
-
-
-def minimize(objective: Objective,
-             maxiter: int = None,
-             samples: int = None,
-             variables: typing.List = None,
-             initial_values: typing.Dict = None,
-             backend: str = None,
-             noise=None,
-             device: str = None,
-             previous: typing.Union[str, list] = None,
-             phoenics_config: typing.Union[str, typing.Dict] = None,
-             file_name: str = None,
-             silent: bool = False,
-             *args,
-             **kwargs) -> PhoenicsResults:
-    """
-
-    Parameters
-    ----------
-    objective: Objective:
-        The tequila objective to optimize
-    initial_values: typing.Dict[typing.Hashable, numbers.Real], optional:
-        Initial values as dictionary of Hashable types (variable keys) and floating point numbers.
-        If given None they will be randomized.
-    variables: typing.List[typing.Hashable], optional:
-         List of Variables to optimize
-    samples: int, optional:
-         samples/shots to take in every run of the quantum circuits (None activates full wavefunction simulation)
-    maxiter: int:
-         how many iterations of phoenics to run.
-         Note that this is NOT identical to the number of times the circuit will run.
-    backend: str, optional:
-         Simulator backend, will be automatically chosen if set to None
-    noise: NoiseModel, optional:
-         a noise model to apply to the circuits of Objective.
-    device: optional:
-        the device from which to (potentially, simulatedly) sample all quantum circuits employed in optimization.
-    previous: optional:
-        Previous phoenics observations. If string, the name of a file from which to load them. Else, a list.
-    phoenics_config: optional:
-        a pre-made phoenics configuration. if str, the name of a file from which to load it; Else, a dictionary.
-        Individual keywords of the 'general' sections can also be passed down as kwargs
-    file_name: str, optional:
-        where to save output to, if save_to_file is True.
-    kwargs: dict:
-        Send down more keywords for single replacements in the phoenics config 'general' section, like e.g. batches=5,
-        boosted=True etc
-    Returns
-    -------
-    PhoenicsResults:
-        the result of an optimization by phoenics.
-    """
-
-    optimizer = OptimizerPhoenics(samples=samples, backend=backend,
-                                  noise=noise,device=device,
-                                  maxiter=maxiter, silent=silent)
-    return optimizer(objective=objective, initial_values=initial_values, variables=variables, previous=previous,
-                     maxiter=maxiter,
-                     phoenics_config=phoenics_config, file_name=file_name, *args, **kwargs)
diff --git a/src/tequila/quantumchemistry/chemistry_tools.py b/src/tequila/quantumchemistry/chemistry_tools.py
index 4d9b2e79..ff1f3ff1 100644
--- a/src/tequila/quantumchemistry/chemistry_tools.py
+++ b/src/tequila/quantumchemistry/chemistry_tools.py
@@ -569,7 +569,7 @@ def _verify_ordering_of(self, trials=100):
         return True
 
     def __init__(self, elems: numpy.ndarray = None, active_indices: list = None, ordering: str = None,
-                 size_full: int = None):
+                 size_full: int = None, verify=False):
         """
         Parameters
         ----------
@@ -611,7 +611,7 @@ def __init__(self, elems: numpy.ndarray = None, active_indices: list = None, ord
         if self.order == 4:
             if ordering is None:
                 ordering = self.identify_ordering()
-            else:
+            elif verify:
                 try: # some RDMs are really sloppy (depends on backend)
                     auto_ordering=self.identify_ordering()
                     if auto_ordering is not ordering:
diff --git a/src/tequila/quantumchemistry/madness_interface.py b/src/tequila/quantumchemistry/madness_interface.py
index e56d58e3..9273046c 100644
--- a/src/tequila/quantumchemistry/madness_interface.py
+++ b/src/tequila/quantumchemistry/madness_interface.py
@@ -122,7 +122,7 @@ def __init__(self, parameters: ParametersQC,
             h = "failed"
             g = "failed"
 
-        if "failed" in h or "failed" in g:
+        if (isinstance(h, str) and "failed" in h) or (isinstance(g, str) and "failed" in g):
             status = "found {}_htensor.npy={}\n".format(name, "failed" not in h)
             status += "found {}_gtensor.npy={}\n".format(name, "failed" not in g)
             try:
diff --git a/src/tequila/quantumchemistry/orbital_optimizer.py b/src/tequila/quantumchemistry/orbital_optimizer.py
index 06b36553..a7bca8d7 100644
--- a/src/tequila/quantumchemistry/orbital_optimizer.py
+++ b/src/tequila/quantumchemistry/orbital_optimizer.py
@@ -118,13 +118,15 @@ def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=N
         print(wrapper)
     if initial_guess is not None:
         if hasattr(initial_guess, "lower"):
-            if "random" in initial_guess.lower():
-                scale = 0.1
+            if "random" or "near_zero" in initial_guess.lower():
+                scale = 1.e-3
+                if "random" in initial_guess.lower():
+                    scale = 1.0
                 loc = 0.0
                 if "scale" in kwargs:
-                    scale = kwargs["scale"]
+                    scale = float(initial_guess.split("scale")[1].split("_")[0].split("=")[1])
                 if "loc" in kwargs:
-                    loc = kwargs["loc"]
+                    loc = float(initial_guess.split("loc")[1].split("_")[0].split("=")[1])
                 initial_guess = numpy.eye(no) + numpy.random.normal(scale=scale, loc=loc, size=no * no).reshape(no, no)
             else:
                 raise Exception("Unknown initial_guess={}".format(initial_guess.lower()))
diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index f2da7b28..9689b042 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -1939,7 +1939,7 @@ def _reset_rdm(rdm):
             self._rdm2 = _assemble_rdm2_spinful(evals_2) if get_rdm2 else self._rdm2
 
         if get_rdm2:
-            rdm2 = NBodyTensor(elems=self.rdm2, ordering="dirac")
+            rdm2 = NBodyTensor(elems=self.rdm2, ordering="dirac", verify=False)
             rdm2.reorder(to=ordering)
             rdm2 = rdm2.elems
             self._rdm2 = rdm2
diff --git a/src/tequila/version.py b/src/tequila/version.py
index f4531a81..a54a7564 100644
--- a/src/tequila/version.py
+++ b/src/tequila/version.py
@@ -1,2 +1,2 @@
-__version__ = "1.9.3"
+__version__ = "1.9.4"
 __author__ = "Tequila Developers "
diff --git a/tests/test_noise_opt.py b/tests/test_noise_opt.py
index 043159ed..a48b9c66 100644
--- a/tests/test_noise_opt.py
+++ b/tests/test_noise_opt.py
@@ -59,20 +59,6 @@ def test_bit_flip_scipy_hessian(p, method):
     result = tq.optimizer_scipy.minimize(objective=O, samples=1, backend=simulator, method=method, noise=NM, tol=1.e-4,
                                          silent=False)
 
-
-@pytest.mark.skipif(len(samplers) == 0, reason="Missing necessary backends")
-@pytest.mark.skipif(not tq.optimizers.has_phoenics, reason="Missing phoenics installation")
-@pytest.mark.parametrize("p", numpy.random.uniform(0.1, .4, 1))
-def test_bit_flip_phoenics(p):
-    simulator = numpy.random.choice(samplers)
-    qubit = 0
-    H = paulis.Qm(qubit)
-    U = gates.Rx(target=qubit, angle=tq.Variable('a'))
-    O = ExpectationValue(U=U, H=H)
-    NM = BitFlip(p, 1)
-    result = tq.optimizers.optimizer_phoenics.minimize(objective=O, maxiter=3, samples=1, backend=simulator, noise=NM)
-
-
 @pytest.mark.skipif(len(samplers) == 0, reason="Missing necessary backends")
 @pytest.mark.skipif(not tq.optimizers.has_gpyopt, reason="Missing gpyopt installation")
 @pytest.mark.parametrize("p", numpy.random.uniform(0.1, .4, 1))
diff --git a/tests/test_phoenics.py b/tests/test_phoenics.py
deleted file mode 100644
index 0cea65eb..00000000
--- a/tests/test_phoenics.py
+++ /dev/null
@@ -1,68 +0,0 @@
-import pytest, numpy
-import tequila as tq
-import multiprocessing as mp
-
-# Get QC backends for parametrized testing
-import select_backends
-simulators = select_backends.get()
-samplers = select_backends.get(sampler=True)
-
-has_phoenics = 'phoenics' in tq.INSTALLED_OPTIMIZERS
-
-@pytest.mark.dependencies
-def test_dependencies():
-    assert 'phoenics' in tq.INSTALLED_OPTIMIZERS
-
-
-@pytest.mark.skipif(condition=not has_phoenics, reason="you don't have phoenics")
-@pytest.mark.parametrize("simulator", simulators)
-def test_execution(simulator):
-    U = tq.gates.Rz(angle="a", target=0) \
-        + tq.gates.X(target=2) \
-        + tq.gates.Ry(angle="b", target=1, control=2) \
-        + tq.gates.Trotterized(angles=["c", "d"],
-                               generators=[-0.25 * tq.paulis.Z(1), tq.paulis.X(0) + tq.paulis.Y(1)], steps=2) \
-        + tq.gates.Trotterized(angles=[1.0, 2.0],
-                               generators=[-0.25 * tq.paulis.Z(1), tq.paulis.X(0) + tq.paulis.Y(1)], steps=2) \
-        + tq.gates.ExpPauli(angle="a", paulistring="X(0)Y(1)Z(2)")
-
-    H = 1.0 * tq.paulis.X(0) + 2.0 * tq.paulis.Y(1) + 3.0 * tq.paulis.Z(2)
-    O = tq.ExpectationValue(U=U, H=H)
-    result = tq.minimize(method="phoenics", objective=O, maxiter=1, backend=simulator)
-
-
-@pytest.mark.skipif(condition=not has_phoenics, reason="you don't have phoenics")
-@pytest.mark.parametrize("simulator", samplers)
-def test_execution_shot(simulator):
-    U = tq.gates.Rz(angle="a", target=0) \
-        + tq.gates.X(target=2) \
-        + tq.gates.Ry(angle="b", target=1, control=2) \
-        + tq.gates.Trotterized(angles=["c", "d"],
-                               generators=[-0.25 * tq.paulis.Z(1), tq.paulis.X(0) + tq.paulis.Y(1)], steps=2) \
-        + tq.gates.Trotterized(angles=[1.0, 2.0],
-                               generators=[-0.25 * tq.paulis.Z(1), tq.paulis.X(0) + tq.paulis.Y(1)], steps=2) \
-        + tq.gates.ExpPauli(angle="a", paulistring="X(0)Y(1)Z(2)")
-    H = 1.0 * tq.paulis.X(0) + 2.0 * tq.paulis.Y(1) + 3.0 * tq.paulis.Z(2)
-    O = tq.ExpectationValue(U=U, H=H)
-    mi = 2
-    result = tq.minimize(method="phoenics", objective=O, maxiter=mi, backend=simulator)
-
-
-@pytest.mark.skipif(condition=not has_phoenics, reason="you don't have phoenics")
-@pytest.mark.parametrize("simulator", simulators)
-def test_one_qubit_wfn(simulator):
-    U = tq.gates.Trotterized(angles=["a"], steps=1, generators=[tq.paulis.Y(0)])
-    H = tq.paulis.X(0)
-    O = tq.ExpectationValue(U=U, H=H)
-    result = tq.minimize(method="phoenics", objective=O, maxiter=8, backend=simulator)
-    assert (numpy.isclose(result.energy, -1.0, atol=1.e-2))
-
-
-@pytest.mark.skipif(condition=not has_phoenics, reason="you don't have phoenics")
-@pytest.mark.parametrize("simulator", samplers)
-def test_one_qubit_shot(simulator):
-    U = tq.gates.Trotterized(angles=["a"], steps=1, generators=[tq.paulis.Y(0)])
-    H = tq.paulis.X(0)
-    O = tq.ExpectationValue(U=U, H=H)
-    result = tq.minimize(method="phoenics", objective=O, maxiter=3, backend=simulator, samples=10000)
-    assert (numpy.isclose(result.energy, -1.0, atol=1.e-1))

From a3d532b7f041d4df932f6104f5e1622cbea6f05e Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 25 Mar 2024 10:56:11 +0100
Subject: [PATCH 2/6] Update README.md

---
 README.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/README.md b/README.md
index 3cec3281..f51e8b75 100644
--- a/README.md
+++ b/README.md
@@ -240,6 +240,10 @@ A.G. Cadavid, I. Montalban, A. Dalal, E. Solano, N.N. Hegade
 Efficient DCQO Algorithm within the Impulse Regime for Portfolio Optimization  
 [arxiv:2308.15475](https://arxiv.org/abs/2308.15475)  
 
+P.W.K. Jensen, E.R. Kjellgren, P. Reinholdt, K.M. Ziems, S. Coriani, J. Kongsted, S. Sauer  
+Quantum Equation of Motion with Orbital Optimization for Computing Molecular Properties in Near-Term Quantum Computing  
+[arxiv:2312.12386](https://arxiv.org/abs/2312.12386)  
+
 Let us know, if you want your research project and/or tutorial to be included in this list!
 
 # Dependencies

From a547f3c046c417f05687f30a7b04ec28362cd356 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 1 Apr 2024 23:45:59 +0200
Subject: [PATCH 3/6] Update README.md

---
 README.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/README.md b/README.md
index f51e8b75..9c350865 100644
--- a/README.md
+++ b/README.md
@@ -240,6 +240,10 @@ A.G. Cadavid, I. Montalban, A. Dalal, E. Solano, N.N. Hegade
 Efficient DCQO Algorithm within the Impulse Regime for Portfolio Optimization  
 [arxiv:2308.15475](https://arxiv.org/abs/2308.15475)  
 
+A. Anand, K. Brown  
+Hamiltonians, groups, graphs and ansätze  
+[arxiv:2312.17146](https://arxiv.org/abs/2312.17146)  
+
 P.W.K. Jensen, E.R. Kjellgren, P. Reinholdt, K.M. Ziems, S. Coriani, J. Kongsted, S. Sauer  
 Quantum Equation of Motion with Orbital Optimization for Computing Molecular Properties in Near-Term Quantum Computing  
 [arxiv:2312.12386](https://arxiv.org/abs/2312.12386)  

From 086ff99b9afd9f809f6d04788bb00f37b56a3d3b Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Fri, 5 Apr 2024 10:57:20 +0200
Subject: [PATCH 4/6] Update version.py

---
 src/tequila/version.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/tequila/version.py b/src/tequila/version.py
index cb9458d9..09ed92e6 100644
--- a/src/tequila/version.py
+++ b/src/tequila/version.py
@@ -1,2 +1,2 @@
-__version__ = "1.9.6-dev"
+__version__ = "1.9.5"
 __author__ = "Tequila Developers "

From 37b3e6a426011be3f735173bfbcb83d54af1f7e7 Mon Sep 17 00:00:00 2001
From: Deadly Artist <ocram.siebenstein@web.de>
Date: Mon, 22 Apr 2024 14:22:24 +0200
Subject: [PATCH 5/6] add givens methods

---
 src/tequila/quantumchemistry/qc_base.py | 194 +++++++++++++++++++++++-
 src/tequila/tools/random_generators.py  |  17 +++
 tests/test_givens.py                    |  60 ++++++++
 3 files changed, 268 insertions(+), 3 deletions(-)
 create mode 100644 tests/test_givens.py

diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index d8a6882b..94d15d80 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -17,7 +17,7 @@
 
 import typing, numpy, numbers
 from itertools import product
-
+import tequila.grouping.fermionic_functions as ff
 
 
 try:
@@ -32,8 +32,7 @@
     except Exception as E:
         raise Exception("{}\nIssue with Tequila Chemistry: Please update openfermion".format(str(E)))
 import warnings
-
-
+OPTIMIZED_ORDERING = "Optimized"
 class QuantumChemistryBase:
     """
     Base Class for tequila chemistry functionality
@@ -2127,6 +2126,58 @@ def perturbative_f12_correction(self, rdm1: numpy.ndarray = None, rdm2: numpy.nd
                                                    n_ri=n_ri, external_info=external_info, **kwargs)
         return correction.compute()
 
+    def n_rotation(self, i, phi):
+        '''
+        Creates a quantum circuit that applies a phase rotation based on phi to both components (up and down) of a given qubit.
+
+        Parameters:
+        - i (int): The index of the qubit to which the rotation will be applied.
+        - phi (float): The rotation angle. The actual rotation applied will be multiplied with -2 for both components.
+
+        Returns:
+        - QCircuit: A quantum circuit object containing the sequence of rotations applied to the up and down components of the specified qubit.
+        '''
+
+        # Generate number operators for the up and down components of the qubit.
+        n_up = self.make_number_op(2*i)
+        n_down = self.make_number_op(2*i+1)
+
+        # Start a new circuit and apply rotations to each component.
+        circuit = gates.GeneralizedRotation(generator = n_up, angle=-2*phi)
+        circuit += gates.GeneralizedRotation(generator = n_down, angle=-2*phi)
+        return circuit
+        
+    def get_givens_circuit(self, unitary, tol = 1e-12, ordering = OPTIMIZED_ORDERING):
+        '''
+        Constructs a quantum circuit from a given real unitary matrix using Givens rotations.
+        
+        This method decomposes a unitary matrix into a series of Givens and Rz (phase) rotations,
+        then constructs and returns a quantum circuit that implements this sequence of rotations.
+
+        Parameters:
+        - unitary (numpy.array): A real unitary matrix representing the transformation to implement.
+        - tol (float): A tolerance threshold below which matrix elements are considered zero.
+        - ordering (list of tuples or 'Optimized'): Custom ordering of indices for Givens rotations or 'Optimized' to generate them automatically.
+
+        Returns:
+        - QCircuit: A quantum circuit implementing the series of rotations decomposed from the unitary.
+        '''
+        # Decompose the unitary matrix into Givens and phase (Rz) rotations.
+        theta_list, phi_list = get_givens_decomposition(unitary, tol, ordering)
+
+        # Initialize an empty quantum circuit.
+        circuit = QCircuit()
+
+        # Add all Rz (phase) rotations to the circuit.
+        for phi in phi_list:
+            circuit += self.n_rotation(phi[1], phi[0])
+
+        # Add all Givens rotations to the circuit.
+        for theta in reversed(theta_list):
+            circuit += self.UR(theta[1], theta[2], theta[0]*2)
+
+        return circuit
+
 
     def print_basis_info(self):
         return self.integral_manager.print_basis_info()
@@ -2147,3 +2198,140 @@ def __str__(self) -> str:
         result += "\nmore information with: self.print_basis_info()\n"
 
         return result
+
+def givens_matrix(n, p, q, theta):
+    '''
+    Construct a complex Givens rotation matrix of dimension n by theta between rows/columns p and q.
+    '''
+    '''
+    Generates a Givens rotation matrix of size n x n to rotate by angle theta in the (p, q) plane. This matrix can be complex
+
+    Parameters:
+    - n (int): The size of the Givens rotation matrix.
+    - p (int): The first index for the rotation plane.
+    - q (int): The second index for the rotation plane.
+    - theta (float): The rotation angle.
+
+    Returns:
+    - numpy.array: The Givens rotation matrix.
+    '''
+    matrix = numpy.eye(n)  # Matrix to hold complex numbers
+    cos_theta = numpy.cos(theta)
+    sin_theta = numpy.sin(theta)
+
+    # Directly assign cosine and sine without complex phase adjustment
+    matrix[p, p] = cos_theta
+    matrix[q, q] = cos_theta
+    matrix[p, q] = sin_theta
+    matrix[q, p] = -sin_theta
+
+    return matrix
+
+def get_givens_decomposition(unitary, tol = 1e-12, ordering = OPTIMIZED_ORDERING, return_diagonal = False):
+    '''
+    Decomposes a real unitary matrix into Givens rotations (theta) and Rz rotations (phi).
+
+    Parameters:
+    - unitary (numpy.array): A real unitary matrix to decompose. It cannot be complex.
+    - tol (float): Tolerance for considering matrix elements as zero. Elements with absolute value less than tol are treated as zero.
+    - ordering (list of tuples or 'Optimized'): Custom ordering of indices for Givens rotations or 'Optimized' to generate them automatically.
+    - return_diagonal (bool): If True, the function also returns the diagonal matrix as part of the output.
+
+    Returns:
+    - list: A list of tuples, each representing a Givens rotation. Each tuple contains the rotation angle theta and indices (i,j) of the rotation.
+    - list: A list of tuples, each representing an Rz rotation. Each tuple contains the rotation angle phi and the index (i) of the rotation.
+    - numpy.array (optional): The diagonal matrix after applying all Givens rotations, returned if return_diagonal is True.
+    '''
+    U = unitary # no need to copy as we don't modify the original
+    U[abs(U) < tol] = 0 # Zeroing out the small elements as per the tolerance level.
+    n = U.shape[0]
+
+    # Determine optimized ordering if specified.
+    if ordering == OPTIMIZED_ORDERING:
+        ordering = ff.depth_eff_order_mf(n)
+
+    theta_list = []
+    phi_list = []
+
+    def calcTheta(U, c, r):
+        '''Calculate and apply the Givens rotation for a specific matrix element.'''
+        t = numpy.arctan2(-U[r,c], U[r-1,c])
+        theta_list.append((t, r, r-1))
+        g = givens_matrix(n,r,r-1,t)
+        U = numpy.dot(g, U)
+
+        return U
+
+    # Apply and store Givens rotations as per the given or computed ordering.
+    if ordering is None:
+        for c in range(n):
+            for r in range(n-1, c, -1):
+                U = calcTheta(U, c, r)
+    else:
+        for r, c in ordering:
+            U = calcTheta(U, c, r)
+    
+    # Calculating the Rz rotations based on the phases of the diagonal elements.
+    # For real elements this means a 180 degree shift, i.e. a sign change.
+    for i in range(n):
+        ph = numpy.angle(U[i,i])
+        phi_list.append((ph, i))
+
+    # Filtering out rotations without significance.
+    theta_list_new = []
+    for i, theta in enumerate(theta_list):
+        if abs(theta[0] % (2*numpy.pi)) > tol:
+            theta_list_new.append(theta)
+    
+    phi_list_new = []
+    for i, phi in enumerate(phi_list):
+        if abs(phi[0]) > tol:
+            phi_list_new.append(phi)
+        
+    if return_diagonal:
+        # Optionally return the resulting diagonal
+        return theta_list_new, phi_list_new, U
+    else:
+        return theta_list_new, phi_list_new
+    
+def reconstruct_matrix_from_givens(n, theta_list, phi_list, to_real_if_possible = True, tol = 1e-12):
+    '''
+    Reconstructs a matrix from given Givens rotations and Rz diagonal rotations.
+    This function is effectively an inverse of get_givens_decomposition, and therefore only works with data in the same format as its output.
+
+    Parameters:
+    - n (int): The size of the unitary matrix to be reconstructed.
+    - theta_list (list of tuples): Each tuple contains (angle, i, j) representing a Givens rotation of `angle` radians, applied to rows/columns `i` and `j`.
+    - phi_list (list of tuples): Each tuple contains (angle, i), representing an Rz rotation by `angle` radians applied to the `i`th diagonal element.
+    - to_real_if_possible (bool): If True, converts the matrix to real if its imaginary part is effectively zero.
+    - tol (float): The tolerance whether to swap a complex rotation for a sign change.
+
+    Returns:
+    - numpy.ndarray: The reconstructed complex or real matrix, depending on the `to_real_if_possible` flag and matrix composition.
+    '''
+    # Start with an identity matrix
+    reconstructed = numpy.eye(n, dtype=complex)
+    
+    # Apply Rz rotations for diagonal elements
+    for phi in phi_list:
+        angle, i = phi
+        # Directly apply a sign flip if the rotation angle is π
+        if numpy.isclose(angle, numpy.pi, atol=tol):
+            reconstructed[i, i] *= -1
+        else:
+            reconstructed[i, i] *= numpy.exp(1j * angle)
+        
+    # Apply Givens rotations in reverse order
+    for theta in reversed(theta_list):
+        angle, i, j = theta
+        g = givens_matrix(n, i, j, angle)
+        reconstructed = numpy.dot(g.conj().T, reconstructed) # Transpose of Givens matrix applied to the left
+
+    # Convert matrix to real if its imaginary part is negligible unless disabled via to_real_if_possible
+    if to_real_if_possible:
+        # Directly apply a sign flip if the rotation angle is π
+        if numpy.all(reconstructed.imag == 0):
+            # Convert to real by taking the real part
+            reconstructed = reconstructed.real
+
+    return reconstructed
diff --git a/src/tequila/tools/random_generators.py b/src/tequila/tools/random_generators.py
index 07f3a7cd..434022fc 100644
--- a/src/tequila/tools/random_generators.py
+++ b/src/tequila/tools/random_generators.py
@@ -2,6 +2,7 @@
 from tequila.circuit import gates
 from tequila.circuit.circuit import QCircuit
 from tequila.hamiltonian.qubit_hamiltonian import QubitHamiltonian
+from scipy.stats import unitary_group, ortho_group
 
 def make_random_circuit(n_qubits: int, rotation_gates: list=['rx', 'ry', 'rz'], n_rotations: int=None,
                         enable_controls: bool=None) -> QCircuit:
@@ -75,3 +76,19 @@ def make_random_hamiltonian(n_qubits: int , paulis: list=['X','Y','Z'], n_ps: in
     
     H = QubitHamiltonian(ham)
     return H
+
+def generate_random_unitary(size, complex = False):
+    '''
+    Generates a random unitary (or furthermore orthogonal if complex is False) matrix of a specified size.
+
+    Parameters:
+    - size (int): The size of the unitary matrix to be generated.
+    - complex (bool, optional): Whether the unitary should be complex.
+
+    Returns:
+    - numpy.ndarray: A randomly generated unitary matrix.
+    '''
+    if complex:
+        return unitary_group.rvs(size)
+    else:
+        return ortho_group.rvs(size)
\ No newline at end of file
diff --git a/tests/test_givens.py b/tests/test_givens.py
new file mode 100644
index 00000000..399c4d6d
--- /dev/null
+++ b/tests/test_givens.py
@@ -0,0 +1,60 @@
+import numpy
+import tequila as tq
+import tequila.quantumchemistry.qc_base as qc
+import tequila.tools.random_generators as rg
+import random
+
+transformations = ["JordanWigner", "ReorderedJordanWigner", "BravyiKitaev", "BravyiKitaevTree"]
+def test_givens_on_molecule():
+    # random size and transformation
+    size = random.randint(2, 10)
+    transformation = random.choice(transformations)
+
+    # dummy one-electron integrals
+    h = numpy.ones(shape=[size,size])
+    # dummy two-electron integrals
+    g = numpy.ones(shape=[size, size, size, size])
+
+    U = rg.generate_random_unitary(size)
+
+    # transformed integrals
+    th = (U.T.dot(h)).dot(U)
+    tg = numpy.einsum("ijkx, xl -> ijkl", g, U, optimize='greedy')
+    tg = numpy.einsum("ijxl, xk -> ijkl", tg, U, optimize='greedy')
+    tg = numpy.einsum("ixkl, xj -> ijkl", tg, U, optimize='greedy')
+    tg = numpy.einsum("xjkl, xi -> ijkl", tg, U, optimize='greedy')
+
+    # original molecule/H
+    mol = tq.Molecule(geometry="He 0.0 0.0 0.0", nuclear_repulsion=0.0, one_body_integrals=h, two_body_integrals=g, basis_set="dummy", transformation=transformation)
+    H = mol.make_hamiltonian()
+    # transformed molecule/H
+    tmol = tq.Molecule(geometry="He 0.0 0.0 0.0", nuclear_repulsion=0.0, one_body_integrals=th, two_body_integrals=tg,basis_set="dummy", transformation=transformation)
+    tH = tmol.make_hamiltonian()
+
+    # transformation in qubit space (this corresponds to the U above)
+    UR = mol.get_givens_circuit(U) # Works!
+
+    # test circuit
+    circuit = rg.make_random_circuit(size)
+
+    # create expectation values and see if they are the same
+    E1 = tq.ExpectationValue(U=circuit, H=tH)
+    E2 = tq.ExpectationValue(U=circuit + UR, H=H)
+
+    result1 = tq.simulate(E1)
+    result2 = tq.simulate(E2)
+    
+    assert numpy.isclose(result1, result2)
+
+def test_givens_decomposition():
+    # random unitary of random size
+    size = random.randint(2, 10)
+    unitary = rg.generate_random_unitary(size)
+
+    # decompose givens
+    theta_list, phi_list = qc.get_givens_decomposition(unitary)
+
+    # reconstruct original unitary from givens
+    reconstructed_matrix = qc.reconstruct_matrix_from_givens(unitary.shape[0], theta_list, phi_list)
+    
+    assert numpy.allclose(unitary, reconstructed_matrix)

From 878084f6d172bc30f35355849d365bf76483c415 Mon Sep 17 00:00:00 2001
From: Deadly Artist <ocram.siebenstein@web.de>
Date: Mon, 22 Apr 2024 14:37:21 +0200
Subject: [PATCH 6/6] move tests

---
 tests/test_chemistry.py | 54 +++++++++++++++++++++++++++++++++++++
 tests/test_givens.py    | 60 -----------------------------------------
 2 files changed, 54 insertions(+), 60 deletions(-)
 delete mode 100644 tests/test_givens.py

diff --git a/tests/test_chemistry.py b/tests/test_chemistry.py
index 86b2dc33..c4369442 100644
--- a/tests/test_chemistry.py
+++ b/tests/test_chemistry.py
@@ -8,6 +8,8 @@
 from tequila.objective import ExpectationValue
 from tequila.quantumchemistry.encodings import known_encodings
 from tequila.simulators.simulator_api import simulate
+import tequila.quantumchemistry.qc_base as qcb
+import tequila.tools.random_generators as rg
 
 HAS_PYSCF = "pyscf" in qc.INSTALLED_QCHEMISTRY_BACKENDS
 HAS_PSI4 = "psi4" in qc.INSTALLED_QCHEMISTRY_BACKENDS
@@ -725,3 +727,55 @@ def test_orbital_optimization_hcb(geometry):
     assert numpy.isclose(opt1.energy,opt2.energy,atol=1.e-5)
     assert time1 < time2
     assert (numpy.isclose(opt1.mo_coeff,opt2.mo_coeff, atol=1.e-5)).all()
+
+@pytest.mark.parametrize("transformation", ["JordanWigner", "ReorderedJordanWigner", "BravyiKitaev", "BravyiKitaevTree"])
+@pytest.mark.parametrize("size", [2, 8])
+def test_givens_on_molecule(size, transformation):
+    # dummy one-electron integrals
+    h = numpy.ones(shape=[size,size])
+    # dummy two-electron integrals
+    g = numpy.ones(shape=[size, size, size, size])
+
+    U = rg.generate_random_unitary(size)
+
+    # transformed integrals
+    th = (U.T.dot(h)).dot(U)
+    tg = numpy.einsum("ijkx, xl -> ijkl", g, U, optimize='greedy')
+    tg = numpy.einsum("ijxl, xk -> ijkl", tg, U, optimize='greedy')
+    tg = numpy.einsum("ixkl, xj -> ijkl", tg, U, optimize='greedy')
+    tg = numpy.einsum("xjkl, xi -> ijkl", tg, U, optimize='greedy')
+
+    # original molecule/H
+    mol = tq.Molecule(geometry="He 0.0 0.0 0.0", nuclear_repulsion=0.0, one_body_integrals=h, two_body_integrals=g, basis_set="dummy", transformation=transformation)
+    H = mol.make_hamiltonian()
+    # transformed molecule/H
+    tmol = tq.Molecule(geometry="He 0.0 0.0 0.0", nuclear_repulsion=0.0, one_body_integrals=th, two_body_integrals=tg,basis_set="dummy", transformation=transformation)
+    tH = tmol.make_hamiltonian()
+
+    # transformation in qubit space (this corresponds to the U above)
+    UR = mol.get_givens_circuit(U) # Works!
+
+    # test circuit
+    circuit = rg.make_random_circuit(size)
+
+    # create expectation values and see if they are the same
+    E1 = tq.ExpectationValue(U=circuit, H=tH)
+    E2 = tq.ExpectationValue(U=circuit + UR, H=H)
+
+    result1 = tq.simulate(E1)
+    result2 = tq.simulate(E2)
+    
+    assert numpy.isclose(result1, result2)
+
+@pytest.mark.parametrize("size", [2, 8])
+def test_givens_decomposition(size):
+    # generate random unitary
+    unitary = rg.generate_random_unitary(size)
+
+    # decompose givens
+    theta_list, phi_list = qcb.get_givens_decomposition(unitary)
+
+    # reconstruct original unitary from givens
+    reconstructed_matrix = qcb.reconstruct_matrix_from_givens(unitary.shape[0], theta_list, phi_list)
+    
+    assert numpy.allclose(unitary, reconstructed_matrix)
diff --git a/tests/test_givens.py b/tests/test_givens.py
deleted file mode 100644
index 399c4d6d..00000000
--- a/tests/test_givens.py
+++ /dev/null
@@ -1,60 +0,0 @@
-import numpy
-import tequila as tq
-import tequila.quantumchemistry.qc_base as qc
-import tequila.tools.random_generators as rg
-import random
-
-transformations = ["JordanWigner", "ReorderedJordanWigner", "BravyiKitaev", "BravyiKitaevTree"]
-def test_givens_on_molecule():
-    # random size and transformation
-    size = random.randint(2, 10)
-    transformation = random.choice(transformations)
-
-    # dummy one-electron integrals
-    h = numpy.ones(shape=[size,size])
-    # dummy two-electron integrals
-    g = numpy.ones(shape=[size, size, size, size])
-
-    U = rg.generate_random_unitary(size)
-
-    # transformed integrals
-    th = (U.T.dot(h)).dot(U)
-    tg = numpy.einsum("ijkx, xl -> ijkl", g, U, optimize='greedy')
-    tg = numpy.einsum("ijxl, xk -> ijkl", tg, U, optimize='greedy')
-    tg = numpy.einsum("ixkl, xj -> ijkl", tg, U, optimize='greedy')
-    tg = numpy.einsum("xjkl, xi -> ijkl", tg, U, optimize='greedy')
-
-    # original molecule/H
-    mol = tq.Molecule(geometry="He 0.0 0.0 0.0", nuclear_repulsion=0.0, one_body_integrals=h, two_body_integrals=g, basis_set="dummy", transformation=transformation)
-    H = mol.make_hamiltonian()
-    # transformed molecule/H
-    tmol = tq.Molecule(geometry="He 0.0 0.0 0.0", nuclear_repulsion=0.0, one_body_integrals=th, two_body_integrals=tg,basis_set="dummy", transformation=transformation)
-    tH = tmol.make_hamiltonian()
-
-    # transformation in qubit space (this corresponds to the U above)
-    UR = mol.get_givens_circuit(U) # Works!
-
-    # test circuit
-    circuit = rg.make_random_circuit(size)
-
-    # create expectation values and see if they are the same
-    E1 = tq.ExpectationValue(U=circuit, H=tH)
-    E2 = tq.ExpectationValue(U=circuit + UR, H=H)
-
-    result1 = tq.simulate(E1)
-    result2 = tq.simulate(E2)
-    
-    assert numpy.isclose(result1, result2)
-
-def test_givens_decomposition():
-    # random unitary of random size
-    size = random.randint(2, 10)
-    unitary = rg.generate_random_unitary(size)
-
-    # decompose givens
-    theta_list, phi_list = qc.get_givens_decomposition(unitary)
-
-    # reconstruct original unitary from givens
-    reconstructed_matrix = qc.reconstruct_matrix_from_givens(unitary.shape[0], theta_list, phi_list)
-    
-    assert numpy.allclose(unitary, reconstructed_matrix)