docs: def density_matrix_from_state_vector(
      state: Sequence,
      indices: Iterable[int] = None
  ) -> np.ndarray:
      r"""Returns the density matrix of the wavefunction.
  
      Calculate the density matrix for the system on the given qubit
      indices, with the qubits not in indices that are present in state
      traced out. If indices is None the full density matrix for state
      is returned. We assume state follows the standard Kronecker
      convention of numpy.kron.
  
      For example:
  
          state = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=np.complex64)
          indices = None
  
      gives us
  
          $$
          \rho = \begin{bmatrix}
                  0.5 & 0.5
                  0.5 & 0.5
              \end{bmatrix}
          $$
  
      Args:
          state: A sequence representing a wave function in which
              the ordering mapping to qubits follows the standard Kronecker
              convention of numpy.kron.
          indices: list containing indices for qubits that you would like
              to include in the density matrix (i.e.) qubits that WON'T
              be traced out. follows the standard Kronecker convention of
              numpy.kron.
  
      Returns:
          A numpy array representing the density matrix.
  
      Raises:
          ValueError: if the size of state is not a power of 2.
          ValueError: if the size of the state represents more than 25 qubits.
          IndexError: if the indices are out of range for the number of qubits
              corresponding to the state.
      """
      n_qubits = _validate_num_qubits(state)
  
      if indices is None:
          return np.outer(state, np.conj(state))
  
      indices = list(indices)
      _validate_indices(n_qubits, indices)
  
      state = np.asarray(state).reshape((2,)*n_qubits)
  
      sum_inds = np.array(range(n_qubits))
      sum_inds[indices] += n_qubits
  
      rho = np.einsum(state, list(range(n_qubits)), np.conj(state),
          sum_inds.tolist(), indices + sum_inds[indices].tolist())
      new_shape = 2**len(indices)
  
      return rho.reshape((new_shape, new_shape))