Python numpy.real_if_close() Examples

The following are 30 code examples of numpy.real_if_close(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module numpy , or try the search function .
Example #1
Source File: fftarma.py    From Splunking-Crime with GNU Affero General Public License v3.0 6 votes vote down vote up
def invpowerspd(self, n):
        '''autocovariance from spectral density

        scaling is correct, but n needs to be large for numerical accuracy
        maybe padding with zero in fft would be faster
        without slicing it returns 2-sided autocovariance with fftshift

        >>> ArmaFft([1, -0.5], [1., 0.4], 40).invpowerspd(2**8)[:10]
        array([ 2.08    ,  1.44    ,  0.72    ,  0.36    ,  0.18    ,  0.09    ,
                0.045   ,  0.0225  ,  0.01125 ,  0.005625])
        >>> ArmaFft([1, -0.5], [1., 0.4], 40).acovf(10)
        array([ 2.08    ,  1.44    ,  0.72    ,  0.36    ,  0.18    ,  0.09    ,
                0.045   ,  0.0225  ,  0.01125 ,  0.005625])
        '''
        hw = self.fftarma(n)
        return np.real_if_close(fft.ifft(hw*hw.conj()), tol=200)[:n] 
Example #2
Source File: simu_hawkes.py    From tick with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def spectral_radius(self):
        """Compute the spectral radius of the matrix of l1 norm of Hawkes
        kernels.

        Notes
        -----
        If the spectral radius is greater that 1, the hawkes process is not
        stable
        """

        get_norm = np.vectorize(lambda kernel: kernel.get_norm())
        norms = get_norm(self.kernels)

        # It might happens that eig returns a complex number but with a
        # negligible complex part, in this case we keep only the real part
        spectral_radius = max(eig(norms)[0])
        spectral_radius = np.real_if_close(spectral_radius)
        return spectral_radius 
Example #3
Source File: param_sweep.py    From scqubits with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def _recast_dressed_eigendata(self, dressed_eigendata):
        """
        Parameters
        ----------
        dressed_eigendata: list of tuple(evals, qutip evecs)

        Returns
        -------
        SpectrumData
        """
        evals_count = self.evals_count
        energy_table = np.empty(shape=(self.param_count, evals_count), dtype=np.float_)
        state_table = []  # for dressed states, entries are Qobj
        for j in range(self.param_count):
            energy_table[j] = np.real_if_close(dressed_eigendata[j][0])
            state_table.append(dressed_eigendata[j][1])
        specdata = storage.SpectrumData(energy_table, system_params={}, param_name=self.param_name,
                                        param_vals=self.param_vals, state_table=state_table)
        return specdata 
Example #4
Source File: test_samples.py    From thewalrus with Apache License 2.0 6 votes vote down vote up
def test_displaced_single_mode_state_hafnian(self, sample_func):
        """Test the sampling routines by comparing the photon number frequencies and the exact
        probability distribution of a single mode coherent state
        """
        n_samples = 1000
        n_cut = 6
        sigma = np.identity(2)
        mean = 10 * np.array([0.1, 0.25])

        samples = sample_func(sigma, samples=n_samples, mean=mean, cutoff=n_cut)

        probs = np.real_if_close(
            np.array([density_matrix_element(mean, sigma, [i], [i]) for i in range(n_cut)])
        )
        freq, _ = np.histogram(samples[:, 0], bins=np.arange(0, n_cut + 1))
        rel_freq = freq / n_samples
        assert np.allclose(
            rel_freq, probs, rtol=rel_tol / np.sqrt(n_samples), atol=rel_tol / np.sqrt(n_samples)
        ) 
Example #5
Source File: test_samples.py    From thewalrus with Apache License 2.0 6 votes vote down vote up
def test_displaced_two_mode_state_hafnian(self, sample_func):
        """Test the sampling routines by comparing the photon number frequencies and the exact
        probability distribution of a two mode coherent state
        """
        n_samples = 1000
        n_cut = 6
        sigma = np.identity(4)
        mean = 5 * np.array([0.1, 0.25, 0.1, 0.25])
        samples = sample_func(sigma, samples=n_samples, mean=mean, cutoff=n_cut)
        # samples = hafnian_sample_classical_state(sigma, mean = mean, samples = n_samples)
        probs = np.real_if_close(
            np.array(
                [
                    [density_matrix_element(mean, sigma, [i, j], [i, j]) for i in range(n_cut)]
                    for j in range(n_cut)
                ]
            )
        )
        freq, _, _ = np.histogram2d(samples[:, 1], samples[:, 0], bins=np.arange(0, n_cut + 1))
        rel_freq = freq / n_samples

        assert np.allclose(
            rel_freq, probs, rtol=rel_tol / np.sqrt(n_samples), atol=rel_tol / np.sqrt(n_samples)
        ) 
Example #6
Source File: test_decompositions_integration.py    From strawberryfields with Apache License 2.0 6 votes vote down vote up
def test_graph_embed(self, setup_eng, tol):
        """Test that embedding a traceless adjacency matrix A
        results in the property Amat/A = c J, where c is a real constant,
        and J is the all ones matrix"""
        N = 3
        eng, prog = setup_eng(3)

        with prog.context as q:
            ops.GraphEmbed(A) | q

        state = eng.run(prog).state
        Amat = eng.backend.circuit.Amat()

        # check that the matrix Amat is constructed to be of the form
        # Amat = [[B^\dagger, 0], [0, B]]
        assert np.allclose(Amat[:N, :N], Amat[N:, N:].conj().T, atol=tol)
        assert np.allclose(Amat[:N, N:], np.zeros([N, N]), atol=tol)
        assert np.allclose(Amat[N:, :N], np.zeros([N, N]), atol=tol)

        ratio = np.real_if_close(Amat[N:, N:] / A)
        ratio /= ratio[0, 0]
        assert np.allclose(ratio, np.ones([N, N]), atol=tol) 
Example #7
Source File: edgeworth.py    From vnpy_crypto with MIT License 6 votes vote down vote up
def __init__(self, cum, name='Edgeworth expanded normal', **kwds):
        if len(cum) < 2:
            raise ValueError("At least two cumulants are needed.")
        self._coef, self._mu, self._sigma = self._compute_coefs_pdf(cum)
        self._herm_pdf = HermiteE(self._coef)
        if self._coef.size > 2:
            self._herm_cdf = HermiteE(-self._coef[1:])
        else:
            self._herm_cdf = lambda x: 0.

        # warn if pdf(x) < 0 for some values of x within 4 sigma 
        r = np.real_if_close(self._herm_pdf.roots())
        r = (r - self._mu) / self._sigma
        if r[(np.imag(r) == 0) & (np.abs(r) < 4)].any():
            mesg = 'PDF has zeros at %s ' % r
            warnings.warn(mesg, RuntimeWarning)

        kwds.update({'name': name,
                     'momtype': 0})   # use pdf, not ppf in self.moment()
        super(ExpandedNormal, self).__init__(**kwds) 
Example #8
Source File: fftarma.py    From vnpy_crypto with MIT License 6 votes vote down vote up
def invpowerspd(self, n):
        '''autocovariance from spectral density

        scaling is correct, but n needs to be large for numerical accuracy
        maybe padding with zero in fft would be faster
        without slicing it returns 2-sided autocovariance with fftshift

        >>> ArmaFft([1, -0.5], [1., 0.4], 40).invpowerspd(2**8)[:10]
        array([ 2.08    ,  1.44    ,  0.72    ,  0.36    ,  0.18    ,  0.09    ,
                0.045   ,  0.0225  ,  0.01125 ,  0.005625])
        >>> ArmaFft([1, -0.5], [1., 0.4], 40).acovf(10)
        array([ 2.08    ,  1.44    ,  0.72    ,  0.36    ,  0.18    ,  0.09    ,
                0.045   ,  0.0225  ,  0.01125 ,  0.005625])
        '''
        hw = self.fftarma(n)
        return np.real_if_close(fft.ifft(hw*hw.conj()), tol=200)[:n] 
Example #9
Source File: mixing.py    From MJHMC with GNU General Public License v2.0 6 votes vote down vote up
def rectify_evecs(eigs):
    """
    eigs: output of linalg.eig
    normalizes evecs by L1 norm, truncates small complex components,
    ensures things are positive
    """
    evecs = eigs[1].T
    l1_norm = np.abs(evecs).sum(axis=1)
    norm_evecs = evecs / l1_norm[:, np.newaxis]
    real_evals = [np.around(np.real_if_close(l), decimals=5) for l in eigs[0]]
    real_evecs = []

    for v in norm_evecs:
        real_v = np.real_if_close(v)
        if (real_v < 0).all():
            real_v *= -1
        real_evecs.append(real_v)

    # skip sorting for now: argsort is pain because numpy will typecase to complex arr
    #    desc_idx = np.argsort(real_evals)[::-1]
    #   return real_evals[desc_idx], real_evecs[desc_idx]
    return real_evals, real_evecs 
Example #10
Source File: edgeworth.py    From Splunking-Crime with GNU Affero General Public License v3.0 6 votes vote down vote up
def __init__(self, cum, name='Edgeworth expanded normal', **kwds):
        if len(cum) < 2:
            raise ValueError("At least two cumulants are needed.")
        self._coef, self._mu, self._sigma = self._compute_coefs_pdf(cum)
        self._herm_pdf = HermiteE(self._coef)
        if self._coef.size > 2:
            self._herm_cdf = HermiteE(-self._coef[1:])
        else:
            self._herm_cdf = lambda x: 0.

        # warn if pdf(x) < 0 for some values of x within 4 sigma 
        r = np.real_if_close(self._herm_pdf.roots())
        r = (r - self._mu) / self._sigma
        if r[(np.imag(r) == 0) & (np.abs(r) < 4)].any():
            mesg = 'PDF has zeros at %s ' % r
            warnings.warn(mesg, RuntimeWarning)

        kwds.update({'name': name,
                     'momtype': 0})   # use pdf, not ppf in self.moment()
        super(ExpandedNormal, self).__init__(**kwds) 
Example #11
Source File: basistools.py    From pyGSTi with Apache License 2.0 6 votes vote down vote up
def stdmx_to_vec(m, basis):
    """
    Convert a matrix in the standard basis to
     a vector in the Pauli basis.

    Parameters
    ----------
    m : numpy array
        The matrix, shape 2x2 (1Q) or 4x4 (2Q)

    Returns
    -------
    numpy array
        The vector, length 4 or 16 respectively.
    """

    assert(len(m.shape) == 2 and m.shape[0] == m.shape[1])
    basis = Basis.cast(basis, m.shape[0]**2)
    v = _np.empty((basis.size, 1))
    for i, mx in enumerate(basis.elements):
        if basis.real:
            v[i, 0] = _np.real(_mt.trace(_np.dot(mx, m)))
        else:
            v[i, 0] = _np.real_if_close(_mt.trace(_np.dot(mx, m)))
    return v 
Example #12
Source File: multislater.py    From pyqmc with MIT License 6 votes vote down vote up
def updateinternals(self, e, epos, mask=None):
        """Update any internals given that electron e moved to epos. mask is a Boolean array 
        which allows us to update only certain walkers"""
        # MAY want to vectorize later if it really hangs here, shouldn't!

        s = int(e >= self._nelec[0])
        if mask is None:
            mask = [True] * epos.configs.shape[0]
        eeff = e - s * self._nelec[0]
        ao = np.real_if_close(
            self._mol.eval_gto(self.pbc_str + "GTOval_sph", epos.configs), tol=1e4
        )
        self._aovals[:, e, :] = ao
        mo = ao.dot(self.parameters[self._coefflookup[s]])

        mo_vals = mo[:, self._det_occup[s]]
        det_ratio, self._inverse[s][mask, :, :, :] = sherman_morrison_ms(
            eeff, self._inverse[s][mask, :, :, :], mo_vals[mask, :]
        )

        self._updateval(det_ratio, s, mask) 
Example #13
Source File: multislater.py    From pyqmc with MIT License 6 votes vote down vote up
def laplacian(self, e, epos):
        """ Compute the laplacian Psi/ Psi. """
        s = int(e >= self._nelec[0])
        ao = np.real_if_close(
            self._mol.eval_gto(self.pbc_str + "GTOval_sph_deriv2", epos.configs)[
                [0, 4, 7, 9]
            ],
            tol=1e4,
        )
        molap = np.dot(
            [ao[0], ao[1:].sum(axis=0)], self.parameters[self._coefflookup[s]]
        )
        molap_vals = self._testrow(e, molap[1][:, self._det_occup[s]])
        testvalue = self._testrow(e, molap[0][:, self._det_occup[s]])

        return molap_vals / testvalue 
Example #14
Source File: test_state_tomography.py    From forest-benchmarking with Apache License 2.0 6 votes vote down vote up
def test_variance_bootstrap(two_q_tomo_fixture):
    qubits = [0, 1]
    results, rho_true = two_q_tomo_fixture
    rho_est = iterative_mle_state_estimate(results=results, qubits=qubits, tol=1e-4)
    purity = np.trace(rho_est @ rho_est)
    purity = np.real_if_close(purity)
    assert purity.imag == 0.0

    faster_tomo_est = partial(iterative_mle_state_estimate, epsilon=.0001, beta=.5, tol=1e-3)

    boot_purity, boot_var = estimate_variance(results=results, qubits=qubits,
                                              tomo_estimator=faster_tomo_est,
                                              functional=dm.purity,  n_resamples=50,
                                              project_to_physical=False)

    np.testing.assert_allclose(purity, boot_purity, atol=2 * np.sqrt(boot_var), rtol=0.01) 
Example #15
Source File: distance_measures.py    From forest-benchmarking with Apache License 2.0 6 votes vote down vote up
def hilbert_schmidt_ip(A: np.ndarray, B: np.ndarray, tol: float = 1000) -> float:
    r"""
    Computes the Hilbert-Schmidt (HS) inner product between two operators A and B.

    This inner product is defined as

    .. math::

        HS = (A|B) = Tr[A^\dagger B]

    where :math:`|B) = vec(B)` and :math:`(A|` is the dual vector to :math:`|A)`.

    :param A: Is a dim by dim positive matrix with unit trace.
    :param B: Is a dim by dim positive matrix with unit trace.
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: HS inner product which is a scalar.
    """
    hs_ip = np.trace(np.matmul(np.transpose(np.conj(A)), B))
    return np.ndarray.item(np.real_if_close(hs_ip, tol)) 
Example #16
Source File: distance_measures.py    From forest-benchmarking with Apache License 2.0 6 votes vote down vote up
def purity(rho: np.ndarray, dim_renorm=False, tol: float = 1000) -> float:
    """
    Calculates the purity :math:`P = tr[ρ^2]` of a quantum state ρ.

    As stated above lower value of the purity depends on the dimension of ρ's Hilbert space. For
    some applications this can be undesirable. For this reason we introduce an optional dimensional
    renormalization flag with the following behavior

    If the dimensional renormalization flag is FALSE (default) then  1/dim ≤ P ≤ 1.
    If the dimensional renormalization flag is TRUE then 0 ≤ P ≤ 1.

    where dim is the dimension of ρ's Hilbert space.

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param dim_renorm: Boolean, default False.
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: P the purity of the state.
    """
    p = np.trace(rho @ rho)
    if dim_renorm:
        dim = rho.shape[0]
        p = (dim / (dim - 1.0)) * (p - 1.0 / dim)
    return np.ndarray.item(np.real_if_close(p, tol)) 
Example #17
Source File: distance_measures.py    From forest-benchmarking with Apache License 2.0 6 votes vote down vote up
def fidelity(rho: np.ndarray, sigma: np.ndarray, tol: float = 1000) -> float:
    r"""
    Computes the fidelity :math:`F(\rho, \sigma)` between two quantum states rho and sigma.

    If the states are pure the expression reduces to

    .. math::

        F(|psi>,|phi>) = |<psi|phi>|^2

    The fidelity obeys :math:`0 ≤ F(\rho, \sigma) ≤ 1`, where
    :math:`F(\rho, \sigma)=1 iff \rho = \sigma`.

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param sigma: Is a dim by dim positive matrix with unit trace.
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: Fidelity which is a scalar.
    """
    sqrt_rho = sqrtm_psd(rho)
    fid = (np.trace(sqrtm_psd(sqrt_rho @ sigma @ sqrt_rho))) ** 2
    return np.ndarray.item(np.real_if_close(fid, tol)) 
Example #18
Source File: distance_measures.py    From forest-benchmarking with Apache License 2.0 6 votes vote down vote up
def impurity(rho: np.ndarray, dim_renorm=False, tol: float = 1000) -> float:
    """
    Calculates the impurity (or linear entropy) :math:`L = 1 - tr[ρ^2]` of a quantum state ρ.

    As stated above the lower value of the impurity depends on the dimension of ρ's Hilbert space.
    For some applications this can be undesirable. For this reason we introduce an optional
    dimensional renormalization flag with the following behavior

    If the dimensional renormalization flag is FALSE (default) then  0 ≤ L ≤ 1/dim.
    If the dimensional renormalization flag is TRUE then 0 ≤ L ≤ 1.

    where dim is the dimension of ρ's Hilbert space.

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param dim_renorm: Boolean, default False.
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: L the impurity of the state.
    """
    imp = 1 - np.trace(rho @ rho)
    if dim_renorm:
        dim = rho.shape[0]
        imp = (dim / (dim - 1.0)) * imp
    return np.ndarray.item(np.real_if_close(imp, tol)) 
Example #19
Source File: mixing.py    From MJHMC with GNU General Public License v2.0 5 votes vote down vote up
def sg(sampler):
    """returns the spectral gap
    t: transition matrix
    """
    while True:
        try:
            t = sampler.get_empirical_transition_matrix()
            w,v = eig(t)
            w_ord = np.sort(w)[::-1]
            if np.around(np.real_if_close(w_ord[0]), decimals=5) != 1:
                raise Exception("no eval with value 1")
            return 1 - np.absolute(w_ord[1])
        except RuntimeError:
            sampler.sample(1000) 
Example #20
Source File: fftarma.py    From Splunking-Crime with GNU Affero General Public License v3.0 5 votes vote down vote up
def spdpoly(self, w, nma=50):
        '''spectral density from MA polynomial representation for ARMA process

        References
        ----------
        Cochrane, section 8.3.3
        '''
        mpoly = np.polynomial.Polynomial(self.arma2ma(nma))
        hw = mpoly(np.exp(1j * w))
        spd = np.real_if_close(hw * hw.conj() * 0.5/np.pi)
        return spd, w 
Example #21
Source File: _qvm.py    From pyquil with Apache License 2.0 5 votes vote down vote up
def expectation(
        self, prep_prog: Program, operator_programs: Optional[Iterable[Program]] = None
    ) -> List[float]:
        """
        Calculate the expectation value of operators given a state prepared by
        prep_program.

        :note: If the execution of ``quil_program`` is **non-deterministic**, i.e., if it includes
            measurements and/or noisy quantum gates, then the final wavefunction from which the
            expectation values are computed itself only represents a stochastically generated
            sample. The expectations returned from *different* ``expectation`` calls *will then
            generally be different*.

        To measure the expectation of a PauliSum, you probably want to
        do something like this::

                progs, coefs = hamiltonian.get_programs()
                expect_coeffs = np.array(cxn.expectation(prep_program, operator_programs=progs))
                return np.real_if_close(np.dot(coefs, expect_coeffs))

        :param prep_prog: Quil program for state preparation.
        :param operator_programs: A list of Programs, each specifying an operator whose
            expectation to compute. Default is a list containing only the empty Program.
        :return: Expectation values of the operators.
        """
        # Developer note: This code is for backwards compatibility. It can't be replaced with
        # ForestConnection._expectation because we've turned off the ability to set
        # `needs_compilation` (that usually indicates the user is doing something iffy like
        # using a noise model with this function)

        if isinstance(operator_programs, Program):
            warnings.warn(
                "You have provided a Program rather than a list of Programs. The results from "
                "expectation will be line-wise expectation values of the operator_programs.",
                SyntaxWarning,
            )

        payload = self._expectation_payload(prep_prog, operator_programs)
        assert self.sync_endpoint is not None
        response = post_json(self.session, self.sync_endpoint + "/qvm", payload)
        return cast(List[float], response.json()) 
Example #22
Source File: test_quantity_non_ufuncs.py    From Carnets with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def test_real_if_close(self):
        q = np.array([1+0j, 0+1j, 1+1j, 0+0j]) * u.m
        out = np.real_if_close(q)
        expected = np.real_if_close(q.value) * u.m
        assert np.all(out == expected) 
Example #23
Source File: distance_measures.py    From forest-benchmarking with Apache License 2.0 5 votes vote down vote up
def infidelity(rho: np.ndarray, sigma: np.ndarray, tol: float = 1000) -> float:
    r"""
    Computes the infidelity, :math:`1 - F(\rho, \sigma)`, between two quantum states rho and sigma
    where :math:`F(\rho, \sigma)` is the fidelity.

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param sigma: Is a dim by dim positive matrix with unit trace.
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: Infidelity which is a scalar.
    """
    return 1 - fidelity(rho, sigma, tol) 
Example #24
Source File: _reference.py    From pyquil with Apache License 2.0 5 votes vote down vote up
def sample_bitstrings(self, n_samples: int, tol_factor: float = 1e8) -> np.ndarray:
        """
        Sample bitstrings from the distribution defined by the wavefunction.

        Qubit 0 is at ``out[:, 0]``.

        :param n_samples: The number of bitstrings to sample
        :param tol_factor: Tolerance to set imaginary probabilities to zero, relative to
            machine epsilon.
        :return: An array of shape (n_samples, n_qubits)
        """
        if self.rs is None:
            raise ValueError(
                "You have tried to perform a stochastic operation without setting the "
                "random state of the simulator. Might I suggest using a PyQVM object?"
            )

        # for np.real_if_close the actual tolerance is (machine_eps * tol_factor),
        # where `machine_epsilon = np.finfo(float).eps`. If we use tol_factor = 1e8, then the
        # overall tolerance is \approx 2.2e-8.
        probabilities = np.real_if_close(np.diagonal(self.density), tol=tol_factor)
        # Next set negative probabilities to zero
        probabilities = [0 if p < 0.0 else p for p in probabilities]
        # Ensure they sum to one
        probabilities = probabilities / np.sum(probabilities)
        possible_bitstrings = all_bitstrings(self.n_qubits)
        inds = self.rs.choice(2 ** self.n_qubits, n_samples, p=probabilities)
        bitstrings = possible_bitstrings[inds, :]
        bitstrings = np.flip(bitstrings, axis=1)  # qubit ordering: 0 on the left.
        return bitstrings 
Example #25
Source File: base.py    From hmmlearn with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def get_stationary_distribution(self):
        """Compute the stationary distribution of states.
        """
        # The stationary distribution is proportional to the left-eigenvector
        # associated with the largest eigenvalue (i.e., 1) of the transition
        # matrix.
        _utils.check_is_fitted(self, "transmat_")
        eigvals, eigvecs = np.linalg.eig(self.transmat_.T)
        eigvec = np.real_if_close(eigvecs[:, np.argmax(eigvals)])
        return eigvec / eigvec.sum() 
Example #26
Source File: demography.py    From msprime with GNU General Public License v3.0 5 votes vote down vote up
def _matrix_exponential(A):
    """
    Returns the matrix exponential of A.
    https://en.wikipedia.org/wiki/Matrix_exponential
    Note: this is not a general purpose method and is only intended for use within
    msprime.
    """
    d, Y = np.linalg.eig(A)
    Yinv = np.linalg.pinv(Y)
    D = np.diag(np.exp(d))
    B = np.matmul(Y, np.matmul(D, Yinv))
    return np.real_if_close(B, tol=1000) 
Example #27
Source File: test_samples.py    From thewalrus with Apache License 2.0 5 votes vote down vote up
def test_displaced_two_mode_squeezed_state_hafnian(self):
        """Test the sampling routines by comparing the photon number frequencies and the exact
        probability distribution of a displaced two mode squeezed vacuum state
        """
        n_samples = 1000
        n_cut = 10
        mean_n = 1
        r = np.arcsinh(np.sqrt(mean_n))
        c = np.cosh(2 * r)
        s = np.sinh(2 * r)
        sigma = np.array([[c, s, 0, 0], [s, c, 0, 0], [0, 0, c, -s], [0, 0, -s, c]])
        mean = 2 * np.array([0.1, 0.25, 0.1, 0.25])
        samples = hafnian_sample_state(sigma, samples=n_samples, mean=mean, cutoff=n_cut)

        probs = np.real_if_close(
            np.array(
                [
                    [density_matrix_element(mean, sigma, [i, j], [i, j]) for i in range(n_cut)]
                    for j in range(n_cut)
                ]
            )
        )
        freq, _, _ = np.histogram2d(samples[:, 1], samples[:, 0], bins=np.arange(0, n_cut + 1))
        rel_freq = freq / n_samples

        assert np.allclose(
            rel_freq, probs, rtol=rel_tol / np.sqrt(n_samples), atol=rel_tol / np.sqrt(n_samples)
        ) 
Example #28
Source File: test_quantum.py    From thewalrus with Apache License 2.0 5 votes vote down vote up
def test_fidelity_is_symmetric(num_modes, hbar, pure, block_diag):
    """Test that the fidelity is symmetric and between 0 and 1"""
    cov1 = random_covariance(num_modes, hbar=hbar, pure=pure, block_diag=block_diag)
    means1 = np.sqrt(2 * hbar) * np.random.rand(2 * num_modes)
    cov2 = random_covariance(num_modes, hbar=hbar, pure=pure, block_diag=block_diag)
    means2 = np.sqrt(2 * hbar) * np.random.rand(2 * num_modes)
    f12 = fidelity(means1, cov1, means2, cov2, hbar=hbar)
    f21 = fidelity(means2, cov2, means1, cov1, hbar=hbar)
    assert np.allclose(f12, f21)
    assert 0 <= np.real_if_close(f12) < 1.0 
Example #29
Source File: network_features.py    From uncertainpy with GNU General Public License v3.0 5 votes vote down vote up
def van_rossum_dist(self, simulation_end, spiketrains):
        """
        Calculate van Rossum distance.

        Parameters
        ----------
        simulation_end : float
            The simulation end time.
        neo_spiketrains : list
            A list of Neo spiketrains.

        Returns
        -------
        time : None
        van_rossum_dist : 2D array
            The van Rossum distance.
        """
        if len(spiketrains) == 0:
            return None, None

        van_rossum_dist = elephant.spike_train_dissimilarity.van_rossum_dist(spiketrains)

        # van_rossum_dist returns 0.j imaginary parts in some cases
        van_rossum_dist = np.real_if_close(van_rossum_dist)
        if np.any(np.iscomplex(van_rossum_dist)):
            return None, None

        return None, van_rossum_dist 
Example #30
Source File: multislater.py    From pyqmc with MIT License 5 votes vote down vote up
def recompute(self, configs):
        """This computes the value from scratch. Returns the logarithm of the wave function as
        (phase,logdet). If the wf is real, phase will be +/- 1."""

        nconf, nelec, ndim = configs.configs.shape
        mycoords = configs.configs.reshape((nconf * nelec, ndim))
        ao = np.real_if_close(
            self._mol.eval_gto(self.pbc_str + "GTOval_sph", mycoords).reshape(
                (nconf, nelec, -1)
            ),
            tol=1e4,
        )

        self._aovals = ao
        self._dets = []
        self._inverse = []
        for s in [0, 1]:
            mo = ao[:, self._nelec[0] * s : self._nelec[0] + self._nelec[1] * s, :].dot(
                self.parameters[self._coefflookup[s]]
            )
            mo_vals = np.swapaxes(mo[:, :, self._det_occup[s]], 1, 2)
            self._dets.append(
                np.array(np.linalg.slogdet(mo_vals))
            )  # Spin, (sign, val), nconf, [ndet_up, ndet_dn]
            self._inverse.append(
                np.linalg.inv(mo_vals)
            )  # spin, Nconf, [ndet_up, ndet_dn], nelec, nelec
        return self.value()