Python numpy.imag() Examples

The following are 30 code examples for showing how to use numpy.imag(). 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: pyGSTi   Author: pyGSTio   File: __init__.py    License: Apache License 2.0 6 votes vote down vote up
def safe_bulk_eval_compact_polys(vtape, ctape, paramvec, dest_shape):
    """Typechecking wrapper for :function:`bulk_eval_compact_polys`.

    The underlying method has two implementations: one for real-valued
    `ctape`, and one for complex-valued. This wrapper will dynamically
    dispatch to the appropriate implementation method based on the
    type of `ctape`. If the type of `ctape` is known prior to calling,
    it's slightly faster to call the appropriate implementation method
    directly; if not.
    """
    if _np.iscomplexobj(ctape):
        ret = bulk_eval_compact_polys_complex(vtape, ctape, paramvec, dest_shape)
        im_norm = _np.linalg.norm(_np.imag(ret))
        if im_norm > 1e-6:
            print("WARNING: norm(Im part) = {:g}".format(im_norm))
    else:
        ret = bulk_eval_compact_polys(vtape, ctape, paramvec, dest_shape)
    return _np.real(ret) 
Example 2
Project: pyGSTi   Author: pyGSTio   File: reportableqty.py    License: Apache License 2.0 6 votes vote down vote up
def infidelity_diff(self, constant_value):
        """
        Returns a ReportableQty that is the (element-wise in the vector case)
        difference between `constant_value` and this one given by:

        `1.0 - Re(conjugate(constant_value) * self )`
        """
        # let diff(x) = 1.0 - Re(const.C * x) = 1.0 - (const.re * x.re + const.im * x.im)
        # so d(diff)/dx.re = -const.re, d(diff)/dx.im = -const.im
        # diff(x + dx) = diff(x) + d(diff)/dx * dx
        # diff(x + dx) - diff(x) =  - (const.re * dx.re + const.im * dx.im)
        v = 1.0 - _np.real(_np.conjugate(constant_value) * self.value)
        if self.has_eb():
            eb = abs(_np.real(constant_value) * _np.real(self.errbar)
                     + _np.imag(constant_value) * _np.real(self.errbar))
            return ReportableQty(v, eb, self.nonMarkovianEBs)
        else:
            return ReportableQty(v) 
Example 3
Project: pyGSTi   Author: pyGSTio   File: circuit.py    License: Apache License 2.0 6 votes vote down vote up
def _num_to_rqc_str(num):
    """Convert float to string to be included in RQC quil DEFGATE block
    (as written by _np_to_quil_def_str)."""
    num = _np.complex(_np.real_if_close(num))
    if _np.imag(num) == 0:
        output = str(_np.real(num))
        return output
    else:
        real_part = _np.real(num)
        imag_part = _np.imag(num)
        if imag_part < 0:
            sgn = '-'
            imag_part = imag_part * -1
        elif imag_part > 0:
            sgn = '+'
        else:
            assert False
        return '{}{}{}i'.format(real_part, sgn, imag_part) 
Example 4
Project: pyGSTi   Author: pyGSTio   File: slowreplib.py    License: Apache License 2.0 6 votes vote down vote up
def __str__(self):
        def fmt(x):
            if abs(_np.imag(x)) > 1e-6:
                if abs(_np.real(x)) > 1e-6: return "(%.3f+%.3fj)" % (x.real, x.imag)
                else: return "(%.3fj)" % x.imag
            else: return "%.3f" % x.real

        termstrs = []
        sorted_keys = sorted(list(self.keys()))
        for k in sorted_keys:
            vinds = self._int_to_vinds(k)
            varstr = ""; last_i = None; n = 0
            for i in sorted(vinds):
                if i == last_i: n += 1
                elif last_i is not None:
                    varstr += "x%d%s" % (last_i, ("^%d" % n) if n > 1 else "")
                last_i = i
            if last_i is not None:
                varstr += "x%d%s" % (last_i, ("^%d" % n) if n > 1 else "")
            #print("DB: vinds = ",vinds, " varstr = ",varstr)
            if abs(self[k]) > 1e-4:
                termstrs.append("%s%s" % (fmt(self[k]), varstr))
        if len(termstrs) > 0:
            return " + ".join(termstrs)
        else: return "0" 
Example 5
Project: pyGSTi   Author: pyGSTio   File: oplessmodel.py    License: Apache License 2.0 6 votes vote down vote up
def bulk_fill_probs(self, mxToFill, evalTree, clipTo=None, check=False, comm=None):
        if False and evalTree.cache:  # TEST (disabled)
            cpolys = evalTree.cache
            ps = _safe_bulk_eval_compact_polys(cpolys[0], cpolys[1], self._paramvec, (evalTree.num_final_elements(),))
            assert(_np.linalg.norm(_np.imag(ps)) < 1e-6)
            ps = _np.real(ps)
            if clipTo is not None: ps = _np.clip(ps, clipTo[0], clipTo[1])
            mxToFill[:] = ps
        else:
            for i, c in enumerate(evalTree):
                cache = evalTree.cache[i] if evalTree.cache else None
                probs = self.probs(c, clipTo, cache)
                elInds = _slct.indices(evalTree.element_indices[i]) \
                    if isinstance(evalTree.element_indices[i], slice) else evalTree.element_indices[i]
                for k, outcome in zip(elInds, evalTree.outcomes[i]):
                    mxToFill[k] = probs[outcome] 
Example 6
Project: pyGSTi   Author: pyGSTio   File: matrixtools.py    License: Apache License 2.0 6 votes vote down vote up
def is_hermitian(mx, TOL=1e-9):
    """
    Test whether mx is a hermitian matrix.

    Parameters
    ----------
    mx : numpy array
        Matrix to test.

    TOL : float, optional
        Tolerance on absolute magitude of elements.

    Returns
    -------
    bool
        True if mx is hermitian, otherwise False.
    """
    (m, n) = mx.shape
    for i in range(m):
        if abs(mx[i, i].imag) > TOL: return False
        for j in range(i + 1, n):
            if abs(mx[i, j] - mx[j, i].conjugate()) > TOL: return False
    return True 
Example 7
Project: pyGSTi   Author: pyGSTio   File: matrixtools.py    License: Apache License 2.0 6 votes vote down vote up
def complex_compare(a, b):
    """
    Comparison function for complex numbers that compares real part, then
    imaginary part.

    Parameters
    ----------
    a,b : complex

    Returns
    -------
    -1 if a < b
     0 if a == b
    +1 if a > b
    """
    if a.real < b.real: return -1
    elif a.real > b.real: return 1
    elif a.imag < b.imag: return -1
    elif a.imag > b.imag: return 1
    else: return 0 
Example 8
Project: pyGSTi   Author: pyGSTio   File: matrixtools.py    License: Apache License 2.0 6 votes vote down vote up
def safereal(A, inplace=False, check=False):
    """
    Returns the real-part of `A` correctly when `A` is either a dense array or
    a sparse matrix
    """
    if check:
        assert(safenorm(A, 'imag') < 1e-6), "Check failed: taking real-part of matrix w/nonzero imaginary part"
    if _sps.issparse(A):
        if _sps.isspmatrix_csr(A):
            if inplace:
                ret = _sps.csr_matrix((_np.real(A.data), A.indices, A.indptr), shape=A.shape, dtype='d')
            else:  # copy
                ret = _sps.csr_matrix((_np.real(A.data).copy(), A.indices.copy(),
                                       A.indptr.copy()), shape=A.shape, dtype='d')
            ret.eliminate_zeros()
            return ret
        else:
            raise NotImplementedError("safereal() doesn't work with %s matrices yet" % str(type(A)))
    else:
        return _np.real(A) 
Example 9
Project: pyGSTi   Author: pyGSTio   File: matrixtools.py    License: Apache License 2.0 6 votes vote down vote up
def safeimag(A, inplace=False, check=False):
    """
    Returns the imaginary-part of `A` correctly when `A` is either a dense array
    or a sparse matrix
    """
    if check:
        assert(safenorm(A, 'real') < 1e-6), "Check failed: taking imag-part of matrix w/nonzero real part"
    if _sps.issparse(A):
        if _sps.isspmatrix_csr(A):
            if inplace:
                ret = _sps.csr_matrix((_np.imag(A.data), A.indices, A.indptr), shape=A.shape, dtype='d')
            else:  # copy
                ret = _sps.csr_matrix((_np.imag(A.data).copy(), A.indices.copy(),
                                       A.indptr.copy()), shape=A.shape, dtype='d')
            ret.eliminate_zeros()
            return ret
        else:
            raise NotImplementedError("safereal() doesn't work with %s matrices yet" % str(type(A)))
    else:
        return _np.imag(A) 
Example 10
def SaveYML(w_um, RefInd, filename, references='', comments=''):
    
    header = np.empty(9, dtype=object)
    header[0] = '# this file is part of refractiveindex.info database'
    header[1] = '# refractiveindex.info database is in the public domain'
    header[2] = '# copyright and related rights waived via CC0 1.0'
    header[3] = ''
    header[4] = 'REFERENCES:' + references
    header[5] = 'COMMENTS:' + comments
    header[6] = 'DATA:'
    header[7] = '  - type: tabulated nk'
    header[8] = '    data: |'
    
    export = np.column_stack((w_um, np.real(RefInd), np.imag(RefInd)))
    np.savetxt(filename, export, fmt='%4.2f %#.4g %#.4g', delimiter=' ', header='\n'.join(header), comments='',newline='\n        ')
    return

###############################################################################

## Wavelengths to sample ## 
Example 11
def SaveYML(w_um, RefInd, filename, references='', comments=''):
    
    header = np.empty(9, dtype=object)
    header[0] = '# this file is part of refractiveindex.info database'
    header[1] = '# refractiveindex.info database is in the public domain'
    header[2] = '# copyright and related rights waived via CC0 1.0'
    header[3] = ''
    header[4] = 'REFERENCES:' + references
    header[5] = 'COMMENTS:' + comments
    header[6] = 'DATA:'
    header[7] = '  - type: tabulated nk'
    header[8] = '    data: |'
    
    export = np.column_stack((w_um, np.real(RefInd), np.imag(RefInd)))
    np.savetxt(filename, export, fmt='%4.2f %#.4g %#.3e', delimiter=' ', header='\n'.join(header), comments='',newline='\n        ')
    return

###############################################################################

## Wavelengths to sample ## 
Example 12
def SaveYML(w_um, RefInd, filename, references='', comments=''):
    
    header = np.empty(9, dtype=object)
    header[0] = '# this file is part of refractiveindex.info database'
    header[1] = '# refractiveindex.info database is in the public domain'
    header[2] = '# copyright and related rights waived via CC0 1.0'
    header[3] = ''
    header[4] = 'REFERENCES:' + references
    header[5] = 'COMMENTS:' + comments
    header[6] = 'DATA:'
    header[7] = '  - type: tabulated nk'
    header[8] = '    data: |'
    
    export = np.column_stack((w_um, np.real(RefInd), np.imag(RefInd)))
    np.savetxt(filename, export, fmt='%4.3f %#.4g %#.3e', delimiter=' ', header='\n'.join(header), comments='',newline='\n        ')
    return

###############################################################################

## Wavelengths to sample ## 
Example 13
def SaveYML(w_um, RefInd, filename, references='', comments=''):
    
    header = np.empty(9, dtype=object)
    header[0] = '# this file is part of refractiveindex.info database'
    header[1] = '# refractiveindex.info database is in the public domain'
    header[2] = '# copyright and related rights waived via CC0 1.0'
    header[3] = ''
    header[4] = 'REFERENCES:' + references
    header[5] = 'COMMENTS:' + comments
    header[6] = 'DATA:'
    header[7] = '  - type: tabulated nk'
    header[8] = '    data: |'
    
    export = np.column_stack((w_um, np.real(RefInd), np.imag(RefInd)))
    np.savetxt(filename, export, fmt='%4.2f %#.4g %#.3e', delimiter=' ', header='\n'.join(header), comments='',newline='\n        ')
    return

###############################################################################

## Wavelengths to sample ## 
Example 14
def SaveYML(w_um, RefInd, filename, references='', comments=''):
    
    header = np.empty(9, dtype=object)
    header[0] = '# this file is part of refractiveindex.info database'
    header[1] = '# refractiveindex.info database is in the public domain'
    header[2] = '# copyright and related rights waived via CC0 1.0'
    header[3] = ''
    header[4] = 'REFERENCES:' + references
    header[5] = 'COMMENTS:' + comments
    header[6] = 'DATA:'
    header[7] = '  - type: tabulated nk'
    header[8] = '    data: |'
    
    export = np.column_stack((w_um, np.real(RefInd), np.imag(RefInd)))
    np.savetxt(filename, export, fmt='%4.2f %#.4g %#.4g', delimiter=' ', header='\n'.join(header), comments='',newline='\n        ')
    return

###############################################################################

## Wavelengths to sample ## 
Example 15
Project: ludwig   Author: uber   File: audio_utils.py    License: Apache License 2.0 6 votes vote down vote up
def get_group_delay(raw_data, sampling_rate_in_hz, window_length_in_s,
                    window_shift_in_s, num_fft_points, window_type):
    X_stft_transform = _get_stft(raw_data, sampling_rate_in_hz,
                                 window_length_in_s, window_shift_in_s,
                                 num_fft_points, window_type=window_type)
    Y_stft_transform = _get_stft(raw_data, sampling_rate_in_hz,
                                 window_length_in_s, window_shift_in_s,
                                 num_fft_points, window_type=window_type,
                                 data_transformation='group_delay')
    X_stft_transform_real = np.real(X_stft_transform)
    X_stft_transform_imag = np.imag(X_stft_transform)
    Y_stft_transform_real = np.real(Y_stft_transform)
    Y_stft_transform_imag = np.imag(Y_stft_transform)
    nominator = np.multiply(X_stft_transform_real,
                            Y_stft_transform_real) + np.multiply(
        X_stft_transform_imag, Y_stft_transform_imag)
    denominator = np.square(np.abs(X_stft_transform))
    group_delay = np.divide(nominator, denominator + 1e-10)
    assert not np.isnan(
        group_delay).any(), 'There are NaN values in group delay'
    return np.transpose(group_delay) 
Example 16
Project: tenpy   Author: tenpy   File: misc.py    License: GNU General Public License v3.0 6 votes vote down vote up
def zero_if_close(a, tol=1.e-15):
    """set real and/or imaginary part to 0 if their absolute value is smaller than `tol`.

    Parameters
    ----------
    a : ndarray
        numpy array to be rounded
    tol : float
        the threashold which values to consider as '0'.
    """
    if a.dtype == np.complex128 or a.dtype == np.complex64:
        ar = np.choose(np.abs(a.real) < tol, [a.real, np.zeros(a.shape)])
        ai = np.choose(np.abs(a.imag) < tol, [a.imag, np.zeros(a.shape)])
        return ar + 1j * ai
    else:
        return np.choose(np.abs(a) < tol, [a, np.zeros_like(a)]) 
Example 17
Project: lambda-packs   Author: ryfeus   File: ltisys.py    License: MIT License 6 votes vote down vote up
def _order_complex_poles(poles):
    """
    Check we have complex conjugates pairs and reorder P according to YT, ie
    real_poles, complex_i, conjugate complex_i, ....
    The lexicographic sort on the complex poles is added to help the user to
    compare sets of poles.
    """
    ordered_poles = np.sort(poles[np.isreal(poles)])
    im_poles = []
    for p in np.sort(poles[np.imag(poles) < 0]):
        if np.conj(p) in poles:
            im_poles.extend((p, np.conj(p)))

    ordered_poles = np.hstack((ordered_poles, im_poles))

    if poles.shape[0] != len(ordered_poles):
        raise ValueError("Complex poles must come with their conjugates")
    return ordered_poles 
Example 18
Project: lambda-packs   Author: ryfeus   File: _ode.py    License: MIT License 6 votes vote down vote up
def _wrap_jac(self, t, y, *jac_args):
        # jac is the complex Jacobian computed by the user-defined function.
        jac = self.cjac(*((t, y[::2] + 1j * y[1::2]) + jac_args))

        # jac_tmp is the real version of the complex Jacobian.  Each complex
        # entry in jac, say 2+3j, becomes a 2x2 block of the form
        #     [2 -3]
        #     [3  2]
        jac_tmp = zeros((2 * jac.shape[0], 2 * jac.shape[1]))
        jac_tmp[1::2, 1::2] = jac_tmp[::2, ::2] = real(jac)
        jac_tmp[1::2, ::2] = imag(jac)
        jac_tmp[::2, 1::2] = -jac_tmp[1::2, ::2]

        ml = getattr(self._integrator, 'ml', None)
        mu = getattr(self._integrator, 'mu', None)
        if ml is not None or mu is not None:
            # Jacobian is banded.  The user's Jacobian function has computed
            # the complex Jacobian in packed format.  The corresponding
            # real-valued version has every other column shifted up.
            jac_tmp = _transform_banded_jac(jac_tmp)

        return jac_tmp 
Example 19
Project: nevergrad   Author: facebookresearch   File: photonics.py    License: MIT License 6 votes vote down vote up
def creneau(k0: float, a0: float, pol: float, e1: float, e2: float, a: float, n: int, x0: float) -> tp.Tuple[np.ndarray, np.ndarray]:
    nmod = int(n / 2)
    alpha = np.diag(a0 + 2 * np.pi * np.arange(-nmod, nmod + 1))
    if pol == 0:
        M = alpha * alpha - k0 * k0 * marche(e1, e2, a, n, x0)
        L, E = np.linalg.eig(M)
        L = np.sqrt(-L + 0j)
        L = (1 - 2 * (np.imag(L) < -1e-15)) * L
        P = np.block([[E], [np.matmul(E, np.diag(L))]])
    else:
        U = marche(1 / e1, 1 / e2, a, n, x0)
        T = np.linalg.inv(U)
        M = (
            np.matmul(
                np.matmul(np.matmul(T, alpha), np.linalg.inv(marche(e1, e2, a, n, x0))),
                alpha,
            )
            - k0 * k0 * T
        )
        L, E = np.linalg.eig(M)
        L = np.sqrt(-L + 0j)
        L = (1 - 2 * (np.imag(L) < -1e-15)) * L
        P = np.block([[E], [np.matmul(np.matmul(U, E), np.diag(L))]])
    return P, L 
Example 20
Project: spectrum_painter   Author: polygon   File: radios.py    License: MIT License 5 votes vote down vote up
def _interleave(self, complex_iq):
        # Interleave I and Q
        intlv = np.zeros(2*complex_iq.size, dtype=np.float32)
        intlv[0::2] = np.real(complex_iq)
        intlv[1::2] = np.imag(complex_iq)
        return intlv 
Example 21
Project: PyRadarMet   Author: nguy   File: attenuation.py    License: GNU General Public License v2.0 5 votes vote down vote up
def abs_coeff(D, lam, m):
    """
    Absorption coefficient of a spherical particle. Unitless.

    Doviak and Zrnic (1993), Eqn 3.14a or Battan (1973), Eqn 6.6

    Parameters
    ----------
    D : float or array
        Particle diameter [m]
    lam : float
        Radar wavelength [m]
    m : float
        Complex refractive index [unitless], in the form 7.14 - 2.89j

    Notes
    -----
    An example from Battan (1973) is for water at 0C m=7.14-2.89j for a
       wavelength of 3.21 cm and for ice m=1.78-0.0024j for
       wavelength range from 1-10 cm.
    See Battan (1973) Ch.4 , Tables 4.1 and 4.2 for values from
       Gunn and East (1954).
    Also see Doviak and Zrnic (1993), Fig. 3.3 caption.
    """
    Km = (m**2 - 1) / (m**2 + 2)
    Qa = (np.pi**2 * np.asarray(D)**3 / lam) * np.imag(-1 * Km)
    return Qa 
Example 22
Project: python-control   Author: python-control   File: margins.py    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def phase_crossover_frequencies(sys):
    """Compute frequencies and gains at intersections with real axis
    in Nyquist plot.

    Call as:
        omega, gain = phase_crossover_frequencies()

    Returns
    -------
    omega: 1d array of (non-negative) frequencies where Nyquist plot
    intersects the real axis

    gain: 1d array of corresponding gains

    Examples
    --------
    >>> tf = TransferFunction([1], [1, 2, 3, 4])
    >>> PhaseCrossoverFrequenies(tf)
    (array([ 1.73205081,  0.        ]), array([-0.5 ,  0.25]))
    """

    # Convert to a transfer function
    tf = xferfcn._convert_to_transfer_function(sys)

    # if not siso, fall back to (0,0) element
    #! TODO: should add a check and warning here
    num = tf.num[0][0]
    den = tf.den[0][0]

    # Compute frequencies that we cross over the real axis
    numj = (1.j)**np.arange(len(num)-1,-1,-1)*num
    denj = (-1.j)**np.arange(len(den)-1,-1,-1)*den
    allfreq = np.roots(np.imag(np.polymul(numj,denj)))
    realfreq = np.real(allfreq[np.isreal(allfreq)])
    realposfreq = realfreq[realfreq >= 0.]

    # using real() to avoid rounding errors and results like 1+0j
    # it would be nice to have a vectorized version of self.evalfr here
    gain = np.real(np.asarray([tf._evalfr(f)[0][0] for f in realposfreq]))

    return realposfreq, gain 
Example 23
Project: pyGSTi   Author: pyGSTio   File: polynomial.py    License: Apache License 2.0 5 votes vote down vote up
def __str__(self):
        def fmt(x):
            if abs(_np.imag(x)) > 1e-6:
                if abs(_np.real(x)) > 1e-6: return "(%.3f+%.3fj)" % (x.real, x.imag)
                else: return "(%.3fj)" % x.imag
            else: return "%.3f" % x.real

        termstrs = []
        coeffs = self.coeffs
        sorted_keys = sorted(list(coeffs.keys()))
        for k in sorted_keys:
            varstr = ""; last_i = None; n = 1
            for i in sorted(k):
                if i == last_i: n += 1
                elif last_i is not None:
                    varstr += "x%d%s" % (last_i, ("^%d" % n) if n > 1 else "")
                    n = 1
                last_i = i
            if last_i is not None:
                varstr += "x%d%s" % (last_i, ("^%d" % n) if n > 1 else "")
            #print("DB: k = ",k, " varstr = ",varstr)
            if abs(coeffs[k]) > 1e-4:
                termstrs.append("%s%s" % (fmt(coeffs[k]), varstr))
        if len(termstrs) > 0:
            return " + ".join(termstrs)
        else: return "0" 
Example 24
Project: pyGSTi   Author: pyGSTio   File: polynomial.py    License: Apache License 2.0 5 votes vote down vote up
def __str__(self):
        def fmt(x):
            if abs(_np.imag(x)) > 1e-6:
                if abs(_np.real(x)) > 1e-6: return "(%.3f+%.3fj)" % (x.real, x.imag)
                else: return "(%.3fj)" % x.imag
            else: return "%.3f" % x.real

        termstrs = []
        sorted_keys = sorted(list(self.keys()))
        for k in sorted_keys:
            varstr = ""; last_i = None; n = 1
            for i in sorted(k):
                if i == last_i: n += 1
                elif last_i is not None:
                    varstr += "x%d%s" % (last_i, ("^%d" % n) if n > 1 else "")
                    n = 1
                last_i = i
            if last_i is not None:
                varstr += "x%d%s" % (last_i, ("^%d" % n) if n > 1 else "")
            #print("DB: k = ",k, " varstr = ",varstr)
            if abs(self[k]) > 1e-4:
                termstrs.append("%s%s" % (fmt(self[k]), varstr))

        self.check_fastpoly()
        if len(termstrs) > 0:
            return " + ".join(termstrs)
        else: return "0" 
Example 25
Project: pyGSTi   Author: pyGSTio   File: reportableqty.py    License: Apache License 2.0 5 votes vote down vote up
def imag(self):
        """ Returns a ReportableQty that is the imaginary part of this one."""
        if self.has_eb():
            return ReportableQty(_np.imag(self.value), _np.imag(self.errbar), self.nonMarkovianEBs)
        else:
            return ReportableQty(_np.imag(self.value)) 
Example 26
Project: pyGSTi   Author: pyGSTio   File: reportableqty.py    License: Apache License 2.0 5 votes vote down vote up
def hermitian_to_real(self):
        """
        Returns a ReportableQty that holds the real matrix
        whose upper/lower triangle contains the real/imaginary parts
        of the corresponding off-diagonal matrix elements of the
        *Hermitian* matrix stored in this ReportableQty.

        This is used for display purposes.  If this object doesn't
        contain a Hermitian matrix, `ValueError` is raised.
        """
        if _np.linalg.norm(self.value - _np.conjugate(self.value).T) > 1e-8:
            raise ValueError("Contained value must be Hermitian!")

        def _convert(A):
            ret = _np.empty(A.shape, 'd')
            for i in range(A.shape[0]):
                ret[i, i] = A[i, i].real
                for j in range(i + 1, A.shape[1]):
                    ret[i, j] = A[i, j].real
                    ret[j, i] = A[i, j].imag
            return ret

        v = _convert(self.value)
        if self.has_eb():
            eb = _convert(self.errbar)
            return ReportableQty(v, eb, self.nonMarkovianEBs)
        else:
            return ReportableQty(v) 
Example 27
Project: pyGSTi   Author: pyGSTio   File: spamvec.py    License: Apache License 2.0 5 votes vote down vote up
def to_vector(self):
        """
        Get the SPAM vector parameters as an array of values.

        Returns
        -------
        numpy array
            The parameters as a 1D array with length num_params().
        """
        if self._evotype == "statevec":
            return _np.concatenate((self.base1D.real, self.base1D.imag), axis=0)
        else:
            return self.base1D 
Example 28
Project: pyGSTi   Author: pyGSTio   File: spamvec.py    License: Apache License 2.0 5 votes vote down vote up
def _construct_vector(self):
        dmDim = self.dmDim

        #  params is an array of length dmDim^2-1 that
        #  encodes a lower-triangular matrix "Lmx" via:
        #  Lmx[i,i] = params[i*dmDim + i] / param-norm  # i = 0...dmDim-2
        #     *last diagonal el is given by sqrt(1.0 - sum(L[i,j]**2))
        #  Lmx[i,j] = params[i*dmDim + j] + 1j*params[j*dmDim+i] (i > j)

        param2Sum = _np.vdot(self.params, self.params)  # or "dot" would work, since params are real
        paramNorm = _np.sqrt(param2Sum)  # also the norm of *all* Lmx els

        for i in range(dmDim):
            self.Lmx[i, i] = self.params[i * dmDim + i] / paramNorm
            for j in range(i):
                self.Lmx[i, j] = (self.params[i * dmDim + j] + 1j * self.params[j * dmDim + i]) / paramNorm

        Lmx_norm = _np.trace(_np.dot(self.Lmx.T.conjugate(), self.Lmx))  # sum of magnitude^2 of all els
        assert(_np.isclose(Lmx_norm, 1.0)), "Violated trace=1 condition!"

        #The (complex, Hermitian) density matrix is build by
        # assuming Lmx is its Cholesky decomp, which makes
        # the density matrix is pos-def.
        density_mx = _np.dot(self.Lmx, self.Lmx.T.conjugate())
        assert(_np.isclose(_np.trace(density_mx), 1.0)), "density matrix must be trace == 1"

        # write density matrix in given basis: = sum_i alpha_i B_i
        # ASSUME that basis is orthogonal, i.e. Tr(Bi^dag*Bj) = delta_ij
        basis_mxs = _np.rollaxis(self.basis_mxs, 2)  # shape (dmDim, dmDim, len(vec))
        vec = _np.array([_np.trace(_np.dot(M.T.conjugate(), density_mx)) for M in basis_mxs])

        #for now, assume Liouville vector should always be real (TODO: add 'real' flag later?)
        assert(_np.linalg.norm(_np.imag(vec)) < IMAG_TOL)
        vec = _np.real(vec)

        self.base1D.flags.writeable = True
        self.base1D[:] = vec[:]  # so shape is (dim,1) - the convention for spam vectors
        self.base1D.flags.writeable = False 
Example 29
Project: pyGSTi   Author: pyGSTio   File: oplessmodel.py    License: Apache License 2.0 5 votes vote down vote up
def bulk_fill_dprobs(self, mxToFill, evalTree, prMxToFill=None, clipTo=None,
                         check=False, comm=None, wrtBlockSize=None,
                         profiler=None, gatherMemLimit=None):

        Np = self.num_params()
        p = self.to_vector()

        if False and evalTree.cache:  # TEST (disabled)
            cpolys = evalTree.cache
            if prMxToFill is not None:
                ps = _safe_bulk_eval_compact_polys(cpolys[0], cpolys[1], p, (evalTree.num_final_elements(),))
                assert(_np.linalg.norm(_np.imag(ps)) < 1e-6)
                ps = _np.real(ps)
                if clipTo is not None: ps = _np.clip(ps, clipTo[0], clipTo[1])
                prMxToFill[:] = ps
            dpolys = _compact_deriv(cpolys[0], cpolys[1], list(range(Np)))
            dps = _safe_bulk_eval_compact_polys(dpolys[0], dpolys[1], p, (evalTree.num_final_elements(), Np))
            mxToFill[:, :] = dps
        else:
            # eps = 1e-6
            for i, c in enumerate(evalTree):
                cache = evalTree.cache[i] if evalTree.cache else None
                probs0 = self.probs(c, clipTo, cache)
                dprobs0 = self.dprobs(c, False, clipTo, cache)
                elInds = _slct.indices(evalTree.element_indices[i]) \
                    if isinstance(evalTree.element_indices[i], slice) else evalTree.element_indices[i]
                for k, outcome in zip(elInds, evalTree.outcomes[i]):
                    if prMxToFill is not None:
                        prMxToFill[k] = probs0[outcome]
                    mxToFill[k, :] = dprobs0[outcome]

                    #Do this to fill mxToFill instead of calling dprobs above as it's a little faster for finite diff?
                    #for j in range(Np):
                    #    p_plus_dp = p.copy()
                    #    p_plus_dp[j] += eps
                    #    self.from_vector(p_plus_dp)
                    #    probs1 = self.probs(c,clipTo,cache)
                    #    mxToFill[k,j] = (probs1[outcome]-probs0[outcome]) / eps
                    #self.from_vector(p) 
Example 30
Project: pyGSTi   Author: pyGSTio   File: matrixtools.py    License: Apache License 2.0 5 votes vote down vote up
def mx_to_string_complex(m, real_width=9, im_width=9, prec=4):
    """
    Generate a "pretty-format" string for a complex-valued matrix.

    Parameters
    ----------
    mx : numpy array
        the matrix (2-D array) to convert.

    real_width : int, opitonal
        the width (in characters) of the real part of each element.

    im_width : int, opitonal
        the width (in characters) of the imaginary part of each element.

    prec : int optional
        the precision (in characters) of each element's real and imaginary parts.

    Returns
    -------
    string
        matrix m as a pretty formated string.
    """
    if len(m.shape) == 1: m = m[None, :]  # so it works w/vectors too
    s = ""; tol = 10**(-prec)
    for i in range(m.shape[0]):
        for j in range(m.shape[1]):
            if abs(m[i, j].real) < tol: s += "{0: {w}.0f}".format(0, w=real_width)
            else: s += "{0: {w}.{p}f}".format(m[i, j].real, w=real_width, p=prec)
            if abs(m[i, j].imag) < tol: s += "{0: >+{w}.0f}j".format(0, w=im_width)
            else: s += "{0: >+{w}.{p}f}j".format(m[i, j].imag, w=im_width, p=prec)
        s += "\n"
    return s