code: def density_matrix_of(self, qubits: List[ops.Qid] = None) -> np.ndarray:
          r
          return density_matrix_from_state_vector(
              self.state_vector(),
              [self.qubit_map[q] for q in qubits] if qubits is not None else None
          )
docs: r"""Returns the density matrix of the state.
  
          Calculate the density matrix for the system on the list, qubits.
          Any qubits not in the list that are present in self.state_vector() will
          be traced out. If qubits is None the full density matrix for
          self.state_vector() is returned, given self.state_vector() follows
          standard Kronecker convention of numpy.kron.
  
          For example:
              self.state_vector() = np.array([1/np.sqrt(2), 1/np.sqrt(2)],
                  dtype=np.complex64)
              qubits = None
              gives us \rho = \begin{bmatrix}
                                  0.5 & 0.5
                                  0.5 & 0.5
                              \end{bmatrix}
  
          Args:
              qubits: list containing qubit IDs that you would like
                  to include in the density matrix (i.e.) qubits that WON'T
                  be traced out.
  
          Returns:
              A numpy array representing the density matrix.
  
          Raises:
              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.