Python numpy.real_if_close() Examples

The following are 30 code examples for showing how to use numpy.real_if_close(). These examples are extracted from open source projects. 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 check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module numpy , or try the search function .

Example 1
Project: pyqmc   Author: WagnerGroup   File: multislater.py    License: 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 2
Project: pyqmc   Author: WagnerGroup   File: multislater.py    License: 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 3
Project: pyGSTi   Author: pyGSTio   File: basistools.py    License: 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 4
Project: vnpy_crypto   Author: birforce   File: fftarma.py    License: 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 5
Project: vnpy_crypto   Author: birforce   File: edgeworth.py    License: 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 6
Project: strawberryfields   Author: XanaduAI   File: test_decompositions_integration.py    License: 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
Project: scqubits   Author: scqubits   File: param_sweep.py    License: 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 8
Project: tick   Author: X-DataInitiative   File: simu_hawkes.py    License: 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 9
Project: thewalrus   Author: XanaduAI   File: test_samples.py    License: 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 10
Project: thewalrus   Author: XanaduAI   File: test_samples.py    License: 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 11
Project: MJHMC   Author: rueberger   File: mixing.py    License: 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 12
Project: Splunking-Crime   Author: nccgroup   File: fftarma.py    License: 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 13
Project: Splunking-Crime   Author: nccgroup   File: edgeworth.py    License: 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 14
Project: forest-benchmarking   Author: rigetti   File: distance_measures.py    License: 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 15
Project: forest-benchmarking   Author: rigetti   File: distance_measures.py    License: 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 16
Project: forest-benchmarking   Author: rigetti   File: distance_measures.py    License: 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 17
Project: forest-benchmarking   Author: rigetti   File: distance_measures.py    License: 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 18
Project: forest-benchmarking   Author: rigetti   File: test_state_tomography.py    License: 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 19
Project: pyqmc   Author: WagnerGroup   File: multislater.py    License: 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() 
Example 20
Project: pyqmc   Author: WagnerGroup   File: multislater.py    License: MIT License 5 votes vote down vote up
def gradient(self, e, epos):
        """ Compute the gradient of the log wave function 
        Note that this can be called even if the internals have not been updated for electron e,
        if epos differs from the current position of electron e."""
        s = int(e >= self._nelec[0])
        aograd = np.real_if_close(
            self._mol.eval_gto("GTOval_sph_deriv1", epos.configs), tol=1e4
        )
        mograd = aograd.dot(self.parameters[self._coefflookup[s]])
        mograd_vals = mograd[:, :, self._det_occup[s]]

        ratios = np.asarray([self._testrow(e, x) for x in mograd_vals])
        return ratios[1:] / ratios[:1] 
Example 21
Project: pyqmc   Author: WagnerGroup   File: multislater.py    License: MIT License 5 votes vote down vote up
def gradient_laplacian(self, e, epos):
        s = int(e >= self._nelec[0])
        ao = np.real_if_close(
            self._mol.eval_gto(self.pbc_str + "GTOval_sph_deriv2", epos.configs)[
                [0, 1, 2, 3, 4, 7, 9]
            ],
            tol=1e4,
        )
        ao = np.concatenate([ao[0:4], ao[4:].sum(axis=0, keepdims=True)])
        mo = np.dot(ao, self.parameters[self._coefflookup[s]])
        mo_vals = mo[:, :, self._det_occup[s]]
        ratios = np.asarray([self._testrow(e, x) for x in mo_vals])
        return ratios[1:-1] / ratios[:1], ratios[-1] / ratios[0] 
Example 22
Project: pyGSTi   Author: pyGSTio   File: termforwardsim.py    License: Apache License 2.0 5 votes vote down vote up
def prs(self, rholabel, elabels, circuit, clipTo, bUseScaling=False, time=None):
        """
        Compute probabilities of a multiple "outcomes" (spam-tuples) for a single
        operation sequence.  The spam tuples may only vary in their effect-label (their
        prep labels must be the same)

        Parameters
        ----------
        rholabel : Label
            The state preparation label.

        elabels : list
            A list of :class:`Label` objects giving the *simplified* effect labels.

        circuit : Circuit or tuple
            A tuple-like object of *simplified* gates (e.g. may include
            instrument elements like 'Imyinst_0')

        clipTo : 2-tuple
          (min,max) to clip returned probability to if not None.
          Only relevant when prMxToFill is not None.

        bUseScaling : bool, optional
          Unused.  Present to match function signature of other calculators.

        time : float, optional
          The *start* time at which `circuit` is evaluated.

        Returns
        -------
        numpy.ndarray
            An array of floating-point probabilities, corresponding to
            the elements of `elabels`.
        """
        assert(time is None), "TermForwardSimulator currently doesn't support time-dependent circuits"
        cpolys = self.prs_as_compact_polys(rholabel, elabels, circuit)
        vals = [_safe_bulk_eval_compact_polys(cpoly[0], cpoly[1], self.paramvec, (1,))[0]
                for cpoly in cpolys]
        ps = _np.array([_np.real_if_close(val) for val in vals])
        if clipTo is not None: ps = _np.clip(ps, clipTo[0], clipTo[1])
        return ps 
Example 23
Project: pyGSTi   Author: pyGSTio   File: lindbladtools.py    License: Apache License 2.0 5 votes vote down vote up
def hamiltonian_to_lindbladian(hamiltonian, sparse=False):
    """
    Construct the Lindbladian corresponding to a given Hamiltonian.

    Mathematically, for a d-dimensional Hamiltonian matrix H, this
    routine constructs the d^2-dimension Lindbladian matrix L whose
    action is given by L(rho) = -1j*2/sqrt(d)*[ H, rho ], where square brackets
    denote the commutator and rho is a density matrix.  L is returned
    as a superoperator matrix that acts on a vectorized density matrices.

    Parameters
    ----------
    hamiltonian : ndarray
      The hamiltonian matrix used to construct the Lindbladian.

    sparse : bool, optional
      Whether to construct a sparse or dense (the default) matrix.

    Returns
    -------
    ndarray or Scipy CSR matrix
    """

    #TODO: there's probably a fast & slick way to so this computation
    #  using vectorization identities
    assert(len(hamiltonian.shape) == 2)
    assert(hamiltonian.shape[0] == hamiltonian.shape[1])
    d = hamiltonian.shape[0]
    if sparse:
        lindbladian = _sps.lil_matrix((d**2, d**2), dtype=hamiltonian.dtype)
    else:
        lindbladian = _np.empty((d**2, d**2), dtype=hamiltonian.dtype)

    for i, rho0 in enumerate(basis_matrices('std', d**2)):  # rho0 == input density mx
        rho1 = _np.sqrt(d) / 2 * (-1j * (_mt.safedot(hamiltonian, rho0) - _mt.safedot(rho0, hamiltonian)))
        lindbladian[:, i] = _np.real_if_close(rho1.flatten()[:, None] if sparse else rho1.flatten())
        # vectorize rho1 & set as linbladian column

    if sparse: lindbladian = lindbladian.tocsr()
    return lindbladian 
Example 24
Project: pyGSTi   Author: pyGSTio   File: rpetools.py    License: Apache License 2.0 5 votes vote down vote up
def extract_theta(model, rpeconfig_inst):
    """
    For a given model, obtain the angle between the estimated "loose axis" and
    the target "loose axis".

    WARNING:  This is a gauge-covariant parameter!  (I think!)  Gauge must be
    fixed prior to estimating.

    Parameters
    ----------
    model : Model
        The model whose loose axis misalignment is to be calculated.

    rpeconfig_inst : rpeconfig object
        Declares which model configuration RPE should be trying to fit;
        determines particular functions and values to be used.

    Returns
    -------
    thetaVal : float
        The value of theta for the input model.
    """
    op_label = rpeconfig_inst.loose_axis_gate_label
    decomp = _decompose_gate_matrix(model.operations[op_label])
    target_axis = rpeconfig_inst.loose_axis_target

    decomp = _decompose_gate_matrix(model.operations[op_label])
    thetaVal = _np.real_if_close([_np.arccos(
        _np.dot(decomp['axis of rotation'], target_axis))])[0]
    if thetaVal > _np.pi / 2:
        thetaVal = _np.pi - thetaVal
    elif thetaVal < -_np.pi / 2:
        thetaVal = _np.pi + thetaVal
    return thetaVal 
Example 25
Project: KAIR   Author: cszn   File: utils_deblur.py    License: MIT License 5 votes vote down vote up
def otf2psf(otf, outsize=None):
    insize = np.array(otf.shape)
    psf = np.fft.ifftn(otf, axes=(0, 1))
    for axis, axis_size in enumerate(insize):
        psf = np.roll(psf, np.floor(axis_size / 2).astype(int), axis=axis)
    if type(outsize) != type(None):
        insize = np.array(otf.shape)
        outsize = np.array(outsize)
        n = max(np.size(outsize), np.size(insize))
        # outsize = postpad(outsize(:), n, 1);
        # insize = postpad(insize(:) , n, 1);
        colvec_out = outsize.flatten().reshape((np.size(outsize), 1))
        colvec_in = insize.flatten().reshape((np.size(insize), 1))
        outsize = np.pad(colvec_out, ((0, max(0, n - np.size(colvec_out))), (0, 0)), mode="constant")
        insize = np.pad(colvec_in, ((0, max(0, n - np.size(colvec_in))), (0, 0)), mode="constant")

        pad = (insize - outsize) / 2
        if np.any(pad < 0):
            print("otf2psf error: OUTSIZE must be smaller than or equal than OTF size")
        prepad = np.floor(pad)
        postpad = np.ceil(pad)
        dims_start = prepad.astype(int)
        dims_end = (insize - postpad).astype(int)
        for i in range(len(dims_start.shape)):
            psf = np.take(psf, range(dims_start[i][0], dims_end[i][0]), axis=i)
    n_ops = np.sum(otf.size * np.log2(otf.shape))
    psf = np.real_if_close(psf, tol=n_ops)
    return psf


# psf2otf copied/modified from https://github.com/aboucaud/pypher/blob/master/pypher/pypher.py 
Example 26
Project: tenpy   Author: tenpy   File: site.py    License: GNU General Public License v3.0 5 votes vote down vote up
def rename_op(self, old_name, new_name):
        """Rename an added operator.

        Parameters
        ----------
        old_name : str
            The old name of the operator.
        new_name : str
            The new name of the operator.
        """
        if old_name == new_name:
            return
        if new_name in self.opnames:
            raise ValueError("new_name already exists")
        old_hc_name = self.hc_ops.get(old_name, None)
        op = getattr(self, old_name)
        need_JW = old_name in self.need_JW_string
        hc_op_name = self.get_hc_op_name(old_name)
        self.remove_op(old_name)
        setattr(self, new_name, op)
        self.opnames.add(new_name)
        if need_JW:
            self.need_JW_string.add(new_name)
        if new_name == 'JW':
            self.JW_exponent = np.real_if_close(np.angle(np.diag(op.to_ndarray())) / np.pi)
        if old_hc_name is not None:
            if old_hc_name == old_name:
                self.hc_ops[new_name] = new_name
            else:
                self.hc_ops[new_name] = old_hc_name
                self.hc_ops[old_hc_name] = new_name 
Example 27
Project: tenpy   Author: tenpy   File: a_mps.py    License: GNU General Public License v3.0 5 votes vote down vote up
def site_expectation_value(self, op):
        """Calculate expectation values of a local operator at each site."""
        result = []
        for i in range(self.L):
            theta = self.get_theta1(i)  # vL i vR
            op_theta = np.tensordot(op, theta, axes=[1, 1])  # i [i*], vL [i] vR
            result.append(np.tensordot(theta.conj(), op_theta, [[0, 1, 2], [1, 0, 2]]))
            # [vL*] [i*] [vR*], [i] [vL] [vR]
        return np.real_if_close(result) 
Example 28
Project: vnpy_crypto   Author: birforce   File: vecm.py    License: MIT License 5 votes vote down vote up
def _estimate_vecm_ml(self):
        y_1_T, delta_y_1_T, y_lag1, delta_x = _endog_matrices(
                self.y, self.exog, self.exog_coint, self.k_ar_diff,
                self.deterministic, self.seasons, self.first_season)
        T = y_1_T.shape[1]

        s00, s01, s10, s11, s11_, _, v = _sij(delta_x, delta_y_1_T, y_lag1)

        beta_tilde = (v[:, :self.coint_rank].T.dot(s11_)).T
        beta_tilde = np.real_if_close(beta_tilde)
        # normalize beta tilde such that eye(r) forms the first r rows of it:
        beta_tilde = np.dot(beta_tilde, inv(beta_tilde[:self.coint_rank]))
        alpha_tilde = s01.dot(beta_tilde).dot(
                inv(beta_tilde.T.dot(s11).dot(beta_tilde)))
        gamma_tilde = (delta_y_1_T - alpha_tilde.dot(beta_tilde.T).dot(y_lag1)
                       ).dot(delta_x.T).dot(inv(np.dot(delta_x, delta_x.T)))
        temp = (delta_y_1_T - alpha_tilde.dot(beta_tilde.T).dot(y_lag1) -
                gamma_tilde.dot(delta_x))
        sigma_u_tilde = temp.dot(temp.T) / T

        return VECMResults(self.y, self.exog, self.exog_coint, self.k_ar,
                           self.coint_rank, alpha_tilde, beta_tilde,
                           gamma_tilde, sigma_u_tilde,
                           deterministic=self.deterministic,
                           seasons=self.seasons, delta_y_1_T=delta_y_1_T,
                           y_lag1=y_lag1, delta_x=delta_x, model=self,
                           names=self.endog_names, dates=self.data.dates,
                           first_season=self.first_season) 
Example 29
Project: vnpy_crypto   Author: birforce   File: fftarma.py    License: MIT License 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 30
Project: garage   Author: rlworkgroup   File: test_bernoulli_mlp_regressor.py    License: MIT License 5 votes vote down vote up
def test_sample_predict(self):
        n_sample = 100
        input_dim = 50
        output_dim = 1
        bmr = BernoulliMLPRegressor(input_shape=(input_dim, ),
                                    output_dim=output_dim)

        xs = np.random.random((input_dim, ))
        p = bmr._f_prob([xs])
        ys = bmr.sample_predict([xs] * n_sample)
        p_predict = np.count_nonzero(ys == 1) / n_sample

        assert np.real_if_close(p, p_predict)