Python numpy.real() Examples

The following are 30 code examples of numpy.real(). 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 also want to check out all available functions/classes of the module numpy , or try the search function .
Example #1
Source Project: pyscf   Author: pyscf   File: linalg_helper.py    License: Apache License 2.0 7 votes vote down vote up
def diagonalize_asymm(H):
    """
    Diagonalize a real, *asymmetric* matrix and return sorted results.

    Return the eigenvalues and eigenvectors (column matrix)
    sorted from lowest to highest eigenvalue.
    """
    E,C = np.linalg.eig(H)
    #if np.allclose(E.imag, 0*E.imag):
    #    E = np.real(E)
    #else:
    #    print "WARNING: Eigenvalues are complex, will be returned as such."

    idx = E.real.argsort()
    E = E[idx]
    C = C[:,idx]

    return E,C 
Example #2
Source Project: simnibs   Author: simnibs   File: grid.py    License: GNU General Public License v3.0 7 votes vote down vote up
def quadrature_cc_1D(N):
    """ Computes the Clenshaw Curtis nodes and weights """
    N = np.int(N)        
    if N == 1:
        knots = 0
        weights = 2
    else:
        n = N - 1
        C = np.zeros((N,2))
        k = 2*(1+np.arange(np.floor(n/2)))
        C[::2,0] = 2/np.hstack((1, 1-k*k))
        C[1,1] = -n
        V = np.vstack((C,np.flipud(C[1:n,:])))
        F = np.real(ifft(V, n=None, axis=0))
        knots = F[0:N,1]
        weights = np.hstack((F[0,0],2*F[1:n,0],F[n,0]))
            
    return knots, weights 
Example #3
Source Project: xrft   Author: xgcm   File: test_xrft.py    License: MIT License 6 votes vote down vote up
def test_dft_real_2d(self):
        """
        Test the real discrete Fourier transform function on one-dimensional
        data. Non-trivial because we need to keep only some of the negative
        frequencies.
        """
        Nx, Ny = 16, 32
        da = xr.DataArray(np.random.rand(Nx, Ny), dims=['x', 'y'],
                          coords={'x': range(Nx), 'y': range(Ny)})
        dx = float(da.x[1] - da.x[0])
        dy = float(da.y[1] - da.y[0])

        daft = xrft.dft(da, real='x')
        npt.assert_almost_equal(daft.values,
                               np.fft.rfftn(da.transpose('y','x')).transpose())
        npt.assert_almost_equal(daft.values,
                               xrft.dft(da, dim=['y'], real='x'))

        actual_freq_x = daft.coords['freq_x'].values
        expected_freq_x = np.fft.rfftfreq(Nx, dx)
        npt.assert_almost_equal(actual_freq_x, expected_freq_x)

        actual_freq_y = daft.coords['freq_y'].values
        expected_freq_y = np.fft.fftfreq(Ny, dy)
        npt.assert_almost_equal(actual_freq_y, expected_freq_y) 
Example #4
Source Project: FRIDA   Author: LCAV   File: point_cloud.py    License: MIT License 6 votes vote down vote up
def classical_mds(self, D):
        ''' 
        Classical multidimensional scaling

        Parameters
        ----------
        D : square 2D ndarray
            Euclidean Distance Matrix (matrix containing squared distances between points
        '''

        # Apply MDS algorithm for denoising
        n = D.shape[0]
        J = np.eye(n) - np.ones((n,n))/float(n)
        G = -0.5*np.dot(J, np.dot(D, J))

        s, U = np.linalg.eig(G)

        # we need to sort the eigenvalues in decreasing order
        s = np.real(s)
        o = np.argsort(s)
        s = s[o[::-1]]
        U = U[:,o[::-1]]

        S = np.diag(s)[0:self.dim,:]
        self.X = np.dot(np.sqrt(S),U.T) 
Example #5
Source Project: transferlearning   Author: jindongwang   File: EasyTL.py    License: MIT License 6 votes vote down vote up
def get_cosine_dist(A, B):
    B = np.reshape(B, (1, -1))
    
    if A.shape[1] == 1:
        A = np.hstack((A, np.zeros((A.shape[0], 1))))
        B = np.hstack((B, np.zeros((B.shape[0], 1))))
    
    aa = np.sum(np.multiply(A, A), axis=1).reshape(-1, 1)
    bb = np.sum(np.multiply(B, B), axis=1).reshape(-1, 1)
    ab = A @ B.T
    
    # to avoid NaN for zero norm
    aa[aa==0] = 1
    bb[bb==0] = 1
    
    D = np.real(np.ones((A.shape[0], B.shape[0])) - np.multiply((1/np.sqrt(np.kron(aa, bb.T))), ab))
    
    return D 
Example #6
Source Project: scarlet   Author: pmelchior   File: interpolation.py    License: MIT License 6 votes vote down vote up
def fft_convolve(*images):
    """Use FFT's to convove an image with a kernel

    Parameters
    ----------
    images: list of array-like
        A list of images to convolve.

    Returns
    -------
    result: array
        The convolution in pixel space of `img` with `kernel`.
    """
    from autograd.numpy.numpy_boxes import ArrayBox

    Images = [np.fft.fft2(np.fft.ifftshift(img)) for img in images]
    if np.any([isinstance(img, ArrayBox) for img in images]):
        Convolved = Images[0]
        for img in Images[1:]:
            Convolved = Convolved * img
    else:
        Convolved = np.prod(Images, 0)
    convolved = np.fft.ifft2(Convolved)
    return np.fft.fftshift(np.real(convolved)) 
Example #7
Source Project: pyscf   Author: pyscf   File: linalg_helper.py    License: Apache License 2.0 6 votes vote down vote up
def __init__(self,matr_multiply,xStart,inPreCon,nroots=1,tol=1e-10):
        self.matrMultiply = matr_multiply
        self.size = xStart.shape[0]
        self.nEigen = min(nroots, self.size)
        self.maxM = min(30, self.size)
        self.maxOuterLoop = 10
        self.tol = tol

        #
        #  Creating initial guess and preconditioner
        #
        self.x0 = xStart.real.copy()

        self.iteration = 0
        self.totalIter = 0
        self.converged = False
        self.preCon = inPreCon.copy()
        #
        #  Allocating other vectors
        #
        self.allocateVecs() 
Example #8
Source Project: pyscf   Author: pyscf   File: linalg_helper.py    License: Apache License 2.0 6 votes vote down vote up
def solveSubspace(self):
        w, v = scipy.linalg.eig(self.subH[:self.currentSize,:self.currentSize])
        idx = w.real.argsort()
        v = v[:,idx]
        w = w[idx].real
#
        imag_norm = np.linalg.norm(w.imag)
        if imag_norm > 1e-12:
            print(" *************************************************** ")
            print(" WARNING  IMAGINARY EIGENVALUE OF NORM %.15g " % (imag_norm))
            print(" *************************************************** ")
        #print "Imaginary norm eigenvectors = ", np.linalg.norm(v.imag)
        #print "Imaginary norm eigenvalue   = ", np.linalg.norm(w.imag)
        #print "eigenvalues = ", w[:min(self.currentSize,7)]
#
        self.sol[:self.currentSize] = v[:,self.ciEig]
        self.evecs[:self.currentSize,:self.currentSize] = v
        self.eigs[:self.currentSize] = w[:self.currentSize]
        self.outeigs[:self.nEigen] = w[:self.nEigen]
        self.cvEig = self.eigs[self.ciEig] 
Example #9
Source Project: pyqmc   Author: WagnerGroup   File: ewald.py    License: MIT License 6 votes vote down vote up
def energy(self, configs):
        r"""
        Compute Coulomb energy for a set of configs.  

        .. math:: E_{\rm Coulomb} &= E_{\rm real+reciprocal}^{ee} 
                + E_{\rm self+charged}^{ee} 
                \\&+ E_{\rm real+reciprocal}^{e\text{-ion}} 
                + E_{\rm self+charged}^{e\text{-ion}} 
                \\&+ E_{\rm real+reciprocal}^{\text{ion-ion}} 
                + E_{\rm self+charged}^{\text{ion-ion}}
        
        Inputs:
            configs: pyqmc PeriodicConfigs object of shape (nconf, nelec, ndim)
        Returns: 
            ee: electron-electron part
            ei: electron-ion part
            ii: ion-ion part
        """
        nelec = configs.configs.shape[1]
        ee, ei = self.ewald_electron(configs)
        ee += self.ee_const(nelec)
        ei += self.ei_const(nelec)
        ii = self.ion_ion + self.ii_const
        return ee, ei, ii 
Example #10
Source Project: python-control   Author: python-control   File: lti.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def damp(self):
        '''Natural frequency, damping ratio of system poles

        Returns
        -------
        wn : array
            Natural frequencies for each system pole
        zeta : array
            Damping ratio for each system pole
        poles : array
            Array of system poles
        '''
        poles = self.pole()

        if isdtime(self, strict=True):
            splane_poles = np.log(poles)/self.dt
        else:
            splane_poles = poles
        wn = absolute(splane_poles)
        Z = -real(splane_poles)/wn
        return wn, Z, poles 
Example #11
Source Project: python-control   Author: python-control   File: rlocus.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def _break_points(num, den):
    """Extract break points over real axis and gains given these locations"""
    # type: (np.poly1d, np.poly1d) -> (np.array, np.array)
    dnum = num.deriv(m=1)
    dden = den.deriv(m=1)
    polynom = den * dnum - num * dden
    real_break_pts = polynom.r
    # don't care about infinite break points
    real_break_pts = real_break_pts[num(real_break_pts) != 0]
    k_break = -den(real_break_pts) / num(real_break_pts)
    idx = k_break >= 0   # only positives gains
    k_break = k_break[idx]
    real_break_pts = real_break_pts[idx]
    if len(k_break) == 0:
        k_break = [0]
        real_break_pts = den.roots
    return k_break, real_break_pts 
Example #12
Source Project: python-control   Author: python-control   File: steering-gainsched.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def control_output(t, x, u, params):
    # Get the controller parameters
    longpole = params.get('longpole', -2.)
    latpole1 = params.get('latpole1', -1/2 + sqrt(-7)/2)
    latpole2 = params.get('latpole2', -1/2 - sqrt(-7)/2)
    l = params.get('wheelbase', 3)
    
    # Extract the system inputs
    ex, ey, etheta, vd, phid = u

    # Determine the controller gains
    alpha1 = -np.real(latpole1 + latpole2)
    alpha2 = np.real(latpole1 * latpole2)

    # Compute and return the control law
    v = -longpole * ex          # Note: no feedfwd (to make plot interesting)
    if vd != 0:
        phi = phid + (alpha1 * l) / vd * ey + (alpha2 * l) / vd * etheta
    else:
        # We aren't moving, so don't turn the steering wheel
        phi = phid
    
    return  np.array([v, phi])

# Define the controller as an input/output system 
Example #13
Source Project: python-control   Author: python-control   File: steering-gainsched.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def control_output(t, x, u, params):
    # Get the controller parameters
    longpole = params.get('longpole', -2.)
    latpole1 = params.get('latpole1', -1/2 + sqrt(-7)/2)
    latpole2 = params.get('latpole2', -1/2 - sqrt(-7)/2)
    l = params.get('wheelbase', 3)
    
    # Extract the system inputs
    ex, ey, etheta, vd, phid = u

    # Determine the controller gains
    alpha1 = -np.real(latpole1 + latpole2)
    alpha2 = np.real(latpole1 * latpole2)

    # Compute and return the control law
    v = -longpole * ex          # Note: no feedfwd (to make plot interesting)
    if vd != 0:
        phi = phid + (alpha1 * l) / vd * ey + (alpha2 * l) / vd * etheta
    else:
        # We aren't moving, so don't turn the steering wheel
        phi = phid
    
    return  np.array([v, phi])

# Define the controller as an input/output system 
Example #14
Source Project: pyGSTi   Author: pyGSTio   File: polynomial.py    License: Apache License 2.0 6 votes vote down vote up
def compact(self, complex_coeff_tape=True):
        """
        Generate a compact form of this polynomial designed for fast evaluation.

        The resulting "tapes" can be evaluated using
        :function:`opcalc.bulk_eval_compact_polys`.

        Parameters
        ----------
        complex_coeff_tape : bool, optional
            Whether the `ctape` returned array is forced to be of complex type.
            If False, the real part of all coefficients is taken (even if they're
            complex).

        Returns
        -------
        vtape, ctape : numpy.ndarray
            These two 1D arrays specify an efficient means for evaluating this
            polynomial.
        """
        if complex_coeff_tape:
            return self._rep.compact_complex()
        else:
            return self._rep.compact_real() 
Example #15
Source 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 #16
Source Project: pyGSTi   Author: pyGSTio   File: reportableqty.py    License: Apache License 2.0 6 votes vote down vote up
def absdiff(self, constant_value, separate_re_im=False):
        """
        Returns a ReportableQty that is the (element-wise in the vector case)
        difference between `constant_value` and this one given by:

        `abs(self - constant_value)`.
        """
        if separate_re_im:
            re_v = _np.fabs(_np.real(self.value) - _np.real(constant_value))
            im_v = _np.fabs(_np.imag(self.value) - _np.imag(constant_value))
            if self.has_eb():
                return (ReportableQty(re_v, _np.fabs(_np.real(self.errbar)), self.nonMarkovianEBs),
                        ReportableQty(im_v, _np.fabs(_np.imag(self.errbar)), self.nonMarkovianEBs))
            else:
                return ReportableQty(re_v), ReportableQty(im_v)

        else:
            v = _np.absolute(self.value - constant_value)
            if self.has_eb():
                return ReportableQty(v, _np.absolute(self.errbar), self.nonMarkovianEBs)
            else:
                return ReportableQty(v) 
Example #17
Source 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 #18
Source 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 #19
Source Project: xrft   Author: xgcm   File: test_xrft.py    License: MIT License 5 votes vote down vote up
def test_dft_real_1d(self, test_data_1d):
        """
        Test the discrete Fourier transform function on one-dimensional data.
        """
        da = test_data_1d
        Nx = len(da)
        dx = float(da.x[1] - da.x[0]) if 'x' in da.dims else 1

        # defaults with no keyword args
        ft = xrft.dft(da, real='x', detrend='constant')
        # check that the frequency dimension was created properly
        assert ft.dims == ('freq_x',)
        # check that the coords are correct
        freq_x_expected = np.fft.rfftfreq(Nx, dx)
        npt.assert_allclose(ft['freq_x'], freq_x_expected)
        # check that a spacing variable was created
        assert ft['freq_x_spacing'] == freq_x_expected[1] - freq_x_expected[0]
        # make sure the function is lazy
        assert isinstance(ft.data, type(da.data))
        # check that the Fourier transform itself is correct
        data = (da - da.mean()).values
        ft_data_expected = np.fft.rfft(data)
        # because the zero frequency component is zero, there is a numerical
        # precision issue. Fixed by setting atol
        npt.assert_allclose(ft_data_expected, ft.values, atol=1e-14)

        with pytest.raises(ValueError):
            xrft.dft(da, real='y', detrend='constant') 
Example #20
Source Project: xrft   Author: xgcm   File: test_xrft.py    License: MIT License 5 votes vote down vote up
def test_power_spectrum(self, dask):
        """Test the power spectrum function"""
        N = 16
        da = xr.DataArray(np.random.rand(2,N,N), dims=['time','y','x'],
                         coords={'time':np.array(['2019-04-18', '2019-04-19'],
                                                dtype='datetime64'),
                                'y':range(N),'x':range(N)}
                         )
        if dask:
            da = da.chunk({'time': 1})
        ps = xrft.power_spectrum(da, dim=['y','x'], window=True, density=False,
                                detrend='constant')
        daft = xrft.dft(da, dim=['y','x'], detrend='constant', window=True)
        npt.assert_almost_equal(ps.values, np.real(daft*np.conj(daft)))
        npt.assert_almost_equal(np.ma.masked_invalid(ps).mask.sum(), 0.)

        ps = xrft.power_spectrum(da, dim=['y'], real='x', window=True,
                                density=False, detrend='constant')
        daft = xrft.dft(da, dim=['y'], real='x', detrend='constant', window=True)
        npt.assert_almost_equal(ps.values, np.real(daft*np.conj(daft)))

        ### Normalized
        ps = xrft.power_spectrum(da, dim=['y','x'], window=True, detrend='constant')
        daft = xrft.dft(da, dim=['y','x'], window=True, detrend='constant')
        test = np.real(daft*np.conj(daft))/N**4
        dk = np.diff(np.fft.fftfreq(N, 1.))[0]
        test /= dk**2
        npt.assert_almost_equal(ps.values, test)
        npt.assert_almost_equal(np.ma.masked_invalid(ps).mask.sum(), 0.)

        ### Remove least-square fit
        ps = xrft.power_spectrum(da, dim=['y','x'],
                                window=True, density=False, detrend='linear'
                                )
        daft = xrft.dft(da, dim=['y','x'], window=True, detrend='linear')
        npt.assert_almost_equal(ps.values, np.real(daft*np.conj(daft)))
        npt.assert_almost_equal(np.ma.masked_invalid(ps).mask.sum(), 0.) 
Example #21
Source Project: xrft   Author: xgcm   File: test_xrft.py    License: MIT License 5 votes vote down vote up
def test_cross_spectrum(self, dask):
        """Test the cross spectrum function"""
        N = 16
        dim = ['x','y']
        da = xr.DataArray(np.random.rand(2,N,N), dims=['time','x','y'],
                          coords={'time':np.array(['2019-04-18', '2019-04-19'],
                                                 dtype='datetime64'),
                                 'x':range(N), 'y':range(N)})
        da2 = xr.DataArray(np.random.rand(2,N,N), dims=['time','x','y'],
                          coords={'time':np.array(['2019-04-18', '2019-04-19'],
                                                 dtype='datetime64'),
                                  'x':range(N), 'y':range(N)})
        if dask:
            da = da.chunk({'time': 1})
            da2 = da2.chunk({'time': 1})

        daft = xrft.dft(da, dim=dim, shift=True, detrend='constant',
                        window=True)
        daft2 = xrft.dft(da2, dim=dim, shift=True, detrend='constant',
                        window=True)
        cs = xrft.cross_spectrum(da, da2, dim=dim, window=True, density=False,
                                detrend='constant')
        npt.assert_almost_equal(cs.values, np.real(daft*np.conj(daft2)))
        npt.assert_almost_equal(np.ma.masked_invalid(cs).mask.sum(), 0.)

        cs = xrft.cross_spectrum(da, da2, dim=dim, shift=True, window=True,
                                detrend='constant')
        test = (daft * np.conj(daft2)).real.values/N**4

        dk = np.diff(np.fft.fftfreq(N, 1.))[0]
        test /= dk**2
        npt.assert_almost_equal(cs.values, test)
        npt.assert_almost_equal(np.ma.masked_invalid(cs).mask.sum(), 0.) 
Example #22
Source 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 #23
Source Project: dc_tts   Author: Kyubyong   File: utils.py    License: Apache License 2.0 5 votes vote down vote up
def griffin_lim(spectrogram):
    '''Applies Griffin-Lim's raw.'''
    X_best = copy.deepcopy(spectrogram)
    for i in range(hp.n_iter):
        X_t = invert_spectrogram(X_best)
        est = librosa.stft(X_t, hp.n_fft, hp.hop_length, win_length=hp.win_length)
        phase = est / np.maximum(1e-8, np.abs(est))
        X_best = spectrogram * phase
    X_t = invert_spectrogram(X_best)
    y = np.real(X_t)

    return y 
Example #24
Source Project: GST-Tacotron   Author: KinglittleQ   File: utils.py    License: MIT License 5 votes vote down vote up
def griffin_lim(spectrogram):
    '''Applies Griffin-Lim's raw.
    '''
    X_best = copy.deepcopy(spectrogram)
    for i in range(hp.n_iter):
        X_t = invert_spectrogram(X_best)
        est = librosa.stft(X_t, hp.n_fft, hp.hop_length, win_length=hp.win_length)
        phase = est / np.maximum(1e-8, np.abs(est))
        X_best = spectrogram * phase
    X_t = invert_spectrogram(X_best)
    y = np.real(X_t)

    return y 
Example #25
Source Project: OpenFermion-Cirq   Author: quantumlib   File: fermionic_simulation.py    License: Apache License 2.0 5 votes vote down vote up
def _canonicalize_weight(w):
    if w == 0:
        return (0, 0)
    if cirq.is_parameterized(w):
        return (cirq.PeriodicValue(abs(w), 2 * sympy.pi), sympy.arg(w))
    period = 2 * np.pi
    return (np.round((w.real % period) if (w == np.real(w)) else
                     (abs(w) % period) * w / abs(w), 8), 0) 
Example #26
Source Project: OpenFermion-Cirq   Author: quantumlib   File: low_rank.py    License: Apache License 2.0 5 votes vote down vote up
def __init__(self,
                 hamiltonian: openfermion.InteractionOperator,
                 truncation_threshold: Optional[float]=1e-8,
                 final_rank: Optional[int]=None,
                 spin_basis=True) -> None:

        self.truncation_threshold = truncation_threshold
        self.final_rank = final_rank

        # Perform the low rank decomposition of two-body operator.
        self.eigenvalues, self.one_body_squares, one_body_correction, _ = (
            openfermion.low_rank_two_body_decomposition(
                hamiltonian.two_body_tensor,
                truncation_threshold=self.truncation_threshold,
                final_rank=self.final_rank,
                spin_basis=spin_basis))

        # Get scaled density-density terms and basis transformation matrices.
        self.scaled_density_density_matrices = []  # type: List[numpy.ndarray]
        self.basis_change_matrices = []            # type: List[numpy.ndarray]
        for j in range(len(self.eigenvalues)):
            density_density_matrix, basis_change_matrix = (
                openfermion.prepare_one_body_squared_evolution(
                    self.one_body_squares[j]))
            self.scaled_density_density_matrices.append(
                    numpy.real(self.eigenvalues[j] * density_density_matrix))
            self.basis_change_matrices.append(basis_change_matrix)

        # Get transformation matrix and orbital energies for one-body terms
        one_body_coefficients = (
                hamiltonian.one_body_tensor + one_body_correction)
        quad_ham = openfermion.QuadraticHamiltonian(one_body_coefficients)
        self.one_body_energies, self.one_body_basis_change_matrix, _ = (
                quad_ham.diagonalizing_bogoliubov_transform()
        )

        super().__init__(hamiltonian) 
Example #27
Source Project: fine-lm   Author: akzaidi   File: yellowfin_test.py    License: MIT License 5 votes vote down vote up
def tune_everything(self, x0squared, c, t, gmin, gmax):
    del t
    # First tune based on dynamic range
    if c == 0:
      dr = gmax / gmin
      mustar = ((np.sqrt(dr) - 1) / (np.sqrt(dr) + 1))**2
      alpha_star = (1 + np.sqrt(mustar))**2/gmax

      return alpha_star, mustar

    dist_to_opt = x0squared
    grad_var = c
    max_curv = gmax
    min_curv = gmin
    const_fact = dist_to_opt * min_curv**2 / 2 / grad_var
    coef = [-1, 3, -(3 + const_fact), 1]
    roots = np.roots(coef)
    roots = roots[np.real(roots) > 0]
    roots = roots[np.real(roots) < 1]
    root = roots[np.argmin(np.imag(roots))]

    assert root > 0 and root < 1 and np.absolute(root.imag) < 1e-6

    dr = max_curv / min_curv
    assert max_curv >= min_curv
    mu = max(((np.sqrt(dr) - 1) / (np.sqrt(dr) + 1))**2, root**2)

    lr_min = (1 - np.sqrt(mu))**2 / min_curv

    alpha_star = lr_min
    mustar = mu

    return alpha_star, mustar 
Example #28
Source Project: transferlearning   Author: jindongwang   File: intra_alignment.py    License: MIT License 5 votes vote down vote up
def getAngle(Ps, Pt, DD):
    
    Q = np.hstack((Ps, scipy.linalg.null_space(Ps.T)))
    dim = Pt.shape[1]
    QPt = Q.T @ Pt
    A, B = QPt[:dim, :], QPt[dim:, :]
    U,V,X,C,S = gsvd(A, B)
    alpha = np.zeros([1, DD])
    for i in range(DD):
        alpha[0][i] = math.sin(np.real(math.acos(C[i][i]*math.pi/180)))
    
    return alpha 
Example #29
Source Project: kss   Author: Kyubyong   File: utils.py    License: Apache License 2.0 5 votes vote down vote up
def griffin_lim(spectrogram):
    '''Applies Griffin-Lim's raw.'''
    X_best = copy.deepcopy(spectrogram)
    for i in range(hp.n_iter):
        X_t = invert_spectrogram(X_best)
        est = librosa.stft(X_t, hp.n_fft, hp.hop_length, win_length=hp.win_length)
        phase = est / np.maximum(1e-8, np.abs(est))
        X_best = spectrogram * phase
    X_t = invert_spectrogram(X_best)
    y = np.real(X_t)

    return y 
Example #30
Source Project: controlpy   Author: markwmuller   File: analysis.py    License: GNU General Public License v3.0 5 votes vote down vote up
def is_hurwitz(A, tolerance = 1e-9):
    '''Test whether the matrix A is Hurwitz (i.e. asymptotically stable).
    
    tolerance defines the minimum distance we should be from the imaginary axis 
     to be considered stable.
    
    '''
    return max(np.real(np.linalg.eig(A)[0])) < -np.abs(tolerance)