Update HEOM guide to use .e_data instead of .expect to retrieve .e_op…

…s values.

doc/guide/heom/bosonic.rst

@@ -133,7 +133,8 @@ The ``result`` is a standard QuTiP results object with the attributes:
 
 - ``times``: the times at which the state was evaluated (i.e. ``tlist``)
 - ``states``: the system states at each time
-- ``expect``: the values of each ``e_ops`` at each time
+- ``expect``: a list with the values of each ``e_ops`` at each time
+- ``e_data``: a dictionary with the values of each ``e_op`` at each time
 - ``ado_states``: see below (an instance of
   :class:`~qutip.solve.nonmarkov.heom.HierarchyADOsState`)
 
@@ -149,7 +150,7 @@ supply it as the initial state of the solver by calling
 
 As with other QuTiP solvers, if expectation operators or functions are supplied
 using ``.run(..., e_ops=[...])`` the expectation values are available in
-``result.expect``.
+``result.expect`` and ``result.e_data``.
 
 Below we run the solver again, but use ``e_ops`` to store the expectation
 values of the population of the system states and the coherence:
@@ -172,8 +173,8 @@ values of the population of the system states and the coherence:
 
     # Plot the results:
     fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8,8))
-    axes.plot(result.times, result.expect["11"], 'b', linewidth=2, label="P11")
-    axes.plot(result.times, result.expect["12"], 'r', linewidth=2, label="P12")
+    axes.plot(result.times, result.e_data["11"], 'b', linewidth=2, label="P11")
+    axes.plot(result.times, result.e_data["12"], 'r', linewidth=2, label="P12")
     axes.set_xlabel(r't', fontsize=28)
     axes.legend(loc=0, fontsize=12)
 
@@ -398,7 +399,7 @@ occurs:
 
     # Plot populations:
     fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8,8))
-    for label, values in result.expect.items():
+    for label, values in result.e_data.items():
         axes.plot(result.times, values, label=label)
     axes.set_xlabel(r't', fontsize=28)
     axes.set_ylabel(r"Population", fontsize=28)

doc/guide/heom/fermionic.rst

@@ -163,8 +163,8 @@ values of the population of the system states:
 
     # Plot the results:
     fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8,8))
-    axes.plot(result.times, result.expect["11"], 'b', linewidth=2, label="P11")
-    axes.plot(result.times, result.expect["22"], 'r', linewidth=2, label="P22")
+    axes.plot(result.times, result.e_data["11"], 'b', linewidth=2, label="P11")
+    axes.plot(result.times, result.e_data["22"], 'r', linewidth=2, label="P22")
     axes.set_xlabel(r't', fontsize=28)
     axes.legend(loc=0, fontsize=12)
 
@@ -244,11 +244,11 @@ baths, we can pass them to ``e_ops`` and plot the results:
     # Plot the results:
     fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8,8))
     axes.plot(
-        result.times, result.expect["left_currents"], 'b',
+        result.times, result.e_data["left_currents"], 'b',
         linewidth=2, label=r"Bath L",
     )
     axes.plot(
-        result.times, result.expect["right_currents"], 'r',
+        result.times, result.e_data["right_currents"], 'r',
         linewidth=2, label="Bath R",
     )
     axes.set_xlabel(r't', fontsize=28)
@@ -292,15 +292,15 @@ Now we can add the steady state currents to the previous plot:
     # Plot the results and steady state currents:
     fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8,8))
     axes.plot(
-        result.times, result.expect["left_currents"], 'b',
+        result.times, result.e_data["left_currents"], 'b',
         linewidth=2, label=r"Bath L",
     )
     axes.plot(
         result.times, [steady_state_current_left] * len(result.times), 'b:',
         linewidth=2, label=r"Bath L (steady state)",
     )
     axes.plot(
-        result.times, result.expect["right_currents"], 'r',
+        result.times, result.e_data["right_currents"], 'r',
         linewidth=2, label="Bath R",
     )
     axes.plot(