Updated readme with a lof of details

+ Fitting Routine is now adapted to MarkovAnalyzer
+ Updated two examples and removed the numeric_spec example

README.md

@@ -7,7 +7,19 @@ so-called quantum polyspectra. Here we refer to the polyspectra as second to fou
 (integration of the stochastic master equation) is implemented via the QuTiP toolbox whereas the 
 calculation of polyspectra from Hamiltonians or measurements traces recorded in the lab is performed 
 as described in [this paper](https://link.aps.org/doi/10.1103/PhysRevB.98.205143) and [this paper](https://arxiv.org/abs/2011.07992) which also shows the utilization of quantum 
-polyspectra to extract Hamiltonian parameters from a quantum measurement. 
+polyspectra to extract Hamiltonian parameters from a quantum measurement.
+
+The quantum polyspectra approach enable the characterization of very general quantum systems that may include
+* Environmental damping
+* Measurement backaction (Zeno effect) and arbitrary measurement strength
+* Coherent quantum dynamics
+* Stochastic in- and out-tunneling
+* Single-photon measurements
+* Additional detector noise
+* Simultaneous measurement of non-commuting observables
+* Incorporation of temperatures
+* Completely automatic analysis of arbitrary measurement traces
+* Covers all limiting case of weak spin noise measurements, strong measurements resulting in quantum jumps, and single photon sampling
 
 [This poster](Examples/Overview_Poster.pdf) provides an overview of the quantum polyspectra approach to quantum system characterization. 
 Here is a brief summary: The analysis of a continuous measurement record 𝑧(𝑡) poses a fundamental challenge in 
@@ -25,26 +37,15 @@ measurements with stochastically arriving single photons [3,4]. [1] Hägele et a
 [2] Sifft et al., PRR 3, 033123 (2021), [3] Sifft et al. PRA 107, 052203, [4] Sifft et al., arXiv:2310.10464
 
 
-The quantum polyspectra approach enable the characterization of very general quantum systems that may include
-* Environmental damping
-* Measurement backaction (Zeno effect) and arbitrary measurement strength
-* Coherent quantum dynamics
-* Stochastic in- and out-tunneling
-* Single-photon measurements
-* Additional detector noise
-* Simultaneous measurement of non-commuting observables
-* Incorporation of temperatures
-* Completely automatic analysis of arbitrary measurement traces
-* Covers all limiting case of weak spin noise measurements, strong measurements resulting in quantum jumps, and single photon sampling
-
 ## Installation
 QuantumCatch is available on `pip` and can be installed with 
 ```bash
 pip install quantumcatch
 ```
 
 ### Installation of Arrayfire
-Besides running on CPU, QuantumCatch offers GPU support for Nvidia and AMD cards, which are highly recommended for quantum systems with more 
+Besides running on CPU, QuantumCatch offers GPU support for Nvidia and AMD cards. Depending on the hardware used, the
+usage of a GPU is highly recommended for quantum systems with more 
 than about 10 states. A comprehensive installation guide for Linux + NVidia GPU can be found [here](https://github.com/MarkusSifft/QuantumCatch/wiki/Installation-Guide). 
 For GPU calculations the high performance library Arrayfire is used. The Python wrapper ([see here](https://github.com/arrayfire/arrayfire-python)) 
 is automatically installed when installing SignalSnap, however, [ArrayFire C/C++ libraries](https://arrayfire.com/download) need to be installed separately. 
@@ -67,7 +68,7 @@ polyspectra to their measured counterparts.
 ### Example: Continuous Measurement of a Qunatum Dot's Occupation
 Here we are demonstrating how to simulate a simple quantum point 
 contact measurement of the charge state of a quantum dot. First we have to 
-import the QuantumPolyspecrta package. We will also import QuTip and NumPy.
+import the QuantumCatch package. We will also import QuTip and NumPy.
 The analysis and generation module are imported automatically.
 
 ```python
@@ -195,7 +196,7 @@ fig = system.plot()
 ![quantum dot spectra](Examples/plots/quantum_dot_spectra.png)
 
 ## Support
-The development of the QuantumPolyspectra package is supported by the working group Spectroscopy of Condensed Matter 
+The development of the QunatumCatch package is supported by the working group Spectroscopy of Condensed Matter 
 of the Faculty of Physics and Astronomy at the Ruhr University Bochum.
 
 ## Dependencies

src/quantumcatch/__init__.py

@@ -32,5 +32,6 @@
 # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 ###############################################################################
 
-from src.quantumcatch.simulation import *
-from src.quantumcatch.telegraph_noise_tools import *
+from .simulation import *
+from .telegraph_noise_tools import *
+from .fitting_tools import *

src/quantumcatch/fitting_tools.py

@@ -34,7 +34,7 @@
 ###############################################################################
 
 
-from src.quantumcatch.simulation import calc_super_A
+from .simulation import calc_super_A
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.ndimage.filters import gaussian_filter
@@ -55,7 +55,7 @@
 
 class FitSystem:
 
- def __init__(self, set_system, m_op, huber_loss=False, huber_delta=1, enable_gpu=False):
+ def __init__(self, set_system, m_op, f_unit='Hz', huber_loss=False, huber_delta=1, enable_gpu=False):
  self.beta_offset = None
  self.set_system = set_system
  self.m_op = m_op
@@ -70,6 +70,7 @@ def __init__(self, set_system, m_op, huber_loss=False, huber_delta=1, enable_gpu
  self.huber_loss = huber_loss
  self.huber_delta = huber_delta
  self.enable_gpu = enable_gpu
+ self.f_unit = f_unit
 
  def s1(self, system, A):
 
@@ -191,6 +192,7 @@ def complete_fit(self, path, params_in, f_min=None, f_max_2=None, f_max_3=None,
  f"The largest number in fit_orders is {max(fit_orders)}.")
 
  self.measurement_spec = load_spec(path)
+ self.f_unit = self.measurement_spec.config.f_unit
  self.show_plot = show_plot
  self.general_weight = general_weight
  self.beta_offset = beta_offset
@@ -211,7 +213,10 @@ def complete_fit(self, path, params_in, f_min=None, f_max_2=None, f_max_3=None,
  for i in range(1, 5):
  self.f_list[i] = self.measurement_spec.freq[i]
  self.s_list[i] = np.real(self.measurement_spec.S[i])
- self.err_list[i] = np.real(self.measurement_spec.S_err[i])
+ if self.measurement_spec.S_err[i] is None:
+ self.err_list[i] = 1e-6 * np.ones_like(self.s_list[i])
+ else:
+ self.err_list[i] = np.real(self.measurement_spec.S_err[i])
 
  if f_min is not None:
  for i in range(1, 5):
@@ -271,77 +276,32 @@ def complete_fit(self, path, params_in, f_min=None, f_max_2=None, f_max_3=None,
  fit_orders = [1, 2, 3, 4]
  self.fit_orders = fit_orders
 
- print('Low Resolution')
+ self.f_list_original = self.f_list.copy()
+ self.s_list_original = self.s_list.copy()
+ self.err_list_original = self.err_list.copy()
 
- f_list_sampled = [data[::2 ** (i + 3)] for i, data in enumerate(self.f_list)]
+ for res in [100, 50, 20, 10, 5, 2, 1]:
 
- s_list_sampled = []
- for i, data in enumerate(self.s_list):
- if i == 0:
- s_list_sampled.append(data[::2 ** (i + 3)])
- else:
- s_list_sampled.append(data[::2 ** (i + 3), ::2 ** (i + 3)])
+ if res == 1:
+ print('Fitting at full resolution')
 
- err_list_sampled = []
- for i, data in enumerate(self.err_list):
- if i == 0:
- err_list_sampled.append(data[::2 ** (i + 3)])
- else:
- err_list_sampled.append(data[::2 ** (i + 3), ::2 ** (i + 3)])
+ for i in range(1, 5):
+ self.f_list[i] = self.f_list_original[i][::res]
+ if i == 2:
+ self.s_list[i] = self.s_list_original[i][::res]
+ self.err_list[i] = self.err_list_original[i][::res]
+ elif i > 2:
+ self.s_list[i] = self.s_list_original[i][::res, ::res]
+ self.err_list[i] = self.err_list_original[i][::res, ::res]
 
- result = self.start_minimizing(fit_params, method, max_nfev, xtol, ftol)
+  result = self.start_minimizing(fit_params, method, max_nfev, xtol, ftol)
 
- for p in result.params:
- fit_params[p].value = result.params[p].value
+  for p in result.params:
+  fit_params[p].value = result.params[p].value
 
- print('Medium Resolution')
-
- f_list_sampled = [data[::2 ** (i + 2)] for i, data in enumerate(self.f_list)]
-
- s_list_sampled = []
- for i, data in enumerate(self.s_list):
- if i == 0:
- s_list_sampled.append(data[::2 ** (i + 2)])
- else:
- s_list_sampled.append(data[::2 ** (i + 2), ::2 ** (i + 2)])
-
- err_list_sampled = []
- for i, data in enumerate(self.err_list):
- if i == 0:
- err_list_sampled.append(data[::2 ** (i + 2)])
- else:
- err_list_sampled.append(data[::2 ** (i + 2), ::2 ** (i + 2)])
-
- result = self.start_minimizing(fit_params, method, max_nfev, xtol, ftol)
-
- for p in result.params:
- fit_params[p].value = result.params[p].value
-
- print('High Resolution')
-
- f_list_sampled = [data[::2 ** (i + 1)] for i, data in enumerate(self.f_list)]
-
- s_list_sampled = []
- for i, data in enumerate(self.s_list):
- if i == 0:
- s_list_sampled.append(data[::2 ** (i + 1)])
- else:
- s_list_sampled.append(data[::2 ** (i + 1), ::2 ** (i + 1)])
-
- err_list_sampled = []
- for i, data in enumerate(self.err_list):
- if i == 0:
- err_list_sampled.append(data[::2 ** (i + 1)])
- else:
- err_list_sampled.append(data[::2 ** (i + 1), ::2 ** (i + 1)])
-
- result = self.start_minimizing(fit_params, method, max_nfev, xtol, ftol) # TODO .._sampled need to be given
-
- for p in result.params:
- fit_params[p].value = result.params[p].value
-
- print('Full Resolution')
- result = self.start_minimizing(fit_params, method, max_nfev, xtol, ftol)
+ errors = {k: result.params[k].stderr for k in result.params.keys()}
+ self.saved_errors[-1] = errors # Update the last element with the final errors
+ self.display_params(result.params.valuesdict().copy(), self.initial_params, errors)
 
  else:
  print('Parameter fit_order must be: (order_wise, resolution_wise)')
@@ -442,22 +402,6 @@ def comp_plot(self, i, params):
 
  fig, ax = plt.subplots(nrows=3, ncols=3, figsize=(21, 16), gridspec_kw={"width_ratios": [1, 1.2, 1.2]})
 
- # fig = plt.figure(figsize=(20, 20))
- #
- # # Create two separate GridSpec objects: one for the first two rows, and one for the last row
- # gs1 = gridspec.GridSpec(2, 3, figure=fig, width_ratios=[1, 1.2, 1.2])
- # gs2 = gridspec.GridSpec(1, 3, figure=fig, width_ratios=[1, 1, 1])
- #
- # # Place them at the same vertical positions
- # gs1.update(left=0.05, right=0.95, wspace=0.2, hspace=0.2, bottom=0.40)
- # gs2.update(left=0.05, right=0.95, wspace=0.2, hspace=0.2, top=0.35)
- #
- # # Create the list of axes
- # ax_list = [plt.subplot(gs) for gs in [*gs1, *gs2]]
- #
- # # Convert the list into a 3x3 NumPy array for easy indexing
- # ax = np.array(ax_list).reshape(3, 3)
-
  plt.rc('text', usetex=False)
  plt.rc('font', size=10)
  plt.rcParams["axes.axisbelow"] = False
@@ -489,8 +433,8 @@ def comp_plot(self, i, params):
  color=[0, 0.5, 0.9], label='meas.')
  c = ax[0, 1].plot(self.f_list[2], fit_list[2], '--k', alpha=0.8, label='fit')
 
- ax[0, 1].set_ylabel(r"$S^{(2)}_z$ (kHz$^{-1}$)", fontdict={'fontsize': 15})
- ax[0, 1].set_xlabel(r"$\omega/ 2 \pi$ (kHz)", fontdict={'fontsize': 15})
+ ax[0, 1].set_ylabel(r"$S^{(2)}_z$ (" + self.f_unit + r"$^{-1}$)", fontdict={'fontsize': 15})
+ ax[0, 1].set_xlabel(r"$\omega/ 2 \pi$ (" + self.f_unit + ")", fontdict={'fontsize': 15})
 
  ax[0, 1].tick_params(axis='both', direction='in', labelsize=14)
  ax[0, 1].legend()
@@ -502,12 +446,12 @@ def comp_plot(self, i, params):
  color=[0, 0.5, 0.9], label='rel. err.')
  relative_measurement_error = sigma * self.err_list[2] / self.s_list[2]
  ax[0, 2].fill_between(self.f_list[2], relative_measurement_error,
- -relative_measurement_error, alpha=0.3)
- ax[0, 2].plot(self.f_list[2], relative_measurement_error, 'k', alpha=0.5)
- ax[0, 2].plot(self.f_list[2], -relative_measurement_error, 'k', alpha=0.5)
+ -relative_measurement_error, alpha=0.1)
+ ax[0, 2].plot(self.f_list[2], relative_measurement_error, 'k', alpha=0.1)
+ ax[0, 2].plot(self.f_list[2], -relative_measurement_error, 'k', alpha=0.1)
 
- ax[0, 2].set_ylabel(r"$S^{(2)}_z$ (kHz$^{-1}$)", fontdict={'fontsize': 15})
- ax[0, 2].set_xlabel(r"$\omega/ 2 \pi$ (kHz)", fontdict={'fontsize': 15})
+ ax[0, 2].set_ylabel(r"$S^{(2)}_z$ (" + self.f_unit + r"$^{-1}$)", fontdict={'fontsize': 15})
+ ax[0, 2].set_xlabel(r"$\omega/ 2 \pi$ (" + self.f_unit + ")", fontdict={'fontsize': 15})
  ax[0, 2].tick_params(axis='both', direction='in', labelsize=14)
  ax[0, 2].set_title(r'rel. err. and fit deviation in $S^{(2)}_z(\omega)$')
  ax[0, 2].legend()
@@ -536,8 +480,9 @@ def comp_plot(self, i, params):
  c = ax[j, 0].pcolormesh(x, y, z_both - np.diag(np.diag(z_both) / 2), cmap=cmap, norm=norm,
  zorder=1)
 
- ax[j, 0].set_ylabel("\n $\omega_2/ 2 \pi$ (kHz)", labelpad=0, fontdict={'fontsize': 15})
- ax[j, 0].set_xlabel(r"$\omega_1 / 2 \pi$ (kHz)", fontdict={'fontsize': 15})
+ ax[j, 0].set_ylabel("\n $\omega_2/ 2 \pi$ (" + self.f_unit + ")", labelpad=0,
+ fontdict={'fontsize': 15})
+ ax[j, 0].set_xlabel(r"$\omega_1 / 2 \pi$ (" + self.f_unit + ")", fontdict={'fontsize': 15})
 
  ax[j, 0].tick_params(axis='both', direction='in', labelsize=14)
  ax[j, 0].set_title('Fit / Measurement')
@@ -556,7 +501,7 @@ def comp_plot(self, i, params):
 
  err_matrix = np.zeros_like(relative_fit_err)
  relative_measurement_error = sigma * self.err_list[i] / self.s_list[i]
- err_matrix[np.abs(relative_fit_err) < relative_measurement_error] = 1
+ err_matrix[np.abs(relative_fit_err) < np.abs(relative_measurement_error)] = 1
 
  relative_fit_err[relative_fit_err > 0.5] = 0 * relative_fit_err[relative_fit_err > 0.5] + 0.5
  relative_fit_err[relative_fit_err < -0.5] = 0 * relative_fit_err[relative_fit_err < -0.5] - 0.5
@@ -569,8 +514,9 @@ def comp_plot(self, i, params):
  c = ax[j, 1].pcolormesh(x, y, relative_fit_err, cmap=cmap, norm=norm, zorder=1)
  ax[j, 1].pcolormesh(x, y, err_matrix, cmap=cmap_sigma, vmin=0, vmax=1, shading='auto')
 
- ax[j, 1].set_ylabel("\n $\omega_2/ 2 \pi$ (kHz)", labelpad=0, fontdict={'fontsize': 15})
- ax[j, 1].set_xlabel(r"$\omega_1 / 2 \pi$ (kHz)", fontdict={'fontsize': 15})
+ ax[j, 1].set_ylabel("\n $\omega_2/ 2 \pi$ (" + self.f_unit + ")", labelpad=0,
+ fontdict={'fontsize': 15})
+ ax[j, 1].set_xlabel(r"$\omega_1 / 2 \pi$ (" + self.f_unit + ")", fontdict={'fontsize': 15})
 
  ax[j, 1].tick_params(axis='both', direction='in', labelsize=14)
  ax[j, 1].set_title('relative error')
@@ -660,7 +606,7 @@ def comp_plot(self, i, params):
  ax[j, 2].plot(self.f_list[i], s_err_diag_n, 'k', alpha=0.1)
 
  ax[j, 2].set_ylabel(r"arcsinh scaled values", fontdict={'fontsize': 15})
- ax[j, 2].set_xlabel(r"$\omega_1/ 2 \pi$ (kHz)", fontdict={'fontsize': 15})
+ ax[j, 2].set_xlabel(r"$\omega_1/ 2 \pi$ (" + self.f_unit + ")", fontdict={'fontsize': 15})
 
  ax[j, 2].tick_params(axis='both', direction='in', labelsize=14)
  ax[j, 2].legend()

src/quantumcatch/simulation.py

@@ -996,7 +996,7 @@ def cmtr(a, b):
  return liou
 
 
-@njit(cache=True)
+@njit(cache=False)
 def calc_super_liou(h_, c_ops):
  """
  Calculates the superoperator form of the Liouvillian by checking one basis state of the density matrix
@@ -1875,7 +1875,7 @@ def numeric_spec(self, t_window_in, measure_op, f_max, power, order, max_samples
  title_in=None, with_noise=False, _normalize=None,
  roll=False, plot_simulation=False, backend='opencl'):
  """
- Method for the automated calculation of the polyspectra from the numerical integration of the SME.
+ Old method! Not functioning! Method for the automated calculation of the polyspectra from the numerical integration of the SME.
  Can be used as an alternative to the analytic quantum polyspectra or to estimated measurement time and
  noise levels of the spectra.
 

src/quantumcatch/telegraph_noise_tools.py

@@ -32,7 +32,7 @@
 # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 ###############################################################################
 
-from src.quantumcatch import System
+from .simulation import System
 from .fitting_tools import FitSystem
 import numpy as np
 import matplotlib.pyplot as plt