docs: def average(self, rho):
          r"""Return the average density matrix of an inhomogeneous ensemble."""
          def marginal(f, rho):
              remaining = len(f.shape)
              if remaining == 0:
                  return rho
              rho = sum([f[i]*rho[i] for i in range(rho.shape[0])])
              f = np.sum(f, 0)
              return marginal(f, rho)
  
          return marginal(self.distribution, rho)