content: """
  Classes that describe environments of open quantum systems
  """
  
  # Required for Sphinx to follow autodoc_type_aliases
  from __future__ import annotations
  
  __all__ = ['BosonicEnvironment',
             'DrudeLorentzEnvironment',
             'UnderDampedEnvironment',
             'OhmicEnvironment',
             'ExponentialBosonicEnvironment',
             'FermionicEnvironment',
             'LorentzianEnvironment',
             'ExponentialFermionicEnvironment',
             'CFExponent',
             'system_terminator']
  
  import abc
  import enum
  from time import time
  from typing import Any, Callable, Literal, Sequence, overload, Union
  import warnings
  
  import numpy as np
  from numpy.typing import ArrayLike
  from scipy.linalg import eigvalsh
  from scipy.interpolate import CubicSpline
  
  try:
      from mpmath import mp
      _mpmath_available = True
  except ModuleNotFoundError:
      _mpmath_available = False
  
  from ..utilities import (n_thermal, iterated_fit, aaa,
                           prony_methods, espira1, espira2)
  from .superoperator import spre, spost
  from .qobj import Qobj
  
  
  class BosonicEnvironment(abc.ABC):
      """
      The bosonic environment of an open quantum system. It is characterized by
      its spectral density and temperature or, equivalently, its power spectrum
      or its two-time auto-correlation function.
  
      Use one of the classmethods :meth:`from_spectral_density`,
      :meth:`from_power_spectrum` or :meth:`from_correlation_function` to
      construct an environment manually from one of these characteristic
      functions, or use a predefined sub-class such as the
      :class:`DrudeLorentzEnvironment`, the :class:`UnderDampedEnvironment` or
      the :class:`OhmicEnvironment`.
  
      Bosonic environments offer various ways to approximate the environment with
      a multi-exponential correlation function, which can be used for example in
      the HEOM solver. The approximated environment is represented as a
      :class:`ExponentialBosonicEnvironment`.
  
      All bosonic environments can be approximated by various fitting methods,
      see the description of :meth:`approximate`.
      Subclasses may offer additional approximation methods such as the
      analytical Matsubara or Pade expansions.
  
      Parameters
      ----------
      T : optional, float
          The temperature of this environment.
      tag : optional, str, tuple or any other object
          An identifier (name) for this environment.
      """
  
      def __init__(self, T: float = None, tag: Any = None):
          self.T = T
          self.tag = tag
  
      @abc.abstractmethod
      def spectral_density(self, w: float | ArrayLike) -> (float | ArrayLike):
          """
          The spectral density of this environment. For negative frequencies,
          a value of zero will be returned. See the Users Guide on
          :ref:`bosonic environments <bosonic environments guide>` for specifics
          on the definitions used by QuTiP.
  
          If no analytical expression for the spectral density is known, it will
          be derived from the power spectrum. In this case, the temperature of
          this environment must be specified.
  
          If no analytical expression for the power spectrum is known either, it
          will be derived from the correlation function via a fast fourier
          transform.
  
          Parameters
          ----------
          w : array_like or float
              The frequencies at which to evaluate the spectral density.
          """
  
          ...
  
      @abc.abstractmethod
      def correlation_function(
          self, t: float | ArrayLike, *, eps: float = 1e-10
      ) -> (float | ArrayLike):
          """
          The two-time auto-correlation function of this environment. See the
          Users Guide on :ref:`bosonic environments <bosonic environments guide>`
          for specifics on the definitions used by QuTiP.
  
          If no analytical expression for the correlation function is known, it
          will be derived from the power spectrum via a fast fourier transform.
  
          If no analytical expression for the power spectrum is known either, it
          will be derived from the spectral density. In this case, the
          temperature of this environment must be specified.
  
          Parameters
          ----------
          t : array_like or float
              The times at which to evaluate the correlation function.
  
          eps : optional, float
              Used in case the power spectrum is derived from the spectral
              density; see the documentation of
              :meth:`BosonicEnvironment.power_spectrum`.
          """
  
          ...
  
      @abc.abstractmethod
      def power_spectrum(
          self, w: float | ArrayLike, *, eps: float = 1e-10
      ) -> (float | ArrayLike):
          """
          The power spectrum of this environment. See the Users Guide on
          :ref:`bosonic environments <bosonic environments guide>` for specifics
          on the definitions used by QuTiP.
  
          If no analytical expression for the power spectrum is known, it will
          be derived from the spectral density. In this case, the temperature of
          this environment must be specified.
  
          If no analytical expression for the spectral density is known either,
          the power spectrum will instead be derived from the correlation
          function via a fast fourier transform.
  
          Parameters
          ----------
          w : array_like or float
              The frequencies at which to evaluate the power spectrum.
  
          eps : optional, float
              To derive the zero-frequency power spectrum from the spectral
              density, the spectral density must be differentiated numerically.
              In that case, this parameter is used as the finite difference in
              the numerical differentiation.
          """
  
          ...
  
      # --- user-defined environment creation
  
      @classmethod
      def from_correlation_function(
          cls,
          C: Callable[[float], complex] | ArrayLike,
          tlist: ArrayLike = None,
          tMax: float = None,
          *,
          T: float = None,
          tag: Any = None,
          args: dict[str, Any] = None,
      ) -> BosonicEnvironment:
          r"""
          Constructs a bosonic environment from the provided correlation
          function. The provided function will only be used for times
          :math:`t \geq 0`. At times :math:`t < 0`, the symmetry relation
          :math:`C(-t) = C(t)^\ast` is enforced.
  
          Parameters
          ----------
          C : callable or array_like
              The correlation function. Can be provided as a Python function or
              as an array. When using a function, the signature should be
  
              ``C(t: array_like, **args) -> array_like``
  
              where ``t`` is time and ``args`` is a dict containing the
              other parameters of the function.
  
          tlist : optional, array_like
              The times where the correlation function is sampled (if it is
              provided as an array).
  
          tMax : optional, float
              Specifies that the correlation function is essentially zero outside
              the interval [-tMax, tMax]. Used for numerical integration
              purposes.
  
          T : optional, float
              Environment temperature. (The spectral density of this environment
              can only be calculated from the correlation function if a
              temperature is provided.)
  
          tag : optional, str, tuple or any other object
              An identifier (name) for this environment.
  
          args : optional, dict
              Extra arguments for the correlation function ``C``.
          """
          return _BosonicEnvironment_fromCF(C, tlist, tMax, T, tag, args)
  
      @classmethod
      def from_power_spectrum(
          cls,
          S: Callable[[float], float] | ArrayLike,
          wlist: ArrayLike = None,
          wMax: float = None,
          *,
          T: float = None,
          tag: Any = None,
          args: dict[str, Any] = None,
      ) -> BosonicEnvironment:
          r"""
          Constructs a bosonic environment with the provided power spectrum.
  
          Parameters
          ----------
          S : callable or array_like
              The power spectrum. Can be provided as a Python function or
              as an array. When using a function, the signature should be
  
              ``S(w: array_like, **args) -> array_like``
  
              where ``w`` is the frequency and ``args`` is a dict containing the
              other parameters of the function.
  
          wlist : optional, array_like
              The frequencies where the power spectrum is sampled (if it is
              provided as an array).
  
          wMax : optional, float
              Specifies that the power spectrum is essentially zero outside the
              interval [-wMax, wMax]. Used for numerical integration purposes.
  
          T : optional, float
              Environment temperature. (The spectral density of this environment
              can only be calculated from the powr spectrum if a temperature is
              provided.)
  
          tag : optional, str, tuple or any other object
              An identifier (name) for this environment.
  
          args : optional, dict
              Extra arguments for the power spectrum ``S``.
          """
          return _BosonicEnvironment_fromPS(S, wlist, wMax, T, tag, args)
  
      @classmethod
      def from_spectral_density(
          cls,
          J: Callable[[float], float] | ArrayLike,
          wlist: ArrayLike = None,
          wMax: float = None,
          *,
          T: float = None,
          tag: Any = None,
          args: dict[str, Any] = None,
      ) -> BosonicEnvironment:
          r"""
          Constructs a bosonic environment with the provided spectral density.
          The provided function will only be used for frequencies
          :math:`\omega > 0`. At frequencies :math:`\omega \leq 0`, the spectral
          density is zero according to the definition used by QuTiP. See the
          Users Guide on :ref:`bosonic environments <bosonic environments guide>`
          for a note on spectral densities with support at negative frequencies.
  
          Parameters
          ----------
          J : callable or array_like
              The spectral density. Can be provided as a Python function or
              as an array. When using a function, the signature should be
  
              ``J(w: array_like, **args) -> array_like``
  
              where ``w`` is the frequency and ``args`` is a tuple containing the
              other parameters of the function.
  
          wlist : optional, array_like
              The frequencies where the spectral density is sampled (if it is
              provided as an array).
  
          wMax : optional, float
              Specifies that the spectral density is essentially zero outside the
              interval [-wMax, wMax]. Used for numerical integration purposes.
  
          T : optional, float
              Environment temperature. (The correlation function and the power
              spectrum of this environment can only be calculated from the
              spectral density if a temperature is provided.)
  
          tag : optional, str, tuple or any other object
              An identifier (name) for this environment.
  
          args : optional, dict
              Extra arguments for the spectral density ``J``.
          """
          return _BosonicEnvironment_fromSD(J, wlist, wMax, T, tag, args)
  
      # --- spectral density, power spectrum, correlation function conversions
  
      def _ps_from_sd(self, w, eps, derivative=None):
          # derivative: value of J'(0)
          if self.T is None:
              raise ValueError(
                  "The temperature must be specified for this operation.")
  
          w = np.asarray(w, dtype=float)
          if self.T == 0:
              return 2 * np.heaviside(w, 0) * self.spectral_density(w)
  
          # at zero frequency, we do numerical differentiation
          # S(0) = 2 J'(0) / beta
          zero_mask = (w == 0)
          nonzero_mask = np.invert(zero_mask)
  
          S = np.zeros_like(w)
          if derivative is None:
              S[zero_mask] = 2 * self.T * self.spectral_density(eps) / eps
          else:
              S[zero_mask] = 2 * self.T * derivative
          S[nonzero_mask] = (
              2 * np.sign(w[nonzero_mask])
              * self.spectral_density(np.abs(w[nonzero_mask]))
              * (n_thermal(w[nonzero_mask], self.T) + 1)
          )
          return S.item() if w.ndim == 0 else S
  
      def _sd_from_ps(self, w):
          w = np.asarray(w, dtype=float)
          J = np.zeros_like(w)
          positive_mask = (w > 0)
          power_spectrum = self.power_spectrum(w[positive_mask])
  
          if self.T is None:
              raise ValueError(
                  "The temperature must be specified for this operation.")
  
          J[positive_mask] = (
              power_spectrum / 2 / (n_thermal(w[positive_mask], self.T) + 1)
          )
          return J.item() if w.ndim == 0 else J
  
      def _ps_from_cf(self, w, tMax):
          w = np.asarray(w, dtype=float)
          if w.ndim == 0:
              wMax = np.abs(w)
          elif len(w) == 0:
              return np.array([])
          else:
              wMax = max(np.abs(w[0]), np.abs(w[-1]))
  
          mirrored_result = _fft(self.correlation_function, wMax, tMax=tMax)
          result = np.real(mirrored_result(-w))
          return result.item() if w.ndim == 0 else result
  
      def _cf_from_ps(self, t, wMax, **ps_kwargs):
          t = np.asarray(t, dtype=float)
          if t.ndim == 0:
              tMax = np.abs(t)
          elif len(t) == 0:
              return np.array([])
          else:
              tMax = max(np.abs(t[0]), np.abs(t[-1]))
  
          result_fct = _fft(lambda w: self.power_spectrum(w, **ps_kwargs),
                            tMax, tMax=wMax)
          result = result_fct(t) / (2 * np.pi)
          return result.item() if t.ndim == 0 else result
  
      # --- fitting
  
      @overload
      def approximate(
          self,
          method: Literal['cf'],
          tlist: ArrayLike,
          target_rmse: float = 2e-5,
          Nr_max: int = 10,
          Ni_max: int = 10,
          guess: list[float] = None,
          lower: list[float] = None,
          upper: list[float] = None,
          sigma: float | ArrayLike = None,
          maxfev: int = None,
          full_ansatz: bool = False,
          combine: bool = True,
          tag: Any = None,
      ) -> tuple[ExponentialBosonicEnvironment, dict[str, Any]]:
          ...
  
      @overload
      def approximate(
          self,
          method: Literal['ps'],
          wlist: ArrayLike,
          target_rmse: float = 5e-6,
          Nmax: int = 5,
          guess: list[float] = None,
          lower: list[float] = None,
          upper: list[float] = None,
          sigma: float | ArrayLike = None,
          maxfev: int = None,
          combine: bool = True,
          tag: Any = None,
      ) -> tuple[ExponentialBosonicEnvironment, dict[str, Any]]:
          ...
  
      @overload
      def approximate(
          self,
          method: Literal["sd"],
          wlist: ArrayLike,
          Nk: int = 1,
          target_rmse: float = 5e-6,
          Nmax: int = 10,
          guess: list[float] = None,
          lower: list[float] = None,
          upper: list[float] = None,
          sigma: float | ArrayLike = None,
          maxfev: int = None,
          combine: bool = True,
          tag: Any = None,
      ) -> tuple[ExponentialBosonicEnvironment, dict[str, Any]]:
          ...
  
      @overload
      def approximate(
          self,
          method: Literal['aaa'],
          wlist: ArrayLike,
          tol: float = 1e-13,
          Nmax: int = 10,
          combine: bool = True,
          tag: Any = None
      ) -> tuple[ExponentialBosonicEnvironment, dict[str, Any]]:
          ...
  
      @overload
      def approximate(
          self,
          method: Literal['prony', 'esprit', 'espira-I', 'espira-II'],
          tlist: ArrayLike,
          separate: bool = False,
          Nr: int = 3,
          Ni: int = 3,
          combine: bool = True,
          tag: Any = None
      ) -> tuple[ExponentialBosonicEnvironment, dict[str, Any]]:
          ...
  
      @property
      def _approximation_methods(self):
          return {
              "cf": (self._approx_by_cf_fit, "Correlation Function NLSQ"),
              "ps": (self._approx_by_ps_fit, "Power Spectrum NLSQ"),
              "sd": (self._approx_by_sd_fit, "Spectral Density NLSQ"),
              "aaa": (self._approx_by_aaa, "Power spectrum AAA"),
              "prony": (self._approx_by_prony, "Correlation Function Prony"),
              "esprit": (self._approx_by_prony, "Correlation Function ESPRIT"),
              "espira-i": (self._approx_by_prony,
                           "Correlation Function ESPIRA-I"),
              "espira-ii": (self._approx_by_prony,
                            "Correlation Function ESPIRA-II"),
          }
  
      def approximate(self, method: str, *args, **kwargs):
          """
          Generates a multi-exponential approximation of this environment.
          The available methods are ``"cf"``, ``"ps"``, ``"sd"``, ``"aaa"``,
          ``"prony"``, ``"esprit"``, ``"espira-I"`` and
          ``"espira-II"``. The methods and the parameters required per method are
          documented in the :ref:`Users Guide<environment approximations api>`.
          """
          dispatch = self._approximation_methods
  
          if method.lower() not in dispatch:
              error_string = (f"Unsupported method: {method}."
                              " The available methods are:\n")
              for key in dispatch:
                  error_string += f" - {dispatch[key][1]} ({key})\n"
              error_string += ("If unsure what fitting method to use, you should"
                               " probably use ESPIRA-I, or Power Spectrum NLSQ."
                               " For more information about the fitting methods,"
                               " see the Users Guide.")
              raise ValueError(error_string)
  
          func = dispatch[method.lower()][0]
          return func(method, *args, **kwargs)
  
      def _approx_by_cf_fit(
          self, method, tlist, target_rmse=2e-5, Nr_max=10, Ni_max=10,
          guess=None, lower=None, upper=None, sigma=None, maxfev=None,
          full_ansatz=False, combine=True, tag=None
      ):
  
          # Process arguments
          if tag is None and self.tag is not None:
              tag = (self.tag, f"{method.upper()} Fit")
  
          if full_ansatz:
              num_params = 4
          else:
              num_params = 3
  
          if target_rmse is None:
              target_rmse = 0
              Nr_min, Ni_min = Nr_max, Ni_max
          else:
              Nr_min, Ni_min = 1, 1
  
          clist = self.correlation_function(tlist)
          if guess is None and lower is None and upper is None:
              guess_re, lower_re, upper_re = _default_guess_cfreal(
                  tlist, np.real(clist), full_ansatz)
              guess_im, lower_im, upper_im = _default_guess_cfimag(
                  np.imag(clist), full_ansatz)
          else:
              guess_re, lower_re, upper_re = guess, lower, upper
              guess_im, lower_im, upper_im = guess, lower, upper
  
          # Fit real part
          start_real = time()
          rmse_real, params_real = iterated_fit(
              _cf_real_fit_model, num_params, tlist, np.real(clist), target_rmse,
              Nr_min, Nr_max, guess=guess_re, lower=lower_re, upper=upper_re,
              sigma=sigma, maxfev=maxfev
          )
          end_real = time()
          fit_time_real = end_real - start_real
  
          # Fit imaginary part
          start_imag = time()
          rmse_imag, params_imag = iterated_fit(
              _cf_imag_fit_model, num_params, tlist, np.imag(clist), target_rmse,
              Ni_min, Ni_max, guess=guess_im, lower=lower_im, upper=upper_im,
              sigma=sigma, maxfev=maxfev
          )
          end_imag = time()
          fit_time_imag = end_imag - start_imag
  
          # Generate summary
          Nr = len(params_real)
          Ni = len(params_imag)
          full_summary = _cf_fit_summary(
              params_real, params_imag, fit_time_real, fit_time_imag,
              Nr, Ni, rmse_real, rmse_imag, n=num_params
          )
  
          fit_info = {"Nr": Nr, "Ni": Ni, "fit_time_real": fit_time_real,
                      "fit_time_imag": fit_time_imag, "rmse_real": rmse_real,
                      "rmse_imag": rmse_imag, "params_real": params_real,
                      "params_imag": params_imag, "summary": full_summary}
  
          # Finally, generate environment and return
          ckAR = []
          vkAR = []
          for term in params_real:
              if full_ansatz:
                  a, b, c, d = term
              else:
                  a, b, c = term
                  d = 0
              ckAR.extend([(a + 1j * d) / 2, (a - 1j * d) / 2])
              vkAR.extend([-b - 1j * c, -b + 1j * c])
  
          ckAI = []
          vkAI = []
          for term in params_imag:
              if full_ansatz:
                  a, b, c, d = term
              else:
                  a, b, c = term
                  d = 0
              ckAI.extend([-1j * (a + 1j * d) / 2, 1j * (a - 1j * d) / 2])
              vkAI.extend([-b - 1j * c, -b + 1j * c])
  
          approx_env = ExponentialBosonicEnvironment(
              ckAR, vkAR, ckAI, vkAI, combine=combine, T=self.T, tag=tag)
          return approx_env, fit_info
  
      def _approx_by_sd_fit(
          self, method, wlist, Nk=1, target_rmse=5e-6, Nmax=10,
          guess=None, lower=None, upper=None, sigma=None, maxfev=None,
          combine=True, tag=None
      ):
  
          # Process arguments
          if tag is None and self.tag is not None:
              tag = (self.tag, f"{method.upper()} Fit")
  
          if target_rmse is None:
              target_rmse = 0
              Nmin = Nmax
          else:
              Nmin = 1
  
          jlist = self.spectral_density(wlist)
          if guess is None and lower is None and upper is None:
              guess, lower, upper = _default_guess_sd(wlist, jlist)
  
          # Fit
          start = time()
          rmse, params = iterated_fit(
              _sd_fit_model, 3, wlist, jlist, target_rmse, Nmin, Nmax,
              guess=guess, lower=lower, upper=upper, sigma=sigma, maxfev=maxfev
          )
          end = time()
          fit_time = end - start
  
          # Generate summary
          N = len(params)
          summary = _fit_summary(
              fit_time, rmse, N, "the spectral density", params
          )
          fit_info = {
              "N": N, "Nk": Nk, "fit_time": fit_time, "rmse": rmse,
              "params": params, "summary": summary}
  
          ckAR, vkAR, ckAI, vkAI = [], [], [], []
          # Finally, generate environment and return
          for a, b, c in params:
              lam = np.sqrt(a + 0j)
              gamma = 2 * b
              w0 = np.sqrt(c**2 + b**2)
  
              env = UnderDampedEnvironment(self.T, lam, gamma, w0)
              coeffs = env._matsubara_params(Nk)
              ckAR.extend(coeffs[0])
              vkAR.extend(coeffs[1])
              ckAI.extend(coeffs[2])
              vkAI.extend(coeffs[3])
  
          approx_env = ExponentialBosonicEnvironment(
              ckAR, vkAR, ckAI, vkAI, combine=combine, T=self.T, tag=tag)
          return approx_env, fit_info
  
      def _approx_by_ps_fit(
          self, method, wlist, target_rmse=5e-6, Nmax=5,
          guess=None, lower=None, upper=None, sigma=None, maxfev=None,
          combine=True, tag=None
      ):
  
          # Process arguments
          if tag is None and self.tag is not None:
              tag = (self.tag, f"{method.upper()} Fit")
  
          if target_rmse is None:
              target_rmse = 0
              Nmin = Nmax
          else:
              Nmin = 1
  
          jlist = self.power_spectrum(wlist)
          if guess is None and lower is None and upper is None:
              guess, lower, upper = _default_guess_ps(wlist, jlist)
  
          # Fit
          start = time()
          rmse, params = iterated_fit(
              _ps_fit_model, 4, wlist, jlist, target_rmse, Nmin, Nmax,
              guess=guess, lower=lower, upper=upper, sigma=sigma, maxfev=maxfev
          )
          end = time()
          fit_time = end - start
  
          # Generate summary
          N = len(params)
          summary = _fit_summary(
              fit_time, rmse, N, "the power spectrum", params,
              columns=['a', 'b', 'c', 'd']
          )
          fit_info = {
              "N": N, "fit_time": fit_time, "rmse": rmse,
              "params": params, "summary": summary}
  
          ck, vk = [], []
          # Finally, generate environment and return
          for a, b, c, d in params:
              ck.append(a+1j*b)
              vk.append(c+1j*d)
          ck = np.array(ck)
          vk = np.array(vk)
          ckAR = np.concatenate((ck/2, ck.conj()/2))
          ckAI = np.concatenate((-1j*ck/2, 1j*ck.conj()/2))
          vkAR = np.concatenate((vk, vk.conj()))
          approx_env = ExponentialBosonicEnvironment(
              ckAR, vkAR, ckAI, vkAR, T=self.T, combine=combine, tag=tag)
          return approx_env, fit_info
  
      def _approx_by_aaa(
          self, method, wlist, tol=1e-13, Nmax=10, combine=True, tag=None
      ):
  
          if tag is None and self.tag is not None:
              tag = (self.tag, f"{method.upper()} Fit")
  
          start = time()
          # The *2 is there because half the poles will be filtered out
          result = aaa(self.power_spectrum, wlist,
                       tol=tol,
                       max_iter=Nmax * 2)
          end = time()
          pol = result['poles']
          res = result['residues']
          mask = np.imag(pol) < 0
  
          new_pols, new_res = pol[mask], res[mask]
  
          vk = 1j * new_pols
          ck = -1j * new_res
          # Create complex conjugates for both vk and ck
          ckAR = np.concatenate((ck/2, ck.conj()/2))
          ckAI = np.concatenate((-1j*ck/2, 1j*ck.conj()/2))
          vkAR = np.concatenate((vk, vk.conj()))
  
          cls = ExponentialBosonicEnvironment(
              ck_real=ckAR, vk_real=vkAR, ck_imag=ckAI,
              vk_imag=vkAR, T=self.T, combine=combine, tag=tag)
          # Generate summary
          N = len(vk)
          fit_time = end - start
          params = [(ck.real[i], ck.imag[i], vk[i].real, vk[i].imag)
                    for i in range(len(ck))]
          summary = _fit_summary(
              fit_time, result['rmse'], N, "the power spectrum", params,
              columns=['ckr', 'cki', 'vkr', 'vki']
          )
          fitinfo = {
              "N": N, "fit_time": fit_time, "rmse": result['rmse'],
              "params": params, "summary": summary}
          return cls, fitinfo
  
      def _approx_by_prony(
          self, method, tlist, separate=False, Nr=3, Ni=None,
          combine=True, tag=None
      ):
          if not separate and Ni is not None:
              raise ValueError("The number of imaginary exponents (Ni) cannot be"
                               " specified if real and imaginary parts are fit"
                               " together (separate=False).")
          if Ni is None:
              Ni = 3  # default value
          if tag is None and self.tag is not None:
              tag = (self.tag, f"{method.upper()} Fit")
  
          def prony(x, n):
              return prony_methods(method, x, n)
  
          def phase_to_exponent(phases):
              return -((len(tlist) - 1) / tlist[-1]) * (
                  np.log(np.abs(phases)) + 1j * np.angle(phases)
              )
  
          methods = {"prony": prony, "esprit": prony,
                     "espira-I": espira1, "espira-II": espira2}
  
          if separate:
              start_real = time()
              rmse_real, params_real = methods[method](
                  self.correlation_function(tlist).real, Nr)
              end_real = time()
  
              start_imag = time()
              rmse_imag, params_imag = methods[method](
                  self.correlation_function(tlist).imag, Ni)
              end_imag = time()
  
              ckAR, phases = params_real.T
              ckAI, phases2 = params_imag.T
              vkAR = phase_to_exponent(phases)
              vkAI = phase_to_exponent(phases2)
              cls = ExponentialBosonicEnvironment(
                  ck_real=ckAR, vk_real=vkAR, ck_imag=ckAI,
                  vk_imag=vkAI, T=self.T, combine=combine, tag=tag)
  
              params_real = [(ckAR[i].real, ckAR[i].imag, vkAR[i].real,
                              vkAR[i].imag) for i in range(len(ckAR))]
              params_imag = [(ckAI[i].real, ckAI[i].imag, vkAI[i].real,
                              vkAI[i].imag) for i in range(len(ckAI))]
              fit_time_real = end_real-start_real
              fit_time_imag = end_imag-start_imag
              full_summary = _cf_fit_summary(
                  params_real, params_imag, fit_time_real, fit_time_imag,
                  Nr, Ni, rmse_real, rmse_imag, n=4)
              fit_info = {
                  "Nr": Nr, "Ni": Ni, "fit_time_real": fit_time_real,
                  "fit_time_imag": fit_time_imag, "rmse_real": rmse_real,
                  "rmse_imag": rmse_imag, "params_real": params_real,
                  "params_imag": params_imag, "summary": full_summary,
              }
  
          else:
              start_real = time()
              rmse_real, params_real = methods[method](
                  self.correlation_function(tlist), Nr)
              end_real = time()
  
              amp, phases = params_real.T
              ck = amp
              vk = phase_to_exponent(phases)
              # Create complex conjugates for both vk and ck
              ckAR = np.concatenate((ck/2, ck.conj()/2))
              ckAI = np.concatenate((-1j*ck/2, 1j*ck.conj()/2))
              vkAR = np.concatenate((vk, vk.conj()))
              cls = ExponentialBosonicEnvironment(
                  ck_real=ckAR, vk_real=vkAR, ck_imag=ckAI,
                  vk_imag=vkAR, T=self.T, combine=combine, tag=tag)
  
              params_real = [(ck[i].real, ck[i].imag, vk[i].real,
                              vk[i].imag) for i in range(len(amp))]
              fit_time_real = end_real-start_real
              full_summary = _fit_summary(fit_time_real, rmse_real, Nr,
                                          "Correlation Function", params_real,
                                          columns=['ckr', 'cki', 'vkr', 'vki'])
              fit_info = {
                  "N": Nr, "fit_time": fit_time_real, "rmse": rmse_real,
                  "params": params_real, "summary": full_summary,
              }
  
          return cls, fit_info
  
      def approx_by_cf_fit(self, *args, **kwargs):
          # TODO remove by 5.3
          warnings.warn('The API has changed. Please use approximate("cf", ...)'
                        ' instead of approx_by_cf_fit(...).', FutureWarning)
          return self.approximate("cf", *args, **kwargs)
  
      def approx_by_sd_fit(self, *args, **kwargs):
          # TODO remove by 5.3
          warnings.warn('The API has changed. Please use approximate("sd", ...)'
                        ' instead of approx_by_sd_fit(...).', FutureWarning)
          return self.approximate("sd", *args, **kwargs)
  
  
  class _BosonicEnvironment_fromCF(BosonicEnvironment):
      def __init__(self, C, tlist, tMax, T, tag, args):
          super().__init__(T, tag)
          self._cf = _complex_interpolation(
              C, tlist, 'correlation function', args)
          if tlist is not None:
              self.tMax = max(np.abs(tlist[0]), np.abs(tlist[-1]))
          else:
              self.tMax = tMax
  
      def correlation_function(self, t, **kwargs):
          t = np.asarray(t, dtype=float)
          result = np.zeros_like(t, dtype=complex)
          positive_mask = (t >= 0)
          non_positive_mask = np.invert(positive_mask)
  
          result[positive_mask] = self._cf(t[positive_mask])
          result[non_positive_mask] = np.conj(
              self._cf(-t[non_positive_mask])
          )
          return result.item() if t.ndim == 0 else result
  
      def spectral_density(self, w):
          return self._sd_from_ps(w)
  
      def power_spectrum(self, w, **kwargs):
          if self.tMax is None:
              raise ValueError('The support of the correlation function (tMax) '
                               'must be specified for this operation.')
          return self._ps_from_cf(w, self.tMax)
  
  
  class _BosonicEnvironment_fromPS(BosonicEnvironment):
      def __init__(self, S, wlist, wMax, T, tag, args):
          super().__init__(T, tag)
          self._ps = _real_interpolation(S, wlist, 'power spectrum', args)
          if wlist is not None:
              self.wMax = max(np.abs(wlist[0]), np.abs(wlist[-1]))
          else:
              self.wMax = wMax
  
      def correlation_function(self, t, **kwargs):
          if self.wMax is None:
              raise ValueError('The support of the power spectrum (wMax) '
                               'must be specified for this operation.')
          return self._cf_from_ps(t, self.wMax)
  
      def spectral_density(self, w):
          return self._sd_from_ps(w)
  
      def power_spectrum(self, w, **kwargs):
          w = np.asarray(w, dtype=float)
          ps = self._ps(w)
          return ps.item() if w.ndim == 0 else self._ps(w)
  
  
  class _BosonicEnvironment_fromSD(BosonicEnvironment):
      def __init__(self, J, wlist, wMax, T, tag, args):
          super().__init__(T, tag)
          self._sd = _real_interpolation(J, wlist, 'spectral density', args)
          if wlist is not None:
  