Repository: qutip/qutip (python)

Issue:
Bug in distributions.py concerning Qobj.data and Qobj.data_as()
### Bug Description

Hi, I think there is an bug in the distributions.py file. The functions `update_psi()` and  `update_rho()` of the class `TwoModeQuadratureCorrelation` both have a `Qobj.data[i,j]` type statement as was the norm before qutip5. I encountered this issue when code from my PC would not run on my MacBook. I could fix it by just replacing the `psi.data[i,0]` with a `psi.data_as()[i,0]` (and analogous with `rho.data`) in the distributions.py file as noted in the change log. 

### Code to Reproduce the Bug

```shell
import qutip as qt
a = qt.tensor(qt.destroy(5), qt.qeye(5))
b = qt.tensor(qt.qeye(5), qt.destroy(5))
z = 0.9
rho=qt.squeezing(a,b,z)*qt.tensor(qt.fock(5,0),qt.fock(5,0))
p=qt.TwoModeQuadratureCorrelation(rho,theta1=0,theta2=0)
p.visualize()
```


### Code Output

```shell
TypeError                                 Traceback (most recent call last)
Cell In[69], line 44
     39 print("P quadrature correlation:", p_corr)
     41 # Ensure the correct dimensions are used by TwoModeQuadratureCorrelation
     42 # Calculating two-mode quadrature correlation (if this is the function you mean)
     43 # This is a placeholder for any specific function you're calling:
---> 44 corr = qt.TwoModeQuadratureCorrelation(two_mode_squeezed_state)
     45 # print("Two-mode quadrature correlation:", corr)
File ~/Library/CloudStorage/OneDrive-ETHZurich/PythonScripts/Qutip/TMS_Kerr_simulations/.venv/lib/python3.12/site-packages/qutip/distributions.py:283, in TwoModeQuadratureCorrelation.__init__(self, state, theta1, theta2, extent, steps)
    280 self.theta1 = theta1
    281 self.theta2 = theta2
--> 283 self.update(state)
File ~/Library/CloudStorage/OneDrive-ETHZurich/PythonScripts/Qutip/TMS_Kerr_simulations/.venv/lib/python3.12/site-packages/qutip/distributions.py:291, in TwoModeQuadratureCorrelation.update(self, state)
    286 """
    287 calculate probability distribution for quadrature measurement
    288 outcomes given a two-mode wavefunction or density matrix
    289 """
    290 if isket(state):
--> 291     self.update_psi(state)
    292 else:
    293     self.update_rho(state)
...
    315         i = state_number_index([N, N], [n1, n2])
--> 316         p += kn1 * kn2 * psi.data[i, 0]
    318 self.data = abs(p) ** 2
TypeError: 'qutip.core.data.dense.Dense' object is not subscriptable
```


### Expected Behaviour

(<Figure size 800x600 with 2 Axes>,
 <Axes: xlabel='$X_1(\\theta_1)$', ylabel='$X_2(\\theta_2)$'>)
![image](https://github.com/user-attachments/assets/4087d55e-5f07-40bd-8b09-8869560817d7)


### Your Environment

```shell
QuTiP Version:      5.0.4
Numpy Version:      2.1.2
Scipy Version:      1.14.1
Cython Version:     None
Matplotlib Version: 3.9.2
Python Version:     3.12.5
Number of CPUs:     11
BLAS Info:          Generic
OPENMP Installed:   False

vs. on PC:

QuTiP Version:      4.7.2
Numpy Version:      1.25.2
Scipy Version:      1.11.1
Cython Version:     None
Matplotlib Version: 3.7.1
Python Version:     3.11.4
Number of CPUs:     24
BLAS Info:          INTEL MKL
OPENMP Installed:   False
```


### Additional Context

_No response_

Pull request:
Resolve incompatibility of TwoModeQuadratureCorrelation class with Qu…
…TiP 5.0.0 and above

Adjust update_psi and update_rho functions of TwoModeQuadratureCorrelation class to be compatible with QuTiP 5.0.0. Change Qobj.data to Qobj.data_as() for both functions.

**Checklist**
Thank you for contributing to QuTiP! Please make sure you have finished the following tasks before opening the PR.

- [x] Include the changelog in a file named: `doc/changes/<PR number>.<type>` 'type' can be one of the following: feature, bugfix, doc, removal, misc, or deprecation (see [here](http://qutip.org/docs/latest/development/contributing.html#changelog-generation) for more information).

Delete this checklist after you have completed all the tasks. If you have not finished them all, you can also open a [Draft Pull Request](https://github.blog/2019-02-14-introducing-draft-pull-requests/) to let the others know this on-going work and keep this checklist in the PR description.

**Description**
Adjust update_psi and update_rho functions of TwoModeQuadratureCorrelation class to be compatible with QuTiP 5.0.0. Change Qobj.data to Qobj.data_as() for both functions.

**Related issues or PRs**
fix #2546

Patch:
diff --git a/doc/changes/2548.bugfix b/doc/changes/2548.bugfix
new file mode 100644
index 000000000..0317cb804
--- /dev/null
+++ b/doc/changes/2548.bugfix
@@ -0,0 +1,1 @@
+Resolve incompatibility of TwoModeQuadratureCorrelation class with QuTiP 5.0.0 and above
\ No newline at end of file
diff --git a/qutip/distributions.py b/qutip/distributions.py
index b291d51c1..498953b92 100644
--- a/qutip/distributions.py
+++ b/qutip/distributions.py
@@ -125,11 +125,11 @@ class Distribution:
             fig, ax = plt.subplots(1, 1, figsize=figsize)
 
         if cmap is None:
-            cmap = mpl.cm.get_cmap('RdBu')
+            cmap = mpl.colormaps['RdBu']
 
-        lim = abs(self.data).max()
+        lim = abs(self.data.real).max()
 
-        cf = ax.contourf(self.xvecs[0], self.xvecs[1], self.data, 100,
+        cf = ax.contourf(self.xvecs[0], self.xvecs[1], self.data.real, 100,
                          norm=mpl.colors.Normalize(-lim, lim),
                          cmap=cmap)
 
@@ -152,12 +152,12 @@ class Distribution:
             ax = Axes3D(fig, azim=-62, elev=25)
 
         if cmap is None:
-            cmap = mpl.cm.get_cmap('RdBu')
+            cmap = mpl.colormaps['RdBu']
 
-        lim = abs(self.data).max()
+        lim = abs(self.data.real).max()
 
         X, Y = np.meshgrid(self.xvecs[0], self.xvecs[1])
-        s = ax.plot_surface(X, Y, self.data,
+        s = ax.plot_surface(X, Y, self.data.real,
                             norm=mpl.colors.Normalize(-lim, lim),
                             rstride=5, cstride=5, cmap=cmap, lw=0.1)
 
@@ -177,7 +177,7 @@ class Distribution:
         if not fig and not ax:
             fig, ax = plt.subplots(1, 1, figsize=figsize)
 
-        p = ax.plot(self.xvecs[0], self.data)
+        p = ax.plot(self.xvecs[0], self.data.real)
 
         if show_xlabel:
             ax.set_xlabel(self.xlabels[0], fontsize=12)
@@ -313,7 +313,7 @@ class TwoModeQuadratureCorrelation(Distribution):
                     sqrt(sqrt(pi) * 2 ** n2 * factorial(n2)) * \
                     exp(-X2 ** 2 / 2.0) * np.polyval(hermite(n2), X2)
                 i = state_number_index([N, N], [n1, n2])
-                p += kn1 * kn2 * psi.data[i, 0]
+                p += kn1 * kn2 * psi.data_as()[i, 0]
 
         self.data = abs(p) ** 2
 
@@ -350,7 +350,7 @@ class TwoModeQuadratureCorrelation(Distribution):
                 for p1 in range(N):
                     for p2 in range(N):
                         j = state_number_index([N, N], [p1, p2])
-                        p += M1[n1, p1] * M2[n2, p2] * rho.data[i, j]
+                        p += M1[n1, p1] * M2[n2, p2] * rho.data_as()[i, j]
 
         self.data = p
 
@@ -381,7 +381,7 @@ class HarmonicOscillatorWaveFunction(Distribution):
                 exp(-self.xvecs[0] ** 2 / 2.0) * \
                 np.polyval(hermite(n), self.xvecs[0])
 
-            self.data += k * psi.data[n, 0]
+            self.data += k * psi.data_as()[n, 0]
 
 
 class HarmonicOscillatorProbabilityFunction(Distribution):
@@ -419,4 +419,4 @@ class HarmonicOscillatorProbabilityFunction(Distribution):
                     exp(-self.xvecs[0] ** 2 / 2.0) * \
                     np.polyval(hermite(n), self.xvecs[0])
 
-                self.data += np.conjugate(k_n) * k_m * rho.data[m, n]
+                self.data += np.conjugate(k_n) * k_m * rho.data_as()[m, n]


Test patch:
diff --git a/qutip/tests/test_visualization.py b/qutip/tests/test_visualization.py
index 636b2a348..ed32021f9 100644
--- a/qutip/tests/test_visualization.py
+++ b/qutip/tests/test_visualization.py
@@ -651,3 +651,69 @@ def test_plot_schmidt_Error():
         fig, ax = qutip.plot_schmidt(state)
     assert str(exc_info.value) == text
     plt.close()
+
+
+
+@pytest.mark.parametrize(
+    "state",  [qutip.basis([2, 2], [0, 0]), qutip.ket2dm(qutip.basis([1, 1], [0, 0]))]
+)
+def test_TwoModeQuadratureCorrelation(state):
+    corr = qutip.TwoModeQuadratureCorrelation(state)
+
+    assert isinstance(corr, qutip.distributions.TwoModeQuadratureCorrelation)
+
+@pytest.mark.parametrize(
+    "state",  [qutip.basis([2, 2], [0, 0]), qutip.ket2dm(qutip.basis([1, 1], [0, 0]))]
+)   
+def test_TwoModeQuadratureCorrelation_plot(state):
+    corr = qutip.TwoModeQuadratureCorrelation(state)
+    
+    fig, ax = corr.visualize()
+    plt.close()
+
+    assert isinstance(fig, mpl.figure.Figure)
+    assert isinstance(ax, mpl.axes.Axes)
+
+
+
+@pytest.mark.parametrize(
+    "state",  [qutip.basis([2, 2], [0, 0]), qutip.ket2dm(qutip.basis([1, 1], [0, 0]))]
+)
+def test_HarmonicOscillatorWaveFunction(state):
+    corr = qutip.HarmonicOscillatorWaveFunction(state)
+
+    assert isinstance(corr, qutip.distributions.HarmonicOscillatorWaveFunction)
+
+@pytest.mark.parametrize(
+    "state",  [qutip.basis([2, 2], [0, 0]), qutip.ket2dm(qutip.basis([1, 1], [0, 0]))]
+)   
+def test_HarmonicOscillatorWaveFunction_plot(state):
+    corr = qutip.HarmonicOscillatorWaveFunction(state)
+    
+    fig, ax = corr.visualize()
+    plt.close()
+
+    assert isinstance(fig, mpl.figure.Figure)
+    assert isinstance(ax, mpl.axes.Axes)
+
+
+
+@pytest.mark.parametrize(
+    "state",  [qutip.basis([2, 2], [0, 0]), qutip.ket2dm(qutip.basis([1, 1], [0, 0]))]
+)
+def test_HarmonicOscillatorProbabilityFunction(state):
+    corr = qutip.HarmonicOscillatorProbabilityFunction(state)
+
+    assert isinstance(corr, qutip.distributions.HarmonicOscillatorProbabilityFunction)
+
+@pytest.mark.parametrize(
+    "state",  [qutip.basis([2, 2], [0, 0]), qutip.ket2dm(qutip.basis([1, 1], [0, 0]))]
+)   
+def test_HarmonicOscillatorProbabilityFunction_plot(state):
+    corr = qutip.HarmonicOscillatorProbabilityFunction(state)
+    
+    fig, ax = corr.visualize()
+    plt.close()
+
+    assert isinstance(fig, mpl.figure.Figure)
+    assert isinstance(ax, mpl.axes.Axes)
\ No newline at end of file