Repository: tequilahub/tequila (language: python)

Issue:
Add Qiskit GPU simulator
This pull request depends on https://github.com/tequilahub/tequila/pull/368.

I added support for running the Qiskit simulator on the GPU by creating a subclass that changes the device name. I used the same structure as for the Qulacs GPU simulator.

Using this requires the `qiskit-aer-gpu` package.

I also commented out the package check for `qulacs-gpu`, because that package is not maintained, instead you have to compile Qulacs with GPU support, see https://github.com/qulacs/qulacs/issues/623#issuecomment-2006123134.



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/ci_backends.yml b/.github/workflows/ci_backends.yml
index 9ae5d9ab..79717dbe 100644
--- a/.github/workflows/ci_backends.yml
+++ b/.github/workflows/ci_backends.yml
@@ -32,7 +32,7 @@ jobs:
             # 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>=0.30" "qulacs" "qibo==0.1.1" "myqlm" "cirq-google"
+            pip install "cirq" "qiskit" "qulacs" "qibo==0.1.1" "myqlm" "cirq-google"
         fi
         pip install -e .
     - name: Lint with flake8
diff --git a/README.md b/README.md
index 9ffaf33a..eecda10c 100644
--- a/README.md
+++ b/README.md
@@ -342,12 +342,7 @@ You can avoid it by downgrading cirq and openfermion
 ```bash
 pip install --upgrade "openfermion<=1.0.0"
 pip install --upgrade "cirq<=0.9.1"
-```  
-
-
-## Qiskit backend
-Qiskit version 0.25 is not yet supported.
-`pip install --upgrade qiskit<0.25` fixes potential issues. If not: Please let us know.
+```
 
 ## Circuit drawing
 Standard graphical circuit representation within a Jupyter environment is often done using `tq.draw`.
diff --git a/src/tequila/apps/unary_state_prep.py b/src/tequila/apps/unary_state_prep.py
index b59a0f98..55ad73e9 100644
--- a/src/tequila/apps/unary_state_prep.py
+++ b/src/tequila/apps/unary_state_prep.py
@@ -89,7 +89,6 @@ def __init__(self, target_space: typing.List[BitString], max_repeat: int = 100,
         simulator.convert_to_numpy = False
         variables = None # {k:k.name.evalf() for k in self._abstract_circuit.extract_variables()}
         wfn = simulator.simulate(initial_state=BitString.from_int(0, nbits=self.n_qubits), variables=variables)
-        wfn.n_qubits = self._n_qubits
         equations = []
         for k in target_space:
             equations.append(wfn[k] - abstract_coefficients[k])
@@ -174,7 +173,7 @@ def __call__(self, wfn: QubitWaveFunction) -> QCircuit:
         :return:
         """
         try:
-            assert (len(wfn) == len(self._target_space))
+            assert wfn.length() == len(self._target_space)
             for key in wfn.keys():
                 try:
                     assert (key in self._target_space)
diff --git a/src/tequila/hamiltonian/paulis.py b/src/tequila/hamiltonian/paulis.py
index 175c5850..d282aa4c 100644
--- a/src/tequila/hamiltonian/paulis.py
+++ b/src/tequila/hamiltonian/paulis.py
@@ -256,7 +256,7 @@ def Projector(wfn, threshold=0.0, n_qubits=None) -> QubitHamiltonian:
 
     """
 
-    wfn = QubitWaveFunction(state=wfn, n_qubits=n_qubits)
+    wfn = QubitWaveFunction.convert_from(n_qubits, wfn)
 
     H = QubitHamiltonian.zero()
     for k1, v1 in wfn.items():
@@ -304,8 +304,9 @@ def KetBra(ket: QubitWaveFunction, bra: QubitWaveFunction, hermitian: bool = Fal
 
     """
     H = QubitHamiltonian.zero()
-    ket = QubitWaveFunction(state=ket, n_qubits=n_qubits)
-    bra = QubitWaveFunction(state=bra, n_qubits=n_qubits)
+
+    ket = QubitWaveFunction.convert_from(n_qubits, ket)
+    bra = QubitWaveFunction.convert_from(n_qubits, bra)
 
     for k1, v1 in bra.items():
         for k2, v2 in ket.items():
diff --git a/src/tequila/quantumchemistry/encodings.py b/src/tequila/quantumchemistry/encodings.py
index d31ae761..c362bd98 100644
--- a/src/tequila/quantumchemistry/encodings.py
+++ b/src/tequila/quantumchemistry/encodings.py
@@ -128,9 +128,9 @@ def map_state(self, state: list, *args, **kwargs) -> list:
         fop = openfermion.FermionOperator(string, 1.0)
         op = self(fop)
         from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
-        wfn = QubitWaveFunction.from_int(0, n_qubits=n_qubits)
+        wfn = QubitWaveFunction.from_basis_state(n_qubits, 0)
         wfn = wfn.apply_qubitoperator(operator=op)
-        assert (len(wfn.keys()) == 1)
+        assert wfn.length() == 1
         key = list(wfn.keys())[0].array
         return key
 
diff --git a/src/tequila/quantumchemistry/orbital_optimizer.py b/src/tequila/quantumchemistry/orbital_optimizer.py
index 5972134a..1fbb843a 100644
--- a/src/tequila/quantumchemistry/orbital_optimizer.py
+++ b/src/tequila/quantumchemistry/orbital_optimizer.py
@@ -208,11 +208,16 @@ def kernel(self, h1, h2, *args, **kwargs):
             H = molecule.make_hardcore_boson_hamiltonian()
         else:
             H = molecule.make_hamiltonian()
+
+        rdm1 = None
+        rdm2 = None
         if self.vqe_solver is not None:
             vqe_solver_arguments = {}
             if self.vqe_solver_arguments is not None:
                 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)
         elif self.circuit is None:
             raise Exception("Orbital Optimizer: Either provide a callable vqe_solver or a circuit")
         else:
@@ -233,8 +238,9 @@ def kernel(self, h1, h2, *args, **kwargs):
             # static ansatz
             U = self.circuit
 
-        rdm1, rdm2 = molecule.compute_rdms(U=U, variables=result.variables, spin_free=True, get_rdm1=True, get_rdm2=True, use_hcb=restrict_to_hcb)
-        rdm2 = self.reorder(rdm2, 'dirac', 'mulliken')
+        if rdm1 is None or rdm2 is None:
+            rdm1, rdm2 = molecule.compute_rdms(U=U, variables=result.variables, spin_free=True, get_rdm1=True, get_rdm2=True, use_hcb=restrict_to_hcb)
+            rdm2 = self.reorder(rdm2, 'dirac', 'mulliken')
         if not self.silent:
             print("{:20} : {}".format("energy", result.energy))
             if len(self.history) > 0:
@@ -259,3 +265,4 @@ def __str__(self):
             else:
                 result += "{:30} : {}\n".format(k, v)
         return result
+
diff --git a/src/tequila/simulators/simulator_api.py b/src/tequila/simulators/simulator_api.py
index 9db3991a..8db872a9 100755
--- a/src/tequila/simulators/simulator_api.py
+++ b/src/tequila/simulators/simulator_api.py
@@ -9,9 +9,10 @@
 from tequila.utils.exceptions import TequilaException, TequilaWarning
 from tequila.simulators.simulator_base import BackendCircuit, BackendExpectationValue
 from tequila.circuit.noise import NoiseModel
+from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
 
-SUPPORTED_BACKENDS = ["qulacs_gpu", "qulacs",'qibo', "qiskit", "cirq", "pyquil", "symbolic", "qlm"]
-SUPPORTED_NOISE_BACKENDS = ["qiskit", 'cirq', 'pyquil'] # qulacs removed in v.1.9
+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
 BackendTypes = namedtuple('BackendTypes', 'CircType ExpValueType')
 INSTALLED_SIMULATORS = {}
 INSTALLED_SAMPLERS = {}
@@ -22,7 +23,6 @@
     from tequila.objective import Objective, Variable
     from tequila.circuit.gates import QCircuit
     import numbers.Real as RealNumber
-    from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
 
 """
 Check which simulators are installed
@@ -43,6 +43,19 @@
     HAS_QISKIT = False
     HAS_QISKIT_NOISE = False
 
+try:
+    pkg_resources.require("qiskit-aer-gpu")
+    from tequila.simulators.simulator_qiskit_gpu import BackendCircuitQiskitGpu, BackendExpectationValueQiskitGpu
+    HAS_QISKIT_GPU = True
+    INSTALLED_SIMULATORS["qiskit_gpu"] = BackendTypes(BackendCircuitQiskitGpu, BackendExpectationValueQiskitGpu)
+    INSTALLED_SAMPLERS["qiskit_gpu"] = BackendTypes(BackendCircuitQiskitGpu, BackendExpectationValueQiskitGpu)
+    from tequila.simulators.simulator_qiskit import HAS_NOISE as HAS_QISKIT_GPU_NOISE
+    if HAS_QISKIT_GPU_NOISE:
+        INSTALLED_NOISE_SAMPLERS["qiskit_gpu"] = BackendTypes(BackendCircuitQiskitGpu, BackendExpectationValueQiskitGpu)
+except (ImportError, DistributionNotFound):
+    HAS_QISKIT_GPU = False
+    HAS_QISKIT_GPU_NOISE = False
+
 HAS_QIBO = True
 try:
     from tequila.simulators.simulator_qibo import BackendCircuitQibo, BackendExpectationValueQibo
@@ -82,8 +95,8 @@
     HAS_QULACS = False
 
 try:
-    pkg_resources.require("qulacs-gpu")
-    import qulacs
+    # pkg_resources.require("qulacs-gpu")
+    from qulacs import QuantumStateGpu
     from tequila.simulators.simulator_qulacs_gpu import BackendCircuitQulacsGpu, BackendExpectationValueQulacsGpu
 
     HAS_QULACS_GPU = True
@@ -350,14 +363,15 @@ def compile_circuit(abstract_circuit: 'QCircuit',
     return CircType(abstract_circuit=abstract_circuit, variables=variables, noise=noise, device=device, *args, **kwargs)
 
 
-def simulate(objective: typing.Union['Objective', 'QCircuit','QTensor'],
+def simulate(objective: typing.Union['Objective', 'QCircuit', 'QTensor'],
              variables: Dict[Union[Variable, Hashable], RealNumber] = None,
              samples: int = None,
              backend: str = None,
              noise: NoiseModel = None,
              device: str = None,
+             initial_state: Union[int, QubitWaveFunction] = 0,
              *args,
-             **kwargs) -> Union[RealNumber, 'QubitWaveFunction']:
+             **kwargs) -> Union[RealNumber, QubitWaveFunction]:
     """Simulate a tequila objective or circuit
 
     Parameters
@@ -375,6 +389,8 @@ def simulate(objective: typing.Union['Objective', 'QCircuit','QTensor'],
         specify a noise model to apply to simulation/sampling
     device:
         a device upon which (or in emulation of which) to sample
+    initial_state: int or QubitWaveFunction:
+        the initial state of the circuit
     *args :
 
     **kwargs :
@@ -394,9 +410,9 @@ def simulate(objective: typing.Union['Objective', 'QCircuit','QTensor'],
                 objective.extract_variables()))
 
     compiled_objective = compile(objective=objective, samples=samples, variables=variables, backend=backend,
-                                 noise=noise,device=device, *args, **kwargs)
+                                 noise=noise, device=device, *args, **kwargs)
 
-    return compiled_objective(variables=variables, samples=samples, *args, **kwargs)
+    return compiled_objective(variables=variables, samples=samples, initial_state=initial_state, *args, **kwargs)
 
 
 def draw(objective, variables=None, backend: str = None, name=None, *args, **kwargs):
diff --git a/src/tequila/simulators/simulator_base.py b/src/tequila/simulators/simulator_base.py
index b765abb2..db67c5b3 100755
--- a/src/tequila/simulators/simulator_base.py
+++ b/src/tequila/simulators/simulator_base.py
@@ -7,6 +7,7 @@
 from tequila import BitString
 from tequila.objective.objective import Variable, format_variable_dictionary
 from tequila.circuit import compiler
+from typing import Union
 
 import numbers, typing, numpy, copy, warnings
 
@@ -107,6 +108,11 @@ class BackendCircuit():
         "cc_max": True
     }
 
+    # Can be overwritten by backends that allow basis state initialization when sampling
+    supports_sampling_initialization: bool = False
+    # Can be overwritten by backends that allow initializing arbitrary states
+    supports_generic_initialization: bool = False
+
     @property
     def n_qubits(self) -> numbers.Integral:
         return len(self.qubit_map)
@@ -328,7 +334,7 @@ def update_variables(self, variables):
         """
         self.circuit = self.create_circuit(abstract_circuit=self.abstract_circuit, variables=variables)
 
-    def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunction:
+    def simulate(self, variables, initial_state: Union[int, QubitWaveFunction] = 0, *args, **kwargs) -> QubitWaveFunction:
         """
         simulate the circuit via the backend.
 
@@ -348,13 +354,12 @@ def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunc
             the wavefunction of the system produced by the action of the circuit on the initial state.
 
         """
+        if isinstance(initial_state, QubitWaveFunction) and not self.supports_generic_initialization:
+            raise TequilaException("Backend does not support arbitrary initial states")
+
         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")
-            initial_state = list(initial_state.keys())[0].integer
 
         all_qubits = list(range(self.abstract_circuit.n_qubits))
         active_qubits = self.qubit_map.keys()
@@ -362,20 +367,30 @@ def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunc
         # Keymap is only necessary if not all qubits are active
         keymap_required = sorted(active_qubits) != all_qubits
 
+        # Combining keymap and general initial states is awkward, because it's not clear what should happen with
+        # different states on non-active qubits. For now, this is simply not allowed.
+        # A better solution might be to check if all components of the initial state differ only on the active qubits.
+        if keymap_required and isinstance(initial_state, QubitWaveFunction):
+            raise TequilaException("Can only set non-basis initial state if all qubits are used")
+
         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)
+        if not isinstance(initial_state, QubitWaveFunction):
+            mapped_initial_state = keymap.inverted(initial_state).integer if keymap_required else int(initial_state)
+        else:
+            mapped_initial_state = initial_state
+
         result = self.do_simulate(variables=variables, initial_state=mapped_initial_state, *args,
                                   **kwargs)
 
         if keymap_required:
-            result.apply_keymap(keymap=keymap, initial_state=initial_state)
+            result = QubitWaveFunction.from_wavefunction(result, keymap, n_qubits=len(all_qubits), initial_state=initial_state)
 
         return result
 
-    def sample(self, variables, samples, read_out_qubits=None, circuit=None, *args, **kwargs):
+    def sample(self, variables, samples, read_out_qubits=None, circuit=None, initial_state=0, *args, **kwargs):
         """
         Sample the circuit. If circuit natively equips paulistrings, sample therefrom.
         Parameters
@@ -395,6 +410,12 @@ def sample(self, variables, samples, read_out_qubits=None, circuit=None, *args,
             The result of sampling, a recreated QubitWaveFunction in the sampled basis.
 
         """
+        if initial_state != 0 and not self.supports_sampling_initialization:
+            raise TequilaException("Backend does not support initial states for sampling")
+
+        if isinstance(initial_state, QubitWaveFunction) and not self.supports_generic_initialization:
+            raise TequilaException("Backend does not support arbitrary initial states")
+
         self.update_variables(variables)
         if read_out_qubits is None:
             read_out_qubits = self.abstract_qubits
@@ -406,7 +427,9 @@ def sample(self, variables, samples, read_out_qubits=None, circuit=None, *args,
             circuit = self.add_measurement(circuit=self.circuit, target_qubits=read_out_qubits)
         else:
             circuit = self.add_measurement(circuit=circuit, target_qubits=read_out_qubits)
-        return self.do_sample(samples=samples, circuit=circuit, read_out_qubits=read_out_qubits, *args, **kwargs)
+
+        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):
         """
@@ -511,7 +534,7 @@ def sample_paulistring(self, samples: int, paulistring, variables, *args,
         E = E / samples * paulistring.coeff
         return E
 
-    def do_sample(self, samples, circuit, noise, abstract_qubits=None, *args, **kwargs) -> QubitWaveFunction:
+    def do_sample(self, samples, circuit, noise, abstract_qubits=None, initial_state=0, *args, **kwargs) -> QubitWaveFunction:
         """
         helper function for sampling. MUST be overwritten by inheritors.
 
@@ -777,7 +800,7 @@ def __call__(self, variables, samples: int = None, *args, **kwargs):
                 raise TequilaException(
                     "BackendExpectationValue received not all variables. Circuit depends on variables {}, you gave {}".format(
                         self._variables, variables))
-        
+
         if samples is None:
             data = self.simulate(variables=variables, *args, **kwargs)
         else:
@@ -856,7 +879,7 @@ def sample(self, variables, samples, *args, **kwargs) -> numpy.array:
             samples = max(1, int(self.abstract_expectationvalue.samples * total_samples))
             suggested = samples
             # samples are not necessarily set (either the user has to set it or some functions like optimize_measurements)
- 
+
         if suggested is not None and suggested != samples:
             warnings.warn("simulating with samples={}, but expectationvalue carries suggested samples={}\nTry calling with samples='auto-total#ofsamples'".format(samples, suggested), TequilaWarning)
 
diff --git a/src/tequila/simulators/simulator_cirq.py b/src/tequila/simulators/simulator_cirq.py
index 5115f90f..7a9505bb 100755
--- a/src/tequila/simulators/simulator_cirq.py
+++ b/src/tequila/simulators/simulator_cirq.py
@@ -173,7 +173,7 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
         simulator = cirq.Simulator()
         backend_result = simulator.simulate(program=self.circuit, param_resolver=self.resolver,
                                             initial_state=initial_state)
-        return QubitWaveFunction.from_array(arr=backend_result.final_state_vector, numbering=self.numbering)
+        return QubitWaveFunction.from_array(array=backend_result.final_state_vector, numbering=self.numbering)
 
     def convert_measurements(self, backend_result: cirq.Result) -> QubitWaveFunction:
         """
@@ -186,18 +186,18 @@ def convert_measurements(self, backend_result: cirq.Result) -> QubitWaveFunction
         Returns
         -------
         QubitWaveFunction:
-            the result of sampling, as a tequila QubitWavefunction.
+            the result of sampling, as a tequila QubitWaveFunction.
 
         """
         assert (len(backend_result.measurements) == 1)
         for key, value in backend_result.measurements.items():
-            counter = QubitWaveFunction()
+            counter = QubitWaveFunction(self.n_qubits, self.numbering)
             for sample in value:
                 binary = BitString.from_array(array=sample.astype(int))
-                if binary in counter._state:
-                    counter._state[binary] += 1
+                if binary in counter.keys():
+                    counter[binary] += 1
                 else:
-                    counter._state[binary] = 1
+                    counter[binary] = 1
             return counter
 
     def do_sample(self, samples, circuit, *args, **kwargs) -> QubitWaveFunction:
diff --git a/src/tequila/simulators/simulator_pyquil.py b/src/tequila/simulators/simulator_pyquil.py
index 15c36bca..8261afc1 100755
--- a/src/tequila/simulators/simulator_pyquil.py
+++ b/src/tequila/simulators/simulator_pyquil.py
@@ -439,7 +439,7 @@ def do_simulate(self, variables, initial_state, *args, **kwargs):
             if val > 0:
                 iprep += pyquil.gates.X(i)
         backend_result = simulator.wavefunction(iprep + self.circuit, memory_map=self.resolver)
-        return QubitWaveFunction.from_array(arr=backend_result.amplitudes, numbering=self.numbering)
+        return QubitWaveFunction.from_array(array=backend_result.amplitudes, numbering=self.numbering)
 
     def do_sample(self, samples, circuit, *args, **kwargs) -> QubitWaveFunction:
         """
@@ -495,7 +495,7 @@ def string_to_array(s):
                     listing.append(int(letter))
             return listing
 
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self.n_qubits, self.numbering)
         bit_dict = {}
         for b in backend_result:
             try:
@@ -505,7 +505,7 @@ def string_to_array(s):
 
         for k, v in bit_dict.items():
             arr = string_to_array(k)
-            result._state[BitString.from_array(arr)] = v
+            result[BitString.from_array(arr)] = v
         return result
 
     def no_translation(self, abstract_circuit):
diff --git a/src/tequila/simulators/simulator_qibo.py b/src/tequila/simulators/simulator_qibo.py
index 027b22cb..4da85352 100755
--- a/src/tequila/simulators/simulator_qibo.py
+++ b/src/tequila/simulators/simulator_qibo.py
@@ -356,20 +356,20 @@ def do_simulate(self, variables, initial_state=None, *args, **kwargs):
         n_qubits = max(self.highest_qubit + 1, self.n_qubits, self.abstract_circuit.max_qubit() + 1)
         if initial_state is not None:
             if isinstance(initial_state, (int, np.int64)):
-                wave = QubitWaveFunction.from_int(i=initial_state, n_qubits=n_qubits)
+                wave = QubitWaveFunction.from_basis_state(n_qubits, initial_state, self.numbering)
             elif isinstance(initial_state, str):
-                wave = QubitWaveFunction.from_string(string=initial_state).to_array()
+                wave = QubitWaveFunction.from_string(initial_state, self.numbering).to_array()
             elif isinstance(initial_state, QubitWaveFunction):
                 wave = initial_state
             elif isinstance(initial_state,np.ndarray):
-                wave = QubitWaveFunction.from_array(initial_state)
+                wave = QubitWaveFunction.from_array(initial_state, self.numbering)
             else:
                 raise TequilaQiboException('could not understand initial state of type {}'.format(type(initial_state)))
             state = wave.to_array()
             result = self.circuit(state)
         else:
             result = self.circuit()
-        back= QubitWaveFunction.from_array(arr=result.numpy())
+        back= QubitWaveFunction.from_array(result.numpy(), self.numbering)
         return back
 
     def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunction:
@@ -398,7 +398,7 @@ def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunc
         if isinstance(initial_state, BitString):
             initial_state = initial_state.integer
         if isinstance(initial_state, QubitWaveFunction):
-            if len(initial_state.keys()) != 1:
+            if initial_state.length() != 1:
                 return self.do_simulate(variables=variables,initial_state=initial_state, *args, **kwargs)
             initial_state = list(initial_state.keys())[0].integer
         if isinstance(initial_state,np.ndarray):
@@ -426,19 +426,18 @@ def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveF
             results transformed to tequila native QubitWaveFunction
         """
 
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self.n_qubits, self.numbering)
         # todo there are faster ways
 
         for k, v in backend_result.frequencies(binary=True).items():
             converted_key = BitString.from_bitstring(other=BitString.from_binary(binary=k))
-            result._state[converted_key] = v
-
+            result[converted_key] = v
 
         if target_qubits is not None:
             mapped_target = [self.qubit_map[q].number for q in target_qubits]
             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)
+            result = QubitWaveFunction.from_wavefunction(result, keymap, len(mapped_target))
 
         return result
 
@@ -521,21 +520,20 @@ def do_sample(self, samples, circuit, noise_model=None, initial_state=None, *arg
         n_qubits = max(self.highest_qubit + 1, self.n_qubits, self.abstract_circuit.max_qubit() + 1)
         if initial_state is not None:
             if isinstance(initial_state, int):
-                wave=QubitWaveFunction.from_int(i=initial_state, n_qubits=n_qubits)
+                wave = QubitWaveFunction.from_basis_state(n_qubits, initial_state, self.numbering)
             elif isinstance(initial_state, str):
-                wave = QubitWaveFunction.from_string(string=initial_state).to_array()
+                wave = QubitWaveFunction.from_string(initial_state, self.numbering).to_array()
             elif isinstance(initial_state, QubitWaveFunction):
                 wave = initial_state
-            elif isinstance(initial_state,np.ndarray):
-                wave = QubitWaveFunction.from_array(arr=initial_state, n_qubits=n_qubits)  # silly but necessary
+            elif isinstance(initial_state, np.ndarray):
+                wave = QubitWaveFunction.from_array(initial_state, self.numbering)  # silly but necessary
             else:
                 raise TequilaQiboException('received an unusable initial state of type {}'.format(type(initial_state)))
-            state=wave.to_array()
-            result = circuit(state,nshots=samples)
+            state = wave.to_array()
+            result = circuit(state, nshots=samples)
         else:
             result = circuit(nshots=samples)
 
-
         back = self.convert_measurements(backend_result=result)
         return back
 
diff --git a/src/tequila/simulators/simulator_qiskit.py b/src/tequila/simulators/simulator_qiskit.py
index 9d2130aa..ee732518 100755
--- a/src/tequila/simulators/simulator_qiskit.py
+++ b/src/tequila/simulators/simulator_qiskit.py
@@ -1,22 +1,28 @@
 from tequila.simulators.simulator_base import BackendCircuit, QCircuit, BackendExpectationValue
+from tequila.utils.bitstrings import reverse_int_bits
 from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
 from tequila import TequilaException, TequilaWarning
 from tequila import BitString, BitNumbering, BitStringLSB
 from tequila.utils.keymap import KeyMapRegisterToSubregister
 from tequila.utils import to_float
-import qiskit, numpy, warnings
+from typing import Union
+import warnings
+import numpy as np
+import qiskit, qiskit_aer, qiskit.providers.fake_provider
+from qiskit import QuantumCircuit
 
-HAS_NOISE=True
+HAS_NOISE = True
 try:
-    from qiskit_aer import noise as qiskitnoise
+    from qiskit_aer import noise as qiskitnoise, AerSimulator
 except:
     HAS_NOISE = False
 
-HAS_IBMQ=True
+HAS_IBMQ = True
 try:
-    from qiskit.providers.ibmq import IBMQBackend
+    from qiskit_ibm_runtime import IBMBackend
 except:
-    HAS_IBMQ=False
+    HAS_IBMQ = False
+
 
 def get_bit_flip(p):
     """
@@ -70,16 +76,24 @@ def get_phase_flip(p):
     'multicontrol': 3
 }
 
-full_basis = ['x', 'y', 'z', 'id', 'u1', 'u2', 'u3', 'h','unitary','sx',
+full_basis = ['x', 'y', 'z', 'id', 'u1', 'u2', 'u3', 'h', 'unitary', 'sx',
               'cx', 'cy', 'cz', 'cu3', 'ccx']
 
+
 def qiskit_device_dict():
-    devices = {}
-    devices.update({str(x).lower():x for x in qiskit.Aer.backends()})
-    devices.update({str(x).lower():x for x in qiskit.test.mock.FakeProvider().backends()})
+    # As of Quiskit Aer 0.15.0, backends() also returns legacy backends which are not AerSimulators that will be
+    # deprecated in the future
+    # TODO: Instead of building a dict with all backends just to search it, search directly with `get_backend`
+    # https://qiskit.github.io/qiskit-aer/stubs/qiskit_aer.AerProvider.html#qiskit_aer.AerProvider.get_backend
+    devices = {backend.name.lower(): backend for backend in qiskit_aer.AerProvider().backends()
+               if isinstance(backend, AerSimulator)}
+
+    # FakeProvider has been removed, see https://github.com/Qiskit/qiskit/issues/10954
+    # devices.update({str(x).lower(): x for x in qiskit.test.mock.FakeProvider().backends()})
 
     return devices
 
+
 class TequilaQiskitException(TequilaException):
     def __str__(self):
         return "Error in qiskit backend:" + self.message
@@ -125,6 +139,8 @@ class BackendCircuitQiskit(BackendCircuit):
         transform a tequila NoiseModel into a qiskit noise model.
 
     """
+    STATEVECTOR_DEVICE_NAME = "aer_simulator_statevector"
+
     compiler_arguments = {
         "trotterized": True,
         "swap": False,
@@ -145,7 +161,10 @@ class BackendCircuitQiskit(BackendCircuit):
     }
 
     numbering = BitNumbering.LSB
-            
+
+    supports_sampling_initialization = True
+    supports_generic_initialization = True
+
     def __init__(self, abstract_circuit: QCircuit, variables, qubit_map=None, noise=None,
                  device=None, *args, **kwargs):
         """
@@ -188,8 +207,9 @@ def __init__(self, abstract_circuit: QCircuit, variables, qubit_map=None, noise=
         if qubit_map is None:
             qubit_map = {q: i for i, q in enumerate(abstract_circuit.qubits)}
         else:
-            warnings.warn("reveived custom qubit_map = {}\n"
-                          "This is not fully integrated with qiskit and might result in unexpected behaviour".format(qubit_map), TequilaWarning)
+            warnings.warn(f"reveived custom qubit_map = {qubit_map}\n"
+                          "This is not fully integrated with qiskit and might result in unexpected behaviour",
+                          TequilaWarning)
 
         n_qubits = max(qubit_map.values()) + 1
 
@@ -201,7 +221,7 @@ def __init__(self, abstract_circuit: QCircuit, variables, qubit_map=None, noise=
 
         self.classical_map = self.make_classical_map(qubit_map=self.qubit_map)
 
-        if noise != None:
+        if noise is not None:
             self.noise_lookup = {
                 'phase damp': qiskitnoise.phase_damping_error,
                 'amplitude damp': qiskitnoise.amplitude_damping_error,
@@ -211,12 +231,12 @@ def __init__(self, abstract_circuit: QCircuit, variables, qubit_map=None, noise=
                 'depolarizing': qiskitnoise.depolarizing_error
             }
 
-            if isinstance(noise, str): #string noise means "use the same noise as the device I tell you to get."
+            if isinstance(noise, str):  # string noise means "use the same noise as the device I tell you to get."
                 try:
                     self.check_device(noise)
                     self.noise_model = qiskitnoise.NoiseModel.from_backend(noise)
                 except TequilaQiskitException:
-                    raise TequilaException("noise init from string requires that noise names a device. Got {}".format(noise))
+                    raise TequilaException(f"noise init from string requires that noise names a device. Got {noise}")
 
             else:
                 self.noise_model = self.noise_model_converter(noise)
@@ -235,7 +255,7 @@ def make_qubit_map(self, qubits: dict = None):
         qubit_map = super().make_qubit_map(qubits=qubits)
         mapped_qubits = [q.number for q in qubit_map.values()]
         for k, v in qubit_map.items():
-            qubit_map[k].instance = self.q [v.number]
+            qubit_map[k].instance = self.q[v.number]
 
         return qubit_map
 
@@ -247,6 +267,20 @@ def make_classical_map(self, qubit_map: dict):
 
         return classical_map
 
+    def add_state_init(self, circuit: QuantumCircuit, initial_state: Union[int, QubitWaveFunction]) -> QuantumCircuit:
+        if initial_state == 0:
+            return circuit
+
+        if isinstance(initial_state, QubitWaveFunction):
+            statevector = initial_state.to_array(self.numbering)
+        else:
+            statevector = np.zeros(2 ** self.n_qubits)
+            statevector[reverse_int_bits(initial_state, self.n_qubits)] = 1.0
+
+        init_circuit = qiskit.QuantumCircuit(self.q, self.c)
+        init_circuit.set_statevector(statevector)
+        return init_circuit.compose(circuit)
+
     def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunction:
         """
         Helper function for performing simulation.
@@ -266,10 +300,12 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
         """
         if self.noise_model is None:
             if self.device is None:
-                qiskit_backend = self.retrieve_device('statevector_simulator')
+                qiskit_backend = self.retrieve_device(self.STATEVECTOR_DEVICE_NAME)
             else:
                 if 'statevector' not in str(self.device):
-                    raise TequilaException('For simulation, only state vector simulators are supported; recieved device={}, you might have forgoten to set the samples keyword - e.g. (device={}, samples=1000). If not set, tequila assumes that full wavefunction simualtion is demanded which is not compatible with qiskit devices or fake devices except for device=statevector'.format(self.device, self.device))
+                    raise TequilaException(
+                        'For simulation, only state vector simulators are supported; recieved device={}, you might have forgoten to set the samples keyword - e.g. (device={}, samples=1000). If not set, tequila assumes that full wavefunction simualtion is demanded which is not compatible with qiskit devices or fake devices except for device=statevector'.format(
+                            self.device, self.device))
                 else:
                     qiskit_backend = self.retrieve_device(self.device)
         else:
@@ -279,22 +315,18 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
         if "optimization_level" in kwargs:
             optimization_level = kwargs['optimization_level']
 
-        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}
+        circuit = self.circuit.assign_parameters(self.resolver)
 
-        circuit = self.circuit.bind_parameters(self.resolver)
+        circuit = self.add_state_init(circuit, initial_state)
 
-        qiskit_job = qiskit_backend.run(circuit,optimization_level=optimization_level,**opts)
+        circuit.save_statevector()
 
-        backend_result = qiskit_job.result()
-        return QubitWaveFunction.from_array(arr=backend_result.get_statevector(circuit), numbering=self.numbering)
+        backend_result = qiskit_backend.run(circuit, optimization_level=optimization_level).result()
 
-    def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubits, *args, **kwargs) -> QubitWaveFunction:
+        return QubitWaveFunction.from_array(array=backend_result.get_statevector(circuit).data, numbering=self.numbering)
+
+    def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubits, initial_state=0, *args,
+                  **kwargs) -> QubitWaveFunction:
         """
         Helper function for performing sampling.
         Parameters
@@ -303,6 +335,8 @@ def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubit
             the circuit from which to sample.
         samples:
             the number of samples to take.
+        initial_state:
+            initial state of the circuit
         args
         kwargs
 
@@ -319,17 +353,18 @@ def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubit
         else:
             qiskit_backend = self.retrieve_device(self.device)
 
-        if isinstance(qiskit_backend,IBMQBackend):
+        if isinstance(qiskit_backend, IBMBackend):
             if self.noise_model is not None:
-                raise TequilaException('Cannot combine backend {} with custom noise models.'.format(str(qiskit_backend)))
-            circuit = circuit.bind_parameters(self.resolver)  # this is necessary in spite of qiskit "fixing" it
+                raise TequilaException(
+                    'Cannot combine backend {} with custom noise models.'.format(str(qiskit_backend)))
+            circuit = circuit.assign_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),
+            return self.convert_measurements(qiskit_backend.run(circuit, shots=samples,
+                                                                optimization_level=optimization_level),
                                              target_qubits=read_out_qubits)
         else:
-            if isinstance(qiskit_backend, qiskit.test.mock.FakeBackend):
-                circuit = circuit.bind_parameters(self.resolver)  # this is necessary in spite of qiskit "fixing" it
+            if isinstance(qiskit_backend, qiskit.providers.fake_provider.FakeBackend):
+                circuit = circuit.assign_parameters(self.resolver)  # this is necessary in spite of qiskit "fixing" it
                 coupling_map = qiskit_backend.configuration().coupling_map
                 from_back = qiskitnoise.NoiseModel.from_backend(qiskit_backend)
                 if self.noise_model is not None:
@@ -343,23 +378,23 @@ 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)
+                return self.convert_measurements(job, target_qubits=read_out_qubits)
             else:
                 if self.noise_model is not None:
                     qiskit_backend.set_options(noise_model=self.noise_model)  # fits better with our methodology.
                     use_basis = full_basis
                 else:
                     use_basis = qiskit_backend.configuration().basis_gates
-                circuit = circuit.bind_parameters(self.resolver)  # this is necessary -- see qiskit-aer issue 1346
+                circuit = circuit.assign_parameters(self.resolver)  # this is necessary -- see qiskit-aer issue 1346
+                circuit = self.add_state_init(circuit, initial_state)
                 circuit = qiskit.transpile(circuit, backend=qiskit_backend,
                                            basis_gates=use_basis,
                                            optimization_level=optimization_level
                                            )
 
                 job = qiskit_backend.run(circuit, shots=samples)
-                return self.convert_measurements(job,
-                                                 target_qubits=read_out_qubits)
+                return self.convert_measurements(job, target_qubits=read_out_qubits)
 
     def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveFunction:
         """
@@ -374,16 +409,17 @@ def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveF
             measurements converted into wave function form.
         """
         qiskit_counts = backend_result.result().get_counts()
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self.n_qubits, self.numbering)
         # todo there are faster ways
         for k, v in qiskit_counts.items():
-            converted_key = BitString.from_bitstring(other=BitStringLSB.from_binary(binary=k))
-            result._state[converted_key] = v
+            # Qiskit uses LSB bitstrings, but from_binary expects MSB
+            converted_key = BitString.from_binary(k[::-1])
+            result[converted_key] = v
         if target_qubits is not None:
             mapped_target = [self.qubit_map[q].number for q in target_qubits]
             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)
+            result = QubitWaveFunction.from_wavefunction(result, keymap, n_qubits=len(target_qubits))
 
         return result
 
@@ -436,7 +472,8 @@ def add_parametrized_gate(self, gate, circuit, *args, **kwargs):
             par = float(gate.parameter)
         if gate.is_controlled():
             if len(gate.control) > 2:
-                raise TequilaQiskitException("multi-controls beyond 2 not yet supported for the qiskit backend. Gate was:\n{}".format(gate) )
+                raise TequilaQiskitException(
+                    "multi-controls beyond 2 not yet supported for the qiskit backend. Gate was:\n{}".format(gate))
             ops[1](circuit)(par, self.qubit(gate.control[0]), self.qubit(gate.target[0]))
         else:
             ops[0](circuit)(par, self.qubit(gate.target[0]))
@@ -590,7 +627,7 @@ def check_device(self, device):
         if device is None:
             return
 
-        elif isinstance(device,qiskit.providers.Backend):
+        elif isinstance(device, qiskit.providers.Backend):
             return
 
         elif isinstance(device, dict):
diff --git a/src/tequila/simulators/simulator_qiskit_gpu.py b/src/tequila/simulators/simulator_qiskit_gpu.py
new file mode 100644
index 00000000..52e461f2
--- /dev/null
+++ b/src/tequila/simulators/simulator_qiskit_gpu.py
@@ -0,0 +1,9 @@
+from tequila.simulators.simulator_qiskit import BackendCircuitQiskit, BackendExpectationValueQiskit
+
+
+class BackendCircuitQiskitGpu(BackendCircuitQiskit):
+    STATEVECTOR_DEVICE_NAME = "aer_simulator_statevector_gpu"
+
+
+class BackendExpectationValueQiskitGpu(BackendExpectationValueQiskit):
+    BackendCircuitType = BackendCircuitQiskitGpu
diff --git a/src/tequila/simulators/simulator_qlm.py b/src/tequila/simulators/simulator_qlm.py
index 1f23f174..2604cf4d 100644
--- a/src/tequila/simulators/simulator_qlm.py
+++ b/src/tequila/simulators/simulator_qlm.py
@@ -395,10 +395,10 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
         if MY_QLM:
             result = PyLinalg().submit(job)
             statevector = get_statevector(result)
-            return QubitWaveFunction.from_array(arr=statevector, numbering=self.numbering)
+            return QubitWaveFunction.from_array(array=statevector, numbering=self.numbering)
 
         result = LinAlg().submit(job)
-        return QubitWaveFunction.from_array(arr=result.statevector, numbering=self.numbering)
+        return QubitWaveFunction.from_array(array=result.statevector, numbering=self.numbering)
 
     def update_variables(self, variables):
         """
@@ -431,7 +431,7 @@ def convert_measurements(self, backend_result) -> QubitWaveFunction:
         QubitWaveFunction:
             measurements converted into wave function form.
         """
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self.n_qubits, self.numbering)
         shots = int(backend_result.meta_data["nbshots"])
         nbits = backend_result[0].qregs[0].length
         for sample in backend_result:
diff --git a/src/tequila/simulators/simulator_qulacs.py b/src/tequila/simulators/simulator_qulacs.py
index 29cfca95..71da737c 100755
--- a/src/tequila/simulators/simulator_qulacs.py
+++ b/src/tequila/simulators/simulator_qulacs.py
@@ -1,9 +1,11 @@
+from typing import Union
+
 import qulacs
 import numbers, numpy
 import warnings
 
 from tequila import TequilaException, TequilaWarning
-from tequila.utils.bitstrings import BitNumbering, BitString, BitStringLSB
+from tequila.utils.bitstrings import BitNumbering, BitString, BitStringLSB, reverse_int_bits
 from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
 from tequila.simulators.simulator_base import BackendCircuit, BackendExpectationValue, QCircuit, change_basis
 from tequila.utils.keymap import KeyMapRegisterToSubregister
@@ -63,6 +65,11 @@ class BackendCircuitQulacs(BackendCircuit):
 
     numbering = BitNumbering.LSB
 
+    quantum_state_class = qulacs.QuantumState
+
+    supports_sampling_initialization = True
+    supports_generic_initialization = True
+
     def __init__(self, abstract_circuit, noise=None, *args, **kwargs):
         """
 
@@ -110,10 +117,17 @@ def __init__(self, abstract_circuit, noise=None, *args, **kwargs):
 
             self.circuit=self.add_noise_to_circuit(noise)
 
-    def initialize_state(self, n_qubits:int=None) -> qulacs.QuantumState:
+    def initialize_state(self, n_qubits: int = None, initial_state: Union[int, QubitWaveFunction] = None) -> qulacs.QuantumState:
         if n_qubits is None:
             n_qubits = self.n_qubits
-        return qulacs.QuantumState(n_qubits)
+
+        state = self.quantum_state_class(n_qubits)
+        if isinstance(initial_state, int):
+            state.set_computational_basis(reverse_int_bits(initial_state, self.n_qubits))
+        elif isinstance(initial_state, QubitWaveFunction):
+            state.load(initial_state.to_array(self.numbering))
+
+        return state
 
     def update_variables(self, variables):
         """
@@ -130,7 +144,7 @@ def update_variables(self, variables):
         for k, angle in enumerate(self.variables):
             self.circuit.set_parameter(k, angle(variables))
 
-    def do_simulate(self, variables, initial_state, *args, **kwargs):
+    def do_simulate(self, variables, initial_state: Union[int, QubitWaveFunction], *args, **kwargs):
         """
         Helper function to perform simulation.
 
@@ -148,12 +162,10 @@ def do_simulate(self, variables, initial_state, *args, **kwargs):
         QubitWaveFunction:
             QubitWaveFunction representing result of the simulation.
         """
-        state = self.initialize_state(self.n_qubits)
-        lsb = BitStringLSB.from_int(initial_state, nbits=self.n_qubits)
-        state.set_computational_basis(BitString.from_binary(lsb.binary).integer)
+        state = self.initialize_state(self.n_qubits, initial_state)
         self.circuit.update_quantum_state(state)
 
-        wfn = QubitWaveFunction.from_array(arr=state.get_vector(), numbering=self.numbering)
+        wfn = QubitWaveFunction.from_array(array=state.get_vector(), numbering=self.numbering)
         return wfn
 
     def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveFunction:
@@ -170,22 +182,17 @@ def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveF
             results transformed to tequila native QubitWaveFunction
         """
 
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self.n_qubits, self.numbering)
         # todo there are faster ways
 
-
         for k in backend_result:
-            converted_key = BitString.from_binary(BitStringLSB.from_int(integer=k, nbits=self.n_qubits).binary)
-            if converted_key in result._state:
-                result._state[converted_key] += 1
-            else:
-                result._state[converted_key] = 1
+            result[k] += 1
 
         if target_qubits is not None:
             mapped_target = [self.qubit_map[q].number for q in target_qubits]
             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)
+            result = QubitWaveFunction.from_wavefunction(result, keymap, n_qubits=len(target_qubits))
 
         return result
 
@@ -211,9 +218,7 @@ def do_sample(self, samples, circuit, noise_model=None, initial_state=0, *args,
         QubitWaveFunction:
             the results of sampling, as a Qubit Wave Function.
         """
-        state = self.initialize_state(self.n_qubits)
-        lsb = BitStringLSB.from_int(initial_state, nbits=self.n_qubits)
-        state.set_computational_basis(BitString.from_binary(lsb.binary).integer)
+        state = self.initialize_state(self.n_qubits, initial_state)
         circuit.update_quantum_state(state)
         sampled = state.sampling(samples)
         return self.convert_measurements(backend_result=sampled, target_qubits=self.measurements)
diff --git a/src/tequila/simulators/simulator_qulacs_gpu.py b/src/tequila/simulators/simulator_qulacs_gpu.py
index 76dcb66e..41999b0b 100755
--- a/src/tequila/simulators/simulator_qulacs_gpu.py
+++ b/src/tequila/simulators/simulator_qulacs_gpu.py
@@ -1,16 +1,18 @@
 import qulacs
+from qulacs_core import QuantumStateGpu
+
 from tequila import TequilaException
 from tequila.simulators.simulator_qulacs import BackendCircuitQulacs, BackendExpectationValueQulacs
 
+
 class TequilaQulacsGpuException(TequilaException):
     def __str__(self):
         return "Error in qulacs gpu backend:" + self.message
 
+
 class BackendCircuitQulacsGpu(BackendCircuitQulacs):
-    def initialize_state(self, n_qubits:int=None) -> qulacs.QuantumState:
-        if n_qubits is None:
-            n_qubits = self.n_qubits
-        return qulacs.QuantumStateGpu(n_qubits)
+    quantum_state_class = QuantumStateGpu
+
 
 class BackendExpectationValueQulacsGpu(BackendExpectationValueQulacs):
     BackendCircuitType = BackendCircuitQulacsGpu
diff --git a/src/tequila/simulators/simulator_symbolic.py b/src/tequila/simulators/simulator_symbolic.py
index ddc3a24f..68003e93 100755
--- a/src/tequila/simulators/simulator_symbolic.py
+++ b/src/tequila/simulators/simulator_symbolic.py
@@ -42,37 +42,30 @@ def update_variables(self, variables):
 
     @classmethod
     def apply_gate(cls, state: QubitWaveFunction, gate: QGate, qubits: dict, variables) -> QubitWaveFunction:
-        result = QubitWaveFunction()
         n_qubits = len(qubits.keys())
+        result = sympy.Integer(1) * QubitWaveFunction(n_qubits)
         for s, v in state.items():
-            s.nbits = n_qubits
-            result += v * cls.apply_on_standard_basis(gate=gate, basisfunction=s, qubits=qubits, variables=variables)
+            result += v * cls.apply_on_standard_basis(gate=gate, basis_state=s, qubits=qubits, variables=variables)
         return result
 
     @classmethod
-    def apply_on_standard_basis(cls, gate: QGate, basisfunction: BitString, qubits:dict, variables) -> QubitWaveFunction:
-
-        basis_array = basisfunction.array
-        if gate.is_controlled():
-            do_apply = True
-            check = [basis_array[qubits[c]] == 1 for c in gate.control]
-            for c in check:
-                if not c:
-                    do_apply = False
-            if not do_apply:
-                return QubitWaveFunction.from_int(basisfunction)
+    def apply_on_standard_basis(cls, gate: QGate, basis_state: BitString, qubits:dict, variables) -> QubitWaveFunction:
+        n_qubits = len(qubits.keys())
+        basis_array = basis_state.array
+        if gate.is_controlled() and not all(basis_array[qubits[c]] == 1 for c in gate.control):
+            return QubitWaveFunction.from_basis_state(n_qubits, basis_state)
 
         if len(gate.target) > 1:
             raise Exception("Multi-targets not supported for symbolic simulators")
 
-        result = QubitWaveFunction()
+        result = sympy.Integer(1) * QubitWaveFunction(n_qubits)
         for tt in gate.target:
             t = qubits[tt]
             qt = basis_array[t]
             a_array = copy.deepcopy(basis_array)
             a_array[t] = (a_array[t] + 1) % 2
-            current_state = QubitWaveFunction.from_int(basisfunction)
-            altered_state = QubitWaveFunction.from_int(BitString.from_array(a_array))
+            current_state = QubitWaveFunction.from_basis_state(n_qubits, basis_state)
+            altered_state = QubitWaveFunction.from_basis_state(n_qubits, BitString.from_array(a_array))
 
             fac1 = None
             fac2 = None
@@ -115,22 +108,21 @@ def do_simulate(self, variables, initial_state: int = None, *args, **kwargs) ->
         count = 0
         for q in self.abstract_circuit.qubits:
             qubits[q] = count
-            count +=1
+            count += 1
 
         n_qubits = len(self.abstract_circuit.qubits)
 
         if initial_state is None:
-            initial_state = QubitWaveFunction.from_int(i=0, n_qubits=n_qubits)
-        elif isinstance(initial_state, int):
-            initial_state = QubitWaveFunction.from_int(initial_state, n_qubits=n_qubits)
+            initial_state = 0
+        initial_state = QubitWaveFunction.from_basis_state(n_qubits, initial_state)
 
         result = initial_state
         for g in self.abstract_circuit.gates:
             result = self.apply_gate(state=result, gate=g, qubits=qubits, variables=variables)
 
-        wfn = QubitWaveFunction()
+        wfn = QubitWaveFunction(n_qubits)
         if self.convert_to_numpy:
-            for k,v in result.items():
+            for k, v in result.items():
                 wfn[k] = complex(v)
         else:
             wfn = result
diff --git a/src/tequila/utils/bitstrings.py b/src/tequila/utils/bitstrings.py
index 92c3aa9c..9ae7cb20 100644
--- a/src/tequila/utils/bitstrings.py
+++ b/src/tequila/utils/bitstrings.py
@@ -178,7 +178,7 @@ def numbering(self) -> BitNumbering:
         return BitNumbering.LSB
 
 
-def _reverse_int_bits(x: int, nbits: int) -> int:
+def reverse_int_bits(x: int, nbits: int) -> int:
     if nbits is None:
         nbits = x.bit_length()
     assert nbits <= 32
@@ -192,8 +192,13 @@ def _reverse_int_bits(x: int, nbits: int) -> int:
 
 
 def initialize_bitstring(integer: int, nbits: int = None, numbering_in: BitNumbering = BitNumbering.MSB,
-                         numbering_out: BitNumbering = BitNumbering.MSB):
-    integer = _reverse_int_bits(integer, nbits) if numbering_in != numbering_out else integer
+                         numbering_out: BitNumbering = None):
+    if numbering_out is None:
+        numbering_out = numbering_in
+
+    if numbering_in != numbering_out:
+        integer = reverse_int_bits(integer, nbits)
+
     if numbering_out == BitNumbering.MSB:
         return BitString.from_int(integer=integer, nbits=nbits)
     else:
diff --git a/src/tequila/version.py b/src/tequila/version.py
index 6d6a8e74..0d2d4787 100644
--- a/src/tequila/version.py
+++ b/src/tequila/version.py
@@ -1,2 +1,2 @@
-__version__ = "1.9.7"
+__version__ = "1.9.8"
 __author__ = "Tequila Developers "
diff --git a/src/tequila/wavefunction/qubit_wavefunction.py b/src/tequila/wavefunction/qubit_wavefunction.py
index 2ad75bd4..e7f48ad9 100644
--- a/src/tequila/wavefunction/qubit_wavefunction.py
+++ b/src/tequila/wavefunction/qubit_wavefunction.py
@@ -1,296 +1,410 @@
 from __future__ import annotations
-import copy
+
 import typing
-from typing import Dict, Union
+from copy import deepcopy
+from math import log2
+from typing import Union, Generator
+
 import numpy
+import numpy as np
+import numpy.typing as npt
 import numbers
 
-from tequila.utils.bitstrings import BitNumbering, BitString, initialize_bitstring
-from tequila import TequilaException
-from tequila.utils.keymap import KeyMapLSB2MSB, KeyMapMSB2LSB
-from tequila.tools import number_to_string
+import sympy
 
-# from __future__ import annotations # can use that in python 3.7+ to get rid of string type hints
+from tequila.utils.bitstrings import BitString, reverse_int_bits
+from tequila import TequilaException, BitNumbering, initialize_bitstring
+from tequila.utils.keymap import KeyMapABC
 
 if typing.TYPE_CHECKING:
-    # don't need those structures, just for convenient type hinting
+    # Don't need those structures, just for convenient type hinting
     from tequila.hamiltonian.qubit_hamiltonian import QubitHamiltonian, PauliString
 
 
 class QubitWaveFunction:
     """
-    Store Wavefunction as dictionary of comp. basis state and complex numbers
-    Use the same structure for Measurments results with int instead of complex numbers (counts)
+    Represents a wavefunction.
+    Amplitudes are either stored in a Numpy array for dense wavefunctions, or a dictionary for sparse wavefunctions.
+    Does not enforce normalization.
     """
 
-    numbering = BitNumbering.MSB
-
-    def apply_keymap(self, keymap, initial_state: BitString = None):
-        self.n_qubits = keymap.n_qubits
-        mapped_state = dict()
-        for k, v in self.state.items():
-            mapped_key=keymap(input_state=k, initial_state=initial_state)
-            if mapped_key in mapped_state:
-                mapped_state[mapped_key] += v
-            else:
-                mapped_state[mapped_key] = v
-
-        self.state = mapped_state
-        return self
+    def __init__(self, n_qubits: int, numbering: BitNumbering = BitNumbering.MSB, dense: bool = False,
+                 init_state: bool = True) -> None:
+        """
+        Initialize a QubitWaveFunction with all amplitudes set to zero.
 
-    @property
-    def n_qubits(self) -> int:
-        if self._n_qubits is None:
-            return self.min_qubits()
+        :param n_qubits: Number of qubits.
+        :param numbering: Whether the first qubit is the most or least significant.
+        :param dense: Whether to store the amplitudes in a Numpy array instead of a dictionary.
+        :param init_state: Whether to initialize the state array.
+            If False, set_state must be called immediately after the constructor.
+        """
+        self._n_qubits: int = n_qubits
+        self._numbering = numbering
+        self._dense = dense
+        if init_state:
+            self._state = np.zeros(2 ** self._n_qubits, dtype=complex) if dense else dict()
         else:
-            return max(self._n_qubits, self.min_qubits())
+            self._state = None
 
-    def min_qubits(self) -> int:
-        if len(self.state) > 0:
-            maxk = max(self.state.keys())
-            return maxk.nbits
+    @classmethod
+    def from_wavefunction(cls, wfn: QubitWaveFunction, keymap: KeyMapABC = None, n_qubits: int = None,
+                          initial_state: BitString = None) -> QubitWaveFunction:
+        """
+        Create a copy of a wavefunction.
+
+        :param wfn: The wavefunction to copy.
+        :param keymap: A keymap to apply to the wavefunction.
+        :param n_qubits: Number of qubits of the new wavefunction.
+            Must not be None if keymap is not None.
+        :param initial_state: Initial state to pass to the keymap.
+        :return: The copied wavefunction.
+        """
+        if keymap is not None:
+            result = QubitWaveFunction(n_qubits, numbering=wfn._numbering, dense=wfn._dense)
+            # Change amplitudes to sympy objects
+            if wfn._dense and wfn._state.dtype == object:
+                result = sympy.Integer(1) * result
+            for index, coeff in wfn.raw_items():
+                key = initialize_bitstring(index, wfn._n_qubits, numbering_in=wfn._numbering, numbering_out=keymap.numbering)
+                key = keymap(key, initial_state)
+                result[key] += coeff
+            return result
         else:
-            return 0
+            return deepcopy(wfn)
 
-    @n_qubits.setter
-    def n_qubits(self, n_qubits):
-        if n_qubits is not None:
-            self._n_qubits = max(n_qubits, self.min_qubits())
-        return self
-
-    @property
-    def state(self):
-        if self._state is None:
-            return dict()
-        else:
-            return self._state
-
-    @state.setter
-    def state(self, other: Dict[BitString, complex]):
-        assert (isinstance(other, dict))
-        self._state = other
-
-    def __init__(self, state: Dict[BitString, complex] = None, n_qubits=None):
-        if state is None:
-            self._state = dict()
-        elif isinstance(state, int):
-            self._state = self.from_int(i=state, n_qubits=n_qubits).state
-        elif isinstance(state, str):
-            self._state = self.from_string(string=state, n_qubits=n_qubits).state
-        elif isinstance(state, numpy.ndarray) or isinstance(state, list):
-            self._state = self.from_array(arr=state, n_qubits=n_qubits).state
-        elif hasattr(state, "state"):
-            self._state = state.state
-        else:
-            self._state = state
-        self._n_qubits = n_qubits
+    @classmethod
+    def from_array(cls, array: npt.NDArray[complex], numbering: BitNumbering = BitNumbering.MSB,
+                   copy: bool = True) -> QubitWaveFunction:
+        """
+        Create a dense wavefunction from a Numpy array.
 
-    def items(self):
-        return self.state.items()
+        :param array: Array of amplitudes.
+        :param numbering: Whether the first qubit is the most or least significant.
+        :param copy: Whether to copy the array or use it directly.
+            If False, the array must not be modified after the constructor.
+        :return: The created wavefunction.
+        """
+        if not log2(len(array)).is_integer():
+            raise ValueError(f"Array length must be a power of 2, received {len(array)}")
+        n_qubits = int(log2(len(array)))
+        result = QubitWaveFunction(n_qubits, numbering, dense=True, init_state=False)
+        result.set_state(array, copy)
+        return result
 
-    def keys(self):
-        return self.state.keys()
+    @classmethod
+    def from_basis_state(cls, n_qubits: int, basis_state: Union[int, BitString],
+                         numbering: BitNumbering = BitNumbering.MSB) -> QubitWaveFunction:
+        """
+        Create a sparse wavefunction that is a basis state.
 
-    def values(self):
-        return self.state.values()
+        :param n_qubits: Number of qubits.
+        :param basis_state: Index of the basis state.
+        :param numbering: Whether the first qubit is the most or least significant.
+        :return: The created wavefunction.
+        """
+        if 2 ** n_qubits <= basis_state:
+            raise ValueError(f"Number of qubits {n_qubits} insufficient for basis state {basis_state}")
+        if isinstance(basis_state, BitString):
+            basis_state = reverse_int_bits(basis_state.integer,
+                                           basis_state.nbits) if numbering != basis_state.numbering else basis_state.integer
+        result = QubitWaveFunction(n_qubits, numbering)
+        result[basis_state] = 1.0
+        return result
 
-    @staticmethod
-    def convert_bitstring(key: Union[BitString, numbers.Integral], n_qubits):
-        if isinstance(key, numbers.Integral):
-            return BitString.from_int(integer=key, nbits=n_qubits)
-        elif isinstance(key, str):
-            return BitString.from_binary(binary=key, nbits=n_qubits)
-        else:
-            return key
+    @classmethod
+    def from_string(cls, string: str, numbering: BitNumbering = BitNumbering.MSB) -> QubitWaveFunction:
+        """
+        Create a sparse wavefunction from a string.
 
-    def __getitem__(self, item: BitString):
-        key = self.convert_bitstring(item, self.n_qubits)
-        return self.state[key]
+        :param string: String representation of the wavefunction.
+        :param numbering: Whether the first qubit is the most or least significant.
+        :return: The created wavefunction.
+        """
+        try:
+            string = string.replace(" ", "")
+            string = string.replace("*", "")
+            terms = string.split(">")[:-1]
+            n_qubits = len(terms[0].split("|")[-1])
+            result = QubitWaveFunction(n_qubits, numbering)
+            for term in terms:
+                coeff, index = term.split("|")
+                coeff = complex(coeff) if coeff != "" else 1.0
+                index = int(index, 2)
+                result[index] = coeff
+            return result
+        except ValueError:
+            raise TequilaException(f"Failed to initialize QubitWaveFunction from string:\n\"{string}\"\n")
 
-    def __call__(self, key, *args, **kwargs) -> numbers.Number:
+    @classmethod
+    def convert_from(cls, n_qubits: int, val: Union[QubitWaveFunction, int, str, numpy.ndarray]):
         """
-        Like getitem but returns zero if key is not there
-
-        Parameters
-        ----------
-        key: bitstring (or int or str)
-        Returns
-        -------
-            Return the amplitude or measurement occurence of a bitstring
+        Convert a value to a QubitWaveFunction.
+        Accepts QubitWaveFunction, int, str, and numpy.ndarray.
+
+        :param n_qubits: Number of qubits.
+        :param val: Value to convert.
+        :return: The converted value.
         """
-        ckey = self.convert_bitstring(key, self.n_qubits)
-        if ckey in self.state:
-            return self.state[ckey]
+        if isinstance(val, QubitWaveFunction):
+            return val
+        elif isinstance(val, int):
+            return cls.from_basis_state(n_qubits=n_qubits, basis_state=val)
+        elif isinstance(val, str):
+            return cls.from_string(val)
+        elif isinstance(val, numpy.ndarray):
+            return cls.from_array(val)
         else:
-            return 0.0
+            raise TequilaException(f"Cannot initialize QubitWaveFunction from type {type(val)}")
 
+    @property
+    def n_qubits(self) -> int:
+        """
+        Returns number of qubits in the wavefunction.
+        """
+        return self._n_qubits
 
+    @property
+    def numbering(self) -> BitNumbering:
+        """
+        Returns the bit numbering of the wavefunction.
+        """
+        return self._numbering
 
-    def __setitem__(self, key: BitString, value: numbers.Number):
-        self._state[self.convert_bitstring(key, self.n_qubits)] = value
-        return self
-
-    def __contains__(self, item: BitString):
-        return self.convert_bitstring(item, self.n_qubits) in self.keys()
-
-    def __len__(self):
-        return len(self.state)
-
-    @classmethod
-    def from_array(cls, arr: numpy.ndarray, keymap=None, threshold: float = 1.e-6,
-                   numbering: BitNumbering = BitNumbering.MSB, n_qubits: int = None):
-        arr = numpy.asarray(arr)
-        assert (len(arr.shape) == 1)
-        state = dict()
-        maxkey = len(arr) - 1
-        maxbit = initialize_bitstring(integer=maxkey, numbering_in=numbering, numbering_out=cls.numbering).nbits
-        for ii, v in enumerate(arr):
-            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)
-
-        return result
+    @property
+    def dense(self) -> bool:
+        """
+        Returns whether the wavefunction is dense.
+        """
+        return self._dense
 
-    @classmethod
-    def from_int(cls, i: int, coeff=1, n_qubits: int = None):
-        if isinstance(i, BitString):
-            return QubitWaveFunction(state={i: coeff}, n_qubits=n_qubits)
+    def to_array(self, out_numbering: BitNumbering = BitNumbering.MSB, copy: bool = True) -> npt.NDArray[complex]:
+        """
+        Returns array of amplitudes.
+
+        :param out_numbering: Whether the first qubit is the most or least significant in the output array indices.
+            For dense wavefunctions, this operation is significantly cheaper when this is the same as the numbering
+            of the wavefunction.
+        :param copy: Whether to copy the array or use it directly for dense Wavefunctions.
+            If False, changes to the array or wavefunction will affect each other.
+        :return: Array of amplitudes.
+        """
+        if self._dense and self._numbering == out_numbering:
+            return self._state.copy() if copy else self._state
         else:
-            return QubitWaveFunction(state={BitString.from_int(integer=i, nbits=n_qubits): coeff}, n_qubits=n_qubits)
+            result = np.zeros(2 ** self._n_qubits, dtype=complex)
+            for k, v in self.raw_items():
+                if self._numbering != out_numbering:
+                    k = reverse_int_bits(k, self._n_qubits)
+                result[k] = v
+            return result
+
+    def set_state(self, value: npt.NDArray[complex], copy: bool = True) -> None:
+        """
+        Sets the state to an array.
+        After this call, the wavefunction will be dense.
 
-    @classmethod
-    def from_string(cls, string: str, n_qubits: int = None):
+        :param value: Array of amplitudes. Length must be 2 ** n_qubits.
+        :param copy: Whether to copy the array or use it directly.
+            If False, changes to the array or wavefunction will affect each other.
         """
-        Complex values like (x+iy)|...> will currently not work, you need to type Real and imaginary separately
-        Or improve this constructor :-)
-        e.g instead of (0.5+1.0j)|0101> do 0.5|0101> + 1.0j|0101>
-        :param paths:
-        :param string:
-        :return:
+        if len(value) != 2 ** self._n_qubits:
+            raise ValueError(f"Wavefunction of {self._n_qubits} qubits must have {2 ** self._n_qubits} amplitudes, "
+                             f"received {len(value)}")
+        self._dense = True
+        if copy:
+            self._state = value.copy()
+        else:
+            self._state = value
+
+    def __getitem__(self, key: Union[int, BitString]) -> complex:
+        if isinstance(key, BitString):
+            key = reverse_int_bits(key.integer, key.nbits) if self._numbering != key.numbering else key.integer
+        return self._state[key] if self._dense else self._state.get(key, 0)
+
+    def __setitem__(self, key: Union[int, BitString], value: complex) -> None:
+        if isinstance(key, BitString):
+            key = reverse_int_bits(key.integer, key.nbits) if self._numbering != key.numbering else key.integer
+        self._state[key] = value
+
+    def __contains__(self, item: Union[int, BitString]) -> bool:
+        if isinstance(item, BitString):
+            item = reverse_int_bits(item.integer, item.nbits) if self._numbering != item.numbering else item.integer
+        return abs(self[item]) > 1e-6
+
+    def raw_items(self) -> Generator[tuple[int, complex]]:
+        """Returns a generator of non-zero amplitudes with integer indices."""
+        return ((k, v) for k, v in (enumerate(self._state) if self._dense else self._state.items()))
+
+    def items(self) -> Generator[tuple[BitString, complex]]:
+        """Returns a generator of non-zero amplitudes with BitString indices."""
+        return ((initialize_bitstring(k, self._n_qubits, self._numbering), v)
+                for k, v in self.raw_items()
+                if isinstance(v, sympy.Basic) or abs(v) > 1e-6)
+
+    def keys(self) -> Generator[BitString]:
+        """Returns a generator of BitString indices of non-zero amplitudes."""
+        return (k for k, v in self.items())
+
+    def values(self) -> Generator[complex]:
+        """Returns a generator of non-zero amplitudes."""
+        return (v for k, v in self.items())
+
+    def __eq__(self, other) -> bool:
+        if not isinstance(other, QubitWaveFunction):
+            return False
+
+        raise TequilaException("Wavefunction equality is not well-defined. Consider using isclose.")
+
+    def isclose(self: QubitWaveFunction,
+                other: QubitWaveFunction,
+                rtol: float = 1e-5,
+                atol: float = 1e-8) -> bool:
         """
-        try:
-            state = dict()
-            string = string.replace(" ", "")
-            string = string.replace("*", "")
-            string = string.replace("+-", "-")
-            string = string.replace("-+", "-")
-            terms = (string + "terminate").split('>')
-            for term in terms:
-                if term == 'terminate':
-                    break
-                tmp = term.split("|")
-                coeff = tmp[0]
-                if coeff == '':
-                    coeff = 1.0
-                else:
-                    coeff = complex(coeff)
-                basis_state = BitString.from_binary(binary=tmp[1])
+        Check if two wavefunctions are close, up to a global phase.
 
-                state[basis_state] = coeff
-        except ValueError:
-            raise TequilaException("Failed to initialize QubitWaveFunction from string:" + string + "\n"
-                                                                                                    "did you try complex values?\n"
-                                                                                                    "currently you need to type real and imaginary parts separately\n"
-                                                                                                    "e.g. instead of (0.5+1.0j)|0101> do 0.5|0101> + 1.0j|0101>")
-        except:
-            raise TequilaException("Failed to initialize QubitWaveFunction from string:" + string)
-        return QubitWaveFunction(state=state, n_qubits=n_qubits)
+        :param other: The other wavefunction.
+        :param rtol: Relative tolerance.
+        :param atol: Absolute tolerance.
+        :return: Whether the wavefunctions are close.
+        """
+        inner = self.inner(other)
+        self_norm = self.norm()
+        other_norm = other.norm()
+        cosine_similarity = inner / (self_norm * other_norm)
 
-    def __repr__(self):
-        result = str()
-        for k, v in self.items():
-            result += number_to_string(number=v) + "|" + str(k.binary) + "> "
-        return result
+        return (np.isclose(abs(cosine_similarity), 1.0, rtol, atol)
+                and np.isclose(self_norm, other_norm, rtol, atol))
 
-    def __eq__(self, other):
-        raise TequilaException("Wavefunction equality is not well-defined. Consider using inner"
-                               + " product equality, wf1.isclose(wf2).")
-
-    def isclose(self :  'QubitWaveFunction',
-                other : 'QubitWaveFunction',
-                rtol : float=1e-5,
-                atol : float=1e-8) -> bool:
-        """Return whether this wavefunction is similar to the target wavefunction."""
-        over1 = complex(self.inner(other))
-        over2 = numpy.sqrt(complex(self.inner(self) * other.inner(other)))
-        # Explicit casts to complex() is required if self or other are sympy
-        # wavefunction with sympy-typed amplitudes
-
-        # Check if the two numbers are equal.
-        return numpy.isclose(over1, over2, rtol=rtol, atol=atol)
-
-    def __add__(self, other):
-        result = QubitWaveFunction(state=copy.deepcopy(self._state))
-        for k, v in other.items():
-            if k in result._state:
-                result._state[k] += v
-            else:
-                result._state[k] = v
-        return result
+    def __add__(self, other: QubitWaveFunction) -> QubitWaveFunction:
+        if self._dense and other._dense and self._numbering == other._numbering:
+            return QubitWaveFunction.from_array(self._state + other._state, self.numbering, copy=False)
+        else:
+            result = QubitWaveFunction.from_wavefunction(self)
+            result += other
+            return result
 
-    def __sub__(self, other):
-        return self + -1.0 * other
+    def __iadd__(self, other: QubitWaveFunction) -> QubitWaveFunction:
+        if self._dense and other._dense and self._numbering == other._numbering:
+            self._state += other._state
+        else:
+            for k, v in other.raw_items():
+                if self._numbering != other._numbering:
+                    k = reverse_int_bits(k, self._n_qubits)
+                self[k] += v
+        return self
+
+    def __sub__(self, other: QubitWaveFunction) -> QubitWaveFunction:
+        if self._dense and other._dense and self._numbering == other._numbering:
+            return QubitWaveFunction.from_array(self._state - other._state, self.numbering, copy=False)
+        else:
+            result = QubitWaveFunction.from_wavefunction(self)
+            result -= other
+            return result
 
-    def __iadd__(self, other):
-        for k, v in other.items():
-            if k in self._state:
-                self._state[k] += v
-            else:
-                self._state[k] = v
+    def __isub__(self, other: QubitWaveFunction) -> QubitWaveFunction:
+        if self._dense and other._dense and self._numbering == other._numbering:
+            self._state -= other._state
+        else:
+            for k, v in other.raw_items():
+                if self._numbering != other._numbering:
+                    k = reverse_int_bits(k, self._n_qubits)
+                self[k] -= v
         return self
 
-    def __rmul__(self, other):
-        result = QubitWaveFunction(state=copy.deepcopy(self._state))
-        for k, v in result._state.items():
-            result._state[k] *= other
-        return result
+    def __rmul__(self, other: complex) -> QubitWaveFunction:
+        if self._dense:
+            return QubitWaveFunction.from_array(other * self._state, self.numbering, copy=False)
+        else:
+            result = QubitWaveFunction.from_wavefunction(self)
+            result *= other
+            return result
 
-    def inner(self, other):
-        # currently very slow and not optimized in any way
-        result = 0.0
-        for k, v in self.items():
-            if k in other._state:
-                result += v.conjugate() * other._state[k]
-        return result
+    def __imul__(self, other: complex) -> QubitWaveFunction:
+        if self._dense:
+            self._state *= other
+        else:
+            for k, v in self.raw_items():
+                self[k] = other * v
+        return self
 
-    def normalize(self):
+    def inner(self, other: QubitWaveFunction) -> complex:
+        """Returns the inner product with another wavefunction."""
+        if self._dense and other._dense and self._numbering == other._numbering:
+            return np.inner(self._state.conjugate(), other._state)
+        else:
+            result = 0
+            for k, v in self.raw_items():
+                if self._numbering != other._numbering:
+                    k = reverse_int_bits(k, self._n_qubits)
+                result += v.conjugate() * other[k]
+            if isinstance(result, sympy.Basic):
+                result = complex(result)
+            return result
+
+    def norm(self) -> float:
+        """Returns the norm of the wavefunction."""
+        return np.sqrt(self.inner(self))
+
+    def normalize(self, inplace: bool = False) -> QubitWaveFunction:
         """
-        NOT AN Inplace operation
-        :return: Normalizes the wavefunction/countrate
+        Normalizes the wavefunction.
+
+        :param inplace: Whether to normalize the wavefunction in place or return a new one.
+        :return: The normalized wavefunction.
         """
-        norm2 = self.inner(other=self)
-        normalized = 1.0 / numpy.sqrt(norm2) * self
-        return normalized
+        norm = self.norm()
+        if inplace:
+            self *= 1.0 / norm
+            return self
+        else:
+            return (1.0 / norm) * QubitWaveFunction.from_wavefunction(self)
+
+    # It would be nice to call this __len__, however for some reason this causes infinite loops
+    # when multiplying wave functions with some types of numbers from the right sight, likely
+    # because the __mul__ implementation of the number tries to perform some sort of array
+    # operation.
+    def length(self):
+        return sum(1 for _ in self.raw_items())
+
+    def __repr__(self):
+        result = str()
+        for index, coeff in self.items():
+            index = index.integer
+            if self.numbering == BitNumbering.LSB:
+                index = reverse_int_bits(index, self._n_qubits)
+            result += f"{coeff} |{index:0{self._n_qubits}b}> "
+        # If the wavefunction contains no states
+        if not result:
+            result = "empty wavefunction"
+        return result
 
-    def compute_expectationvalue(self, operator: 'QubitHamiltonian') -> numbers.Real:
+    def compute_expectationvalue(self, operator: QubitHamiltonian) -> numbers.Real:
         tmp = self.apply_qubitoperator(operator=operator)
         E = self.inner(other=tmp)
-        if hasattr(E, "imag") and numpy.isclose(E.imag, 0.0, atol=1.e-6):
+        if hasattr(E, "imag") and np.isclose(E.imag, 0.0, atol=1.e-6):
             return float(E.real)
         else:
             return E
 
-    def apply_qubitoperator(self, operator: 'QubitHamiltonian'):
+    def apply_qubitoperator(self, operator: QubitHamiltonian) -> QubitWaveFunction:
         """
         Inefficient function which computes the action of a QubitHamiltonian on this wfn
         :param operator: QubitOperator
         :return: resulting Qubitwavefunction
         """
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self.n_qubits, self._numbering)
         for ps in operator.paulistrings:
             result += self.apply_paulistring(paulistring=ps)
-        result = result.simplify()
         return result
 
-    def apply_paulistring(self, paulistring: 'PauliString'):
+    def apply_paulistring(self, paulistring: PauliString) -> QubitWaveFunction:
         """
         Inefficient function which computes action of a single paulistring
         :param paulistring: PauliString
         :return: Expectation Value
         """
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self._n_qubits, self._numbering)
         for k, v in self.items():
             arr = k.array
             c = v
@@ -306,17 +420,3 @@ def apply_paulistring(self, paulistring: 'PauliString'):
                     raise TequilaException("unknown pauli: " + str(p))
             result[BitString.from_array(array=arr)] = c
         return paulistring.coeff * result
-
-    def to_array(self):
-        result = numpy.zeros(shape=2 ** self.n_qubits, dtype=complex)
-        for k, v in self.items():
-            result[int(k)] = v
-        return result
-
-    def simplify(self, threshold = 1.e-8):
-        state = {}
-        for k, v in self.state.items():
-            if not numpy.isclose(v, 0.0, atol=threshold):
-                state[k] = v
-        return QubitWaveFunction(state=state)
-


Test patch:
diff --git a/tests/test_binary_pauli.py b/tests/test_binary_pauli.py
index 499fe43f..ff3bdcc3 100644
--- a/tests/test_binary_pauli.py
+++ b/tests/test_binary_pauli.py
@@ -169,7 +169,7 @@ def test_commuting_groups():
 
 def prepare_cov_dict(H):
     eigenValues, eigenVectors = np.linalg.eigh(H.to_matrix())
-    wfn0 = tq.QubitWaveFunction(eigenVectors[:,0])
+    wfn0 = tq.QubitWaveFunction.from_array(eigenVectors[:, 0])
     terms = BinaryHamiltonian.init_from_qubit_hamiltonian(H).binary_terms
     cov_dict = {}
     for term1 in terms:
diff --git a/tests/test_circuits.py b/tests/test_circuits.py
index 05a9fb86..f93be5d9 100644
--- a/tests/test_circuits.py
+++ b/tests/test_circuits.py
@@ -11,6 +11,7 @@
 from tequila.simulators.simulator_api import simulate
 from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
 
+
 def test_qubit_map():
 
     for gate in [Rx, Ry, Rz]:
@@ -82,7 +83,7 @@ def test_conventions():
 
 
 def strip_sympy_zeros(wfn: QubitWaveFunction):
-    result = QubitWaveFunction()
+    result = QubitWaveFunction(wfn.n_qubits, wfn.numbering)
     for k, v in wfn.items():
         if v != 0:
             result[k] = v
@@ -130,7 +131,7 @@ def test_basic_gates():
     cos = sympy.cos
     sin = sympy.sin
     exp = sympy.exp
-    BS = QubitWaveFunction.from_int
+    def BS(state): return QubitWaveFunction.from_basis_state(n_qubits=1, basis_state=state)
     angle = sympy.pi
     gates = [X(0), Y(0), Z(0), Rx(target=0, angle=angle), Ry(target=0, angle=angle), Rz(target=0, angle=angle), H(0)]
     results = [
@@ -310,15 +311,17 @@ def test_unitary_gate_u_u3(gate, theta, phi, lambd):
 
 def test_swap():
     U = X(0)
-    U += SWAP(0,2)
+    U += SWAP(0, 2)
     wfn = simulate(U)
     wfnx = simulate(X(2))
     assert numpy.isclose(numpy.abs(wfn.inner(wfnx))**2,1.0)
 
     U = X(2)
-    U += SWAP(0,2, power=2.0)
+    U += SWAP(0,2, power=3.0)
     wfn = simulate(U)
-    wfnx = simulate(X(0))
+    V = X(0)
+    V.n_qubits = 3
+    wfnx = simulate(V)
     assert numpy.isclose(numpy.abs(wfn.inner(wfnx))**2,1.0)
 
     U = X(0)+X(3)
@@ -345,9 +348,13 @@ def test_givens():
     U = X(0)
     U += Givens(0, 1, angle=-numpy.pi/4)
     wfn = simulate(U)
-    wfnx = simulate(X(0))
+    V = X(0)
+    V.n_qubits = 2
+    wfnx = simulate(V)
     assert numpy.isclose(numpy.abs(wfn.inner(wfnx))**2, 0.5)
-    wfnx = simulate(X(1))
+    V = X(1)
+    V.n_qubits = 2
+    wfnx = simulate(V)
     assert numpy.isclose(numpy.abs(wfn.inner(wfnx))**2, 0.5)
     
     U = X(0)
diff --git a/tests/test_hamiltonian_arithmetic.py b/tests/test_hamiltonian_arithmetic.py
index 42870a56..2fffa75a 100644
--- a/tests/test_hamiltonian_arithmetic.py
+++ b/tests/test_hamiltonian_arithmetic.py
@@ -84,16 +84,16 @@ def test_initialization():
 def test_ketbra():
     ket = QubitWaveFunction.from_string("1.0*|00> + 1.0*|11>").normalize()
     operator = paulis.KetBra(ket=ket, bra="|00>")
-    result = operator*QubitWaveFunction.from_int(0, n_qubits=2)
-    assert(result.isclose(ket))
+    result = operator * QubitWaveFunction.from_basis_state(n_qubits=2, basis_state=0)
+    assert result.isclose(ket)
 
 @pytest.mark.parametrize("n_qubits", [1,2,3,5])
 def test_ketbra_random(n_qubits):
     ket = numpy.random.uniform(0.0, 1.0, 2**n_qubits)
-    bra = QubitWaveFunction.from_int(0, n_qubits=n_qubits)
+    bra = QubitWaveFunction.from_basis_state(n_qubits, 0)
     operator = paulis.KetBra(ket=ket, bra=bra)
     result = operator * bra
-    assert(result.isclose(QubitWaveFunction.from_array(ket)))
+    assert result.isclose(QubitWaveFunction.from_array(ket))
 
 def test_paulistring_conversion():
     X1 = QubitHamiltonian.from_string("X0", openfermion_format=True)
@@ -196,7 +196,7 @@ def test_projectors(qubits):
     real = numpy.random.uniform(0.0,1.0,2**qubits)
     imag = numpy.random.uniform(0.0,1.0,2**qubits)
     array = real + 1.j*imag
-    wfn = QubitWaveFunction.from_array(arr=array)
+    wfn = QubitWaveFunction.from_array(array=array)
     P = paulis.Projector(wfn=wfn.normalize())
     assert(P.is_hermitian())
     assert(wfn.apply_qubitoperator(P).isclose(wfn))
@@ -304,17 +304,17 @@ def test_simple_trace_out():
 
     H1 = QubitHamiltonian.from_string("1.0*X(0)*Z(1)*Z(5)*X(100)Y(50)")
     H2 = QubitHamiltonian.from_string("1.0*X(0)X(100)Y(50)")
-    assert H2 == H1.trace_out_qubits(qubits=[1,5], states=[QubitWaveFunction.from_string("1.0*|1>")]*2)
+    assert H2 == H1.trace_out_qubits(qubits=[1,5], states=[QubitWaveFunction.from_string("1.0*|1>")] * 2)
 
     H1 = QubitHamiltonian.from_string("1.0*X(0)*Z(1)*X(100)Y(50)")
     H2 = QubitHamiltonian.from_string("-1.0*X(0)X(100)Y(50)")
-    assert H2 == H1.trace_out_qubits(qubits=[1,5], states=[QubitWaveFunction.from_string("1.0*|1>")]*2)
+    assert H2 == H1.trace_out_qubits(qubits=[1,5], states=[QubitWaveFunction.from_string("1.0*|1>")] * 2)
 
 @pytest.mark.parametrize("theta", numpy.random.uniform(0.0, 6.0, 10))
 def test_trace_out_xy(theta):
     a = numpy.sin(theta)
     b = numpy.cos(theta)
-    state = QubitWaveFunction.from_array([a,b])
+    state = QubitWaveFunction.from_array([a, b])
 
     H1 = QubitHamiltonian.from_string("1.0*X(0)*X(1)*X(100)")
     H2 = QubitHamiltonian.from_string("1.0*X(0)*X(100)")
diff --git a/tests/test_mappings.py b/tests/test_mappings.py
index b151201d..5659831a 100644
--- a/tests/test_mappings.py
+++ b/tests/test_mappings.py
@@ -12,15 +12,15 @@
 
 
 def test_keymaps():
-    initial_state = 0
+    initial_state = BitString.from_int(0)
 
-    small = QubitWaveFunction.from_int(i=int("0b1111", 2))
-    large = QubitWaveFunction.from_int(i=int("0b01010101", 2))
-    large.n_qubits = 8
+    small = QubitWaveFunction.from_basis_state(4, basis_state=int("0b1111", 2))
+    large = QubitWaveFunction.from_basis_state(8, basis_state=int("0b01010101", 2))
 
     keymap = KeyMapSubregisterToRegister(register=[0, 1, 2, 3, 4, 5, 6, 7], subregister=[1, 3, 5, 7])
 
-    assert (small.apply_keymap(keymap=keymap, initial_state=initial_state).isclose(large))
+    mapped = QubitWaveFunction.from_wavefunction(small, keymap, 8, initial_state)
+    assert mapped.isclose(large)
 
 
 def test_endianness():
@@ -47,6 +47,7 @@ def test_endianness():
         assert (i1 == BitString.from_bitstring(i22))
         assert (i2 == BitString.from_bitstring(i11))
 
+
 @pytest.mark.skipif(condition="cirq" not in INSTALLED_SAMPLERS or "qiskit" not in INSTALLED_SAMPLERS, reason="need cirq and qiskit for cross validation")
 def test_endianness_simulators():
     tests = ["000111",
@@ -73,7 +74,7 @@ def test_endianness_simulators():
         print("counts_cirq  =", counts_cirq)
         print("counts_qiskit=", counts_qiskit)
         assert (counts_cirq.isclose(counts_qiskit))
-        assert (wfn_cirq.state == counts_cirq.state)
+        assert (wfn_cirq.to_array() == counts_cirq.to_array())
 
 
 @pytest.mark.parametrize("backend", INSTALLED_SAMPLERS)
@@ -83,7 +84,7 @@ def test_paulistring_sampling(backend, case):
     H = QubitHamiltonian.from_paulistrings(PauliString.from_string(case[0]))
     U = gates.X(target=1) + gates.X(target=3) + gates.X(target=5)
     E = ExpectationValue(H=H, U=U)
-    result = simulate(E,backend=backend, samples=1)
+    result = simulate(E, backend=backend, samples=1)
     assert isclose(result, case[1], 1.e-4)
 
 
@@ -93,12 +94,13 @@ def test_paulistring_sampling_2(backend, case):
     H = QubitHamiltonian.from_paulistrings(PauliString.from_string(case[0]))
     U = gates.H(target=1) + gates.H(target=3) + gates.X(target=5) + gates.H(target=5)
     E = ExpectationValue(H=H, U=U)
-    result = simulate(E,backend=backend, samples=1)
-    assert (isclose(result, case[1], 1.e-4))
+    result = simulate(E, backend=backend, samples=1)
+    assert isclose(result, case[1], 1.e-4)
+
 
 @pytest.mark.parametrize("array", [numpy.random.uniform(0.0,1.0,i) for i in [2,4,8,16]])
 def test_qubitwavefunction_array(array):
     print(array)
-    wfn = QubitWaveFunction.from_array(arr=array)
+    wfn = QubitWaveFunction.from_array(array=array)
     array2 = wfn.to_array()
-    assert numpy.allclose(array,array2)
+    assert numpy.allclose(array, array2)
diff --git a/tests/test_simulator_backends.py b/tests/test_simulator_backends.py
index 7b004b3f..367aad40 100644
--- a/tests/test_simulator_backends.py
+++ b/tests/test_simulator_backends.py
@@ -4,11 +4,13 @@
 
 import importlib
 import numpy
+import numpy as np
 import pytest
 import random
 import numbers
 import tequila as tq
 import tequila.simulators.simulator_api
+from tequila import BitString, QubitWaveFunction
 
 """
 Warn if Simulators are not installed
@@ -40,7 +42,7 @@ def teardown_function(function):
 @pytest.mark.dependencies
 def test_dependencies():
     for package in tequila.simulators.simulator_api.SUPPORTED_BACKENDS:
-        if package != "qulacs_gpu":
+        if package not in ["qulacs_gpu", "qiskit_gpu"]:
             assert (package in tq.simulators.simulator_api.INSTALLED_BACKENDS)
 
 
@@ -275,7 +277,7 @@ def test_wfn_multitarget(simulator):
 @pytest.mark.parametrize("simulator", tequila.simulators.simulator_api.INSTALLED_SIMULATORS.keys())
 def test_wfn_multi_control(simulator):
     # currently no compiler, so that test can not succeed
-    if simulator == 'qiskit':
+    if simulator in ["qiskit", "qiskit_gpu"]:
         return
     ac = tq.gates.X([0, 1, 2])
     ac += tq.gates.Ry(target=[0], control=[1, 2], angle=2.3 / 2)
@@ -351,8 +353,27 @@ def test_initial_state_from_integer(simulator, initial_state):
         U += tq.gates.X(target=i) + tq.gates.X(target=i)
 
     wfn = tq.simulate(U, initial_state=initial_state, backend=simulator)
-    assert (initial_state in wfn)
-    assert (numpy.isclose(wfn[initial_state], 1.0))
+    assert BitString.from_int(initial_state, 6) in wfn
+    assert numpy.isclose(wfn[BitString.from_int(initial_state, 6)], 1.0)
+
+
+@pytest.mark.parametrize("simulator", tequila.simulators.simulator_api.INSTALLED_SIMULATORS.keys())
+def test_initial_state_from_wavefunction(simulator):
+    if not tequila.simulators.simulator_api.INSTALLED_SIMULATORS[simulator][0].supports_generic_initialization:
+        return
+
+    U = tq.gates.H(target=0)
+
+    state = QubitWaveFunction.from_array(np.array([1.0, 1.0])).normalize()
+    result = tq.simulate(U, initial_state=state, backend=simulator)
+    assert result.isclose(QubitWaveFunction.from_basis_state(n_qubits=1, basis_state=0))
+    result = tq.simulate(U, initial_state=state, backend=simulator, samples=100)
+    assert result.isclose(QubitWaveFunction.from_array(np.array([100.0, 0.0])))
+
+    state = QubitWaveFunction.from_array(np.array([1.0, -1.0])).normalize()
+    result = tq.simulate(U, initial_state=state, backend=simulator, samples=100)
+    assert result.isclose(QubitWaveFunction.from_array(np.array([0.0, 100.0])))
+
 
 
 @pytest.mark.parametrize("backend", tequila.simulators.simulator_api.INSTALLED_SIMULATORS.keys())
@@ -388,7 +409,7 @@ def test_sampling_read_out_qubits(backend):
     U = tq.gates.X(0)
     U += tq.gates.Z(1)
 
-    wfn = tq.QubitWaveFunction(2)
+    wfn = tq.QubitWaveFunction.from_basis_state(n_qubits=2, basis_state=BitString.from_int(2))
 
     result = tq.simulate(U, backend=backend, samples=1, read_out_qubits=[0, 1])
     assert (numpy.isclose(numpy.abs(wfn.inner(result)) ** 2, 1.0, atol=1.e-4))
diff --git a/tests/test_symbolic_simulator.py b/tests/test_symbolic_simulator.py
index 0da79126..6b5f867e 100644
--- a/tests/test_symbolic_simulator.py
+++ b/tests/test_symbolic_simulator.py
@@ -6,14 +6,14 @@
 import numpy
 import pytest
 
-@pytest.mark.parametrize("paulis", [(gates.X,paulis.X), (gates.Y, paulis.Y), (gates.Z,paulis.Z)])
-@pytest.mark.parametrize("qubit", [0,1,2])
-@pytest.mark.parametrize("init", [0,1])
+@pytest.mark.parametrize("paulis", [(gates.X, paulis.X), (gates.Y, paulis.Y), (gates.Z, paulis.Z)])
+@pytest.mark.parametrize("qubit", [0, 1, 2])
+@pytest.mark.parametrize("init", [0, 1])
 def test_pauli_gates(paulis, qubit, init):
-    iwfn = QubitWaveFunction.from_int(i=init, n_qubits=qubit+1)
+    iwfn = QubitWaveFunction.from_basis_state(n_qubits=qubit + 1, basis_state=init)
+    iwfn = iwfn.apply_qubitoperator(paulis[1](qubit))
     wfn = simulate(paulis[0](qubit), initial_state=init)
-    iwfn=iwfn.apply_qubitoperator(paulis[1](qubit))
-    assert(iwfn.isclose(wfn))
+    assert iwfn.isclose(wfn)
 
 @pytest.mark.parametrize("rot", [(gates.Rx,paulis.X), (gates.Ry, paulis.Y), (gates.Rz,paulis.Z)])
 @pytest.mark.parametrize("angle", numpy.random.uniform(0.0, 2*numpy.pi, 5))
@@ -22,7 +22,7 @@ def test_pauli_gates(paulis, qubit, init):
 def test_rotations(rot, qubit, angle, init):
     pauli = rot[1](qubit)
     gate = rot[0](target=qubit, angle=angle)
-    iwfn = QubitWaveFunction.from_int(i=init, n_qubits=qubit+1)
+    iwfn = QubitWaveFunction.from_basis_state(basis_state=init, n_qubits=qubit + 1)
     wfn = simulate(gate, initial_state=init)
     test= numpy.cos(-angle/2.0)*iwfn + 1.0j*numpy.sin(-angle/2.0)* iwfn.apply_qubitoperator(pauli)
     assert(wfn.isclose(test))
@@ -31,7 +31,7 @@ def test_rotations(rot, qubit, angle, init):
 @pytest.mark.parametrize("init", [0,1])
 def test_hadamard(qubit, init):
     gate = gates.H(target=qubit)
-    iwfn = QubitWaveFunction.from_int(i=init, n_qubits=qubit+1)
+    iwfn = QubitWaveFunction.from_basis_state(basis_state=init, n_qubits=qubit + 1)
     wfn = simulate(gate, initial_state=init)
     test= 1.0/numpy.sqrt(2)*(iwfn.apply_qubitoperator(paulis.Z(qubit)) + iwfn.apply_qubitoperator(paulis.X(qubit)))
     assert(wfn.isclose(test))
@@ -47,11 +47,8 @@ def test_controls(target, control, gate):
     assert(wfn0.isclose(wfn1))
 
     c0 = gates.QCircuit()
+    c0.n_qubits = max(target, control) + 1
     c1 = gate(target=target, control=control)
     wfn0 = simulate(c0, initial_state=0)
     wfn1 = simulate(c1, initial_state=0)
     assert(wfn0.isclose(wfn1))
-
-
-
-
diff --git a/tests/test_unary_state_prep.py b/tests/test_unary_state_prep.py
index 0d7ee4ba..a89fded4 100644
--- a/tests/test_unary_state_prep.py
+++ b/tests/test_unary_state_prep.py
@@ -21,7 +21,7 @@ def test_unary_states(target_space: list):
     coeff = 1.0 / numpy.sqrt(qubits)  # fails for the 3-Qubit Case because the wrong sign is picked in the solution
     coeffs = [coeff for i in range(qubits)]
 
-    wfn = QubitWaveFunction()
+    wfn = QubitWaveFunction(qubits)
     for i, c in enumerate(coeffs):
         wfn += c * QubitWaveFunction.from_string("1.0|" + target_space[i] + ">")