Convert all Qobj to CSR

pulser-core/pulser/noise_model.py

@@ -357,7 +357,7 @@ def _check_eff_noise(
                 operator = np.array(op, dtype=complex)
             except TypeError as e1:
                 try:
-                    operator = np.array(op.data_as("ndarray"), dtype=complex)
+                    operator = np.array(op.to("Dense").data_as("ndarray"), dtype=complex)
                 except AttributeError as e2:
                     raise TypeError(
                         f"Operator {op!r} is not castable to a Numpy array."

pulser-simulation/pulser_simulation/hamiltonian.py

@@ -157,6 +157,7 @@ def _build_collapse_operators(
         if "eff_noise" in config.noise_types:
             for id, rate in enumerate(config.eff_noise_rates):
                 op = np.array(config.eff_noise_opers[id])
+                print(op)
                 basis_dim = len(eigenbasis)
                 op_shape = (basis_dim, basis_dim)
                 if op.shape != op_shape:
@@ -366,10 +367,14 @@ def _build_operator(
                         operator = self.op_matrix[operator]
                     except KeyError:
                         raise ValueError(f"{operator} is not a valid operator")
+                elif isinstance(operator, qutip.Qobj):
+                    operator = operator.to(qutip.core.data.CSR)
+                else:
+                    operator = qutip.Qobj(operator).to(qutip.core.data.CSR)
                 for qubit in qubits:
                     k = self._qid_index[qubit]
                     op_list[k] = operator
-        return qutip.tensor(list(map(qutip.Qobj, op_list)))
+        return qutip.tensor(op_list)
 
     def build_operator(self, operations: Union[list, tuple]) -> qutip.Qobj:
         """Creates an operator with non-trivial actions on some qubits.
@@ -443,8 +448,9 @@ def _get_basis_op_matrices(
     ) -> tuple[dict[States, qutip.Qobj], dict[str, qutip.Qobj]]:
         """Determine basis and projector operators."""
         dim = len(eigenbasis)
-        basis = {b: qutip.basis(dim, i) for i, b in enumerate(eigenbasis)}
-        op_matrix = {"I": qutip.qeye(dim)}
+        with qutip.CoreOptions(default_dtype="CSR"):
+            basis = {b: qutip.basis(dim, i) for i, b in enumerate(eigenbasis)}
+            op_matrix = {"I": qutip.qeye(dim)}
         for proj0 in eigenbasis:
             for proj1 in eigenbasis:
                 proj_name = "sigma_" + proj0 + proj1
@@ -522,7 +528,8 @@ def make_interaction_term(masked: bool = False) -> qutip.Qobj:
                     return 0 * self.build_operator([("I", "global")])
 
             # make interaction term
-            dipole_interaction = cast(qutip.Qobj, 0)
+            with qutip.CoreOptions(default_dtype="CSR"):
+                dipole_interaction = cast(qutip.Qobj, 0)
             for q1, q2 in itertools.combinations(self._qdict.keys(), r=2):
                 if (
                     self._bad_atoms[q1]

pulser-simulation/pulser_simulation/simulation.py

@@ -579,14 +579,22 @@ def _run_solver() -> CoherentResults:
                     self.initial_state,
                     self._eval_times_array,
                     self._hamiltonian._collapse_ops,
-                    options=dict(progress_bar=p_bar, **options),
+                    options=dict(
+                        progress_bar=p_bar,
+                        normalize_output=False,
+                        **options
+                    ),
                 )
             else:
                 result = qutip.sesolve(
                     self._hamiltonian._hamiltonian,
                     self.initial_state,
                     self._eval_times_array,
-                    options=dict(progress_bar=p_bar, **options),
+                    options=dict(
+                        progress_bar=p_bar,
+                        normalize_output=False,
+                        **options
+                    ),
                 )
             results = [
                 QutipResult(

tests/test_qutip_backend.py

@@ -87,3 +87,11 @@ def test_with_default_noise(sequence):
     new_results = backend.run()
     assert isinstance(new_results, NoisyResults)
     assert backend._sim_obj.config == SimConfig.from_noise_model(spam_noise)
+
+proj = [[0, 0], [0, 1]]
+
+@pytest.mark.parametrize("collapse_op", [qutip.sigmax(), qutip.Qobj(proj), np.array(proj), proj])
+def test_collapse_op(sequence, collapse_op):
+    noise_model = pulser.NoiseModel(eff_noise_opers=[collapse_op], eff_noise_rates=[0.1])
+    backend = QutipBackend(sequence, config=pulser.EmulatorConfig(noise_model=noise_model))
+    assert [op.type == qutip.core.data.CSR for op in backend._sim_obj._hamiltonian._collapse_ops]

tests/test_simulation.py

@@ -171,6 +171,7 @@ def test_initialization_and_construction_of_hamiltonian(seq, mod_device):
     )
 
     assert isinstance(sim._hamiltonian._hamiltonian, qutip.QobjEvo)
+    assert sim._hamiltonian._hamiltonian(0).dtype == qutip.core.data.CSR
 
     assert not seq.is_parametrized()
     with pytest.warns(UserWarning, match="returns a copy of itself"):
@@ -812,6 +813,7 @@ def test_noises_rydberg(matrices, noise, result, n_collapse_ops):
             eff_noise_rates=[0.1 if "leakage" in noise else 0.025],
         ),
     )
+    assert [op.type == qutip.core.data.CSR for op in sim._hamiltonian._collapse_ops]
     res = sim.run()
     res_samples = res.sample_final_state()
     assert res_samples == Counter(result)
@@ -832,6 +834,7 @@ def test_relaxation_noise():
 
     sim = QutipEmulator.from_sequence(seq)
     sim.add_config(SimConfig(noise="relaxation", relaxation_rate=0.1))
+    assert [op.type == qutip.core.data.CSR for op in sim._hamiltonian._collapse_ops]
     res = sim.run()
     start_samples = res.sample_state(1)
     ryd_pop = start_samples["1"]
@@ -897,7 +900,7 @@ def test_noises_digital(matrices, noise, result, n_collapse_ops, seq_digital):
             eff_noise_rates=[0.1 if "leakage" in noise else 0.025],
         ),
     )
-
+    assert [op.type == qutip.core.data.CSR for op in sim._hamiltonian._collapse_ops]
     with pytest.raises(
         ValueError,
         match="'relaxation' noise requires addressing of the 'ground-rydberg'",
@@ -935,7 +938,7 @@ def test_noises_digital(matrices, noise, result, n_collapse_ops, seq_digital):
         ("eff_noise", {"111": 958, "110": 19, "011": 12, "101": 11}, 2),
         (
             "relaxation",
-            {"000": 421, "010": 231, "001": 172, "100": 171, "101": 5},
+            {"000": 420, "010": 231, "001": 173, "100": 171, "101": 5},
             1,
         ),
         (("dephasing", "relaxation"), res_deph_relax, 3),
@@ -988,6 +991,7 @@ def test_noises_all(matrices, reg, noise, result, n_collapse_ops, seq):
             eff_noise_rates=[0.2, 0.2],
         ),
     )
+    assert [op.type == qutip.core.data.CSR for op in sim._hamiltonian._collapse_ops]
     with pytest.raises(
         ValueError,
         match="Incompatible shape for effective noise operator n°0.",