code: def random_density(qubits: Union[int, Qubits]) -> Density:
      
      N, qubits = qubits_count_tuple(qubits)
      size = (2**N, 2**N)
      ginibre_ensemble = (np.random.normal(size=size) +
                          1j * np.random.normal(size=size)) / np.sqrt(2.0)
      matrix = ginibre_ensemble @ np.transpose(np.conjugate(ginibre_ensemble))
      matrix /= np.trace(matrix)
  
      return Density(matrix, qubits=qubits)
docs: Returns: A randomly sampled Density from the Hilbert–Schmidt
                  ensemble of quantum states
  
      Ref: "Induced measures in the space of mixed quantum states" Karol
          Zyczkowski, Hans-Juergen Sommers, J. Phys. A34, 7111-7125 (2001)
          https://arxiv.org/abs/quant-ph/0012101