Convert a pure state to a density matrix
510881
def asdensity(self) -> 'Density':
        """"""
        matrix = bk.outer(self.tensor, bk.conj(self.tensor))
        return Density(matrix, self.qubits, self._memory)
362371
def state_to_wavefunction(state: State) -> pyquil.Wavefunction:
    """"""
    # TODO: qubits?
    amplitudes = state.vec.asarray()

    # pyQuil labels states backwards.
    amplitudes = amplitudes.transpose()
    amplitudes = amplitudes.reshape([amplitudes.size])
    return pyquil.Wavefunction(amplitudes)
652344
def as_df():
    """"""
    records = []
    for sym, ele in vars(_this).items():
        if sym not in ["Element", "Isotope"] and not sym.startswith("_"):
            for k, v in vars(ele).items():
                if k.startswith("_") and k[1].isdigit():
                    records.append({kk: vv for kk, vv in vars(v).items() if not kk.startswith("_")})
    return _DF.from_records(records)
671310
def dense_output(t_current, t_old, h_current, rcont):
    
    # initialization
    s = (t_current - t_old) / h_current
    s1 = 1.0 - s

    return rcont[0] + s * (rcont[1] + s1 * (
            rcont[2] + s * (rcont[3] + s1 * (rcont[4] + s * (rcont[5] + s1 * (rcont[6] + s * rcont[7]))))))
1032023
def from_dict(d):
        
        return ISA(
            qubits=sorted([Qubit(id=int(qid),
                                 type=q.get("type", DEFAULT_QUBIT_TYPE),
                                 dead=q.get("dead", False))
                           for qid, q in d["1Q"].items()],
                          key=lambda qubit: qubit.id),
            edges=sorted([Edge(targets=[int(q) for q in eid.split('-')],
                               type=e.get("type", DEFAULT_EDGE_TYPE),
                               dead=e.get("dead", False))
                          for eid, e in d["2Q"].items()],
                         key=lambda edge: edge.targets),
        )
1280216
def evolve(self, rho: Density) -> Density:
        """"""
        N = rho.qubit_nb
        qubits = rho.qubits

        indices = list([qubits.index(q) for q in self.qubits]) + \
            list([qubits.index(q) + N for q in self.qubits])

        tensor = bk.tensormul(self.tensor, rho.tensor, indices)
        return Density(tensor, qubits, rho.memory)
902151
def relabel(self, qubits: Qubits) -> 'Density':
        """"""
        return Density(self.vec.tensor, qubits, self._memory)
381519
def join_densities(*densities: Density) -> Density:
    """"""
    vectors = [rho.vec for rho in densities]
    vec = reduce(outer_product, vectors)

    memory = dict(ChainMap(*[rho.memory for rho in densities]))  # TESTME
    return Density(vec.tensor, vec.qubits, memory)
76587
def _one_body_mapping(a_i, a_j, threshold=0.000001):
        
        pauli_list = []
        for alpha in range(2):
            for beta in range(2):
                pauli_prod = Pauli.sgn_prod(a_i[alpha], a_j[beta])
                coeff = 1.0/4 * pauli_prod[1] * np.power(-1j, alpha) * np.power(1j, beta)
                pauli_term = [coeff, pauli_prod[0]]
                if np.absolute(pauli_term[0]) > threshold:
                    pauli_list.append(pauli_term)
        return Operator(paulis=pauli_list)
225628
def do_gate_matrix(self, matrix: np.ndarray,
                       qubits: Sequence[int]) -> 'AbstractQuantumSimulator':
        
        unitary = lifted_gate_matrix(matrix=matrix, qubit_inds=qubits, n_qubits=self.n_qubits)
        self.density = unitary.dot(self.density).dot(np.conj(unitary).T)
        return self
887153
def flatten(self) -> bk.BKTensor:
        """"""
        N = self.qubit_nb
        R = self.rank
        return bk.reshape(self.tensor, [2**N]*R)
169336
def from_pure(cls, z):
        
        return cls(cls._key, {z: 1.0}, {z: 1.0}, pyxray.element_symbol(z))
591657
def ph2full(ptrans, htrans):
    
    n_pstates = len(ptrans)
    n_hstates = len(htrans[0, 0])
    N = n_pstates * n_hstates
    trans = np.zeros((N, N))
    for pidx in range(n_pstates):
        for hidx in range(n_hstates):
            trans[pidx * n_hstates + hidx] = (ptrans[pidx, :, np.newaxis] * htrans[pidx, :, hidx]).flatten()

    return trans
661758
def to_matrix(self, n_qubits=-1):
        """"""
        if n_qubits == -1:
            n_qubits = self.max_n() + 1
        mat = I.to_matrix(n_qubits)
        for op in self.ops:
            if op.is_identity:
                continue
            mat = mat @ op.to_matrix(n_qubits)
        return mat * self.coeff
717532
def _probs(density_matrix: np.ndarray, indices: List[int],
    num_qubits: int) -> List[float]:
    """"""
    # Only diagonal elements matter.
    all_probs = np.diagonal(
        np.reshape(density_matrix, (2 ** num_qubits, 2 ** num_qubits)))
    # Shape into a tensor
    tensor = np.reshape(all_probs, [2] * num_qubits)

    # Calculate the probabilities for measuring the particular results.
    probs = [
        np.sum(np.abs(tensor[linalg.slice_for_qubits_equal_to(indices, b)]))
        for b in range(2 ** len(indices))]

    # To deal with rounding issues, ensure that the probabilities sum to 1.
    probs /= np.sum(probs) # type: ignore
    return probs
189069
def fidelity(rho0: Density, rho1: Density) -> float:
    
    assert rho0.qubit_nb == rho1.qubit_nb   # FIXME

    rho1 = rho1.permute(rho0.qubits)
    op0 = asarray(rho0.asoperator())
    op1 = asarray(rho1.asoperator())

    fid = np.real((np.trace(sqrtm(sqrtm(op0) @ op1 @ sqrtm(op0)))) ** 2)
    fid = min(fid, 1.0)
    fid = max(fid, 0.0)     # Correct for rounding errors

    return fid
909636
def add_densities(density1, density2):
    
    return {spin: np.array(density1[spin]) + np.array(density2[spin])
            for spin in density1.keys()}