docs: def inverse(self, rhov, time_derivative=False):
          r"""Fold a vector into a matrix.
  
          The input of this function can be a numpy array or a sympy Matrix.
  
          If the input is understood to represent the time derivative of a
          density matrix, then the flag time_derivative must be set to True.
  
          >>> unfolding = Unfolding(2, real=True, lower_triangular=True,
          ...                       normalized=True)
          >>> rhos = np.array([[0.6, 1+2j], [1-2j, 0.4]])
          >>> print(rhos == unfolding.inverse(unfolding(rhos)))
          [[ True  True]
           [ True  True]]
  
          >>> from fast import define_density_matrix
          >>> from sympy import pprint
          >>> rho = define_density_matrix(2)
          >>> pprint(unfolding.inverse(unfolding(rho)), use_unicode=False)
          [      -rho22 + 1         re(rho21) - I*im(rho21)]
          [                                                ]
          [re(rho21) + I*im(rho21)           rho22         ]
  
  
          >>> rhops = np.array([[0.0, 0.0],
          ...                   [0.0, 0.0]])
  
          >>> print(unfolding.inverse(unfolding(rhops), True))
          [[-0.-0.j  0.-0.j]
           [ 0.+0.j  0.+0.j]]
  
          """
          Ne = self.Ne
          Nrho = self.Nrho
          IJ = self.IJ
  
          if isinstance(rhov, np.ndarray):
              rho = np.zeros((Ne, Ne), complex)
              numeric = True
          elif isinstance(rhov, sympy.Matrix):
              rho = sympy.zeros(Ne, Ne)
              numeric = False
  
          for mu in range(Nrho):
              s, i, j = IJ(mu)
              if numeric:
                  if s == 1:
                      rho[i, j] += rhov[mu]
                  elif s == -1:
                      rho[i, j] += 1j*rhov[mu]
                  elif s == 0:
                      rho[i, j] += rhov[mu]
              else:
                  if s == 1:
                      rho[i, j] += rhov[mu]
                  elif s == -1:
                      rho[i, j] += sympy.I*rhov[mu]
                  elif s == 0:
                      rho[i, j] += rhov[mu]
  
          if self.lower_triangular:
              for i in range(Ne):
                  for j in range(i):
                      rho[j, i] = rho[i, j].conjugate()
  
          if self.normalized:
              if time_derivative:
                  rho[0, 0] = -sum([rho[i, i] for i in range(1, Ne)])
              else:
                  rho[0, 0] = 1-sum([rho[i, i] for i in range(1, Ne)])
  
          return rho