Python numpy.roots() 使用实例

The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.

Example 1

def phormants(x, Fs):
    N = len(x)
    w = numpy.hamming(N)

    # Apply window and high pass filter.
    x1 = x * w   
    x1 = lfilter([1], [1., 0.63], x1)
    
    # Get LPC.    
    ncoeff = 2 + Fs / 1000
    A, e, k = lpc(x1, ncoeff)    
    #A, e, k = lpc(x1, 8)

    # Get roots.
    rts = numpy.roots(A)
    rts = [r for r in rts if numpy.imag(r) >= 0]

    # Get angles.
    angz = numpy.arctan2(numpy.imag(rts), numpy.real(rts))

    # Get frequencies.    
    frqs = sorted(angz * (Fs / (2 * math.pi)))

    return frqs 

Example 2

def lpf(x, cutoff, fs, order=5):
    """
    low pass filters signal with Butterworth digital
    filter according to cutoff frequency

    filter uses Gustafsson’s method to make sure
    forward-backward filt == backward-forward filt

    Note that edge effects are expected

    Args:
        x      (array): signal data (numpy array)
        cutoff (float): cutoff frequency (Hz)
        fs       (int): sample rate (Hz)
        order    (int): order of filter (default 5)

    Returns:
        filtered (array): low pass filtered data
    """
    nyquist = fs / 2
    b, a = butter(order, cutoff / nyquist)
    if not np.all(np.abs(np.roots(a)) < 1):
        raise PsolaError('Filter with cutoff at {} Hz is unstable given '
                         'sample frequency {} Hz'.format(cutoff, fs))
    filtered = filtfilt(b, a, x, method='gust')
    return filtered 

Example 3

def estimate_time_constant(y, p=2, sn=None, lags=5, fudge_factor=1.):
    """
    Estimate AR model parameters through the autocovariance function

    Parameters
    ----------
    y : array, shape (T,)
        One dimensional array containing the fluorescence intensities with
        one entry per time-bin.
    p : positive integer
        order of AR system
    sn : float
        sn standard deviation, estimated if not provided.
    lags : positive integer
        number of additional lags where he autocovariance is computed
    fudge_factor : float (0< fudge_factor <= 1)
        shrinkage factor to reduce bias

    Returns
    -------
    g : estimated coefficients of the AR process
    """

    if sn is None:
        sn = GetSn(y)

    lags += p
    xc = axcov(y, lags)
    xc = xc[:, np.newaxis]

    A = scipy.linalg.toeplitz(xc[lags + np.arange(lags)],
                              xc[lags + np.arange(p)]) - sn**2 * np.eye(lags, p)
    g = np.linalg.lstsq(A, xc[lags + 1:])[0]
    gr = np.roots(np.concatenate([np.array([1]), -g.flatten()]))
    gr = (gr + gr.conjugate()) / 2.
    gr[gr > 1] = 0.95 + np.random.normal(0, 0.01, np.sum(gr > 1))
    gr[gr < 0] = 0.15 + np.random.normal(0, 0.01, np.sum(gr < 0))
    g = np.poly(fudge_factor * gr)
    g = -g[1:]

    return g.flatten() 

Example 4

def find_closest(t, v, t0, v0):
    """ Find the closest point on the curve f = a + b/x
    to the given point (t,v)
    """
    a = v0
    b = v0*t0
    # Solve for intersection points
    eqn_coefs = [1/b, -t/b, 0, v-a, -b]
    tis = np.roots(eqn_coefs)
    tis = tis[abs(tis.imag/tis.real)<0.01].real # We care only real solutions
    tis = tis[tis>0] # and positive ones
    # Choose the shortest among solutions
    ds = abs(tis-t)*np.sqrt(1 + np.power(tis,4)/(b*b)) # Distance from solutions to given point (t,v)
    idx = np.argmin(ds)
    ti = tis[idx]
    vi = a + b/ti
    return ti, vi 

Example 5

def get_Z(self, T, p):

        kappa = 0.37464 + 1.54226*self.acentric - 0.26992*self.acentric**2

        Tr = T/self.T_crit
        alpha = (1 + kappa*(1 - Tr**0.5))**2

        A = alpha * self.a * p / self.R**2 / T**2
        B = self.b * p / self.R / T

        # Solve the cubic equation for compressibility factor z
        coeffs = [1, -(1-B), A-2*B-3*B**2, -(A*B-B**2-B**3)]
        #print(coeffs)
        roots = np.roots(coeffs)

        real_roots = roots[np.isreal(roots)].real
        valid_roots = real_roots[real_roots > p*self.b/self.R/T]

        return valid_roots 

Example 6

def fit_cubic(y0, y1, g0, g1):
    """Fit cubic polynomial to function values and derivatives at x = 0, 1.

    Returns position and function value of minimum if fit succeeds. Fit does
    not succeeds if

    1. polynomial doesn't have extrema or
    2. maximum is from (0,1) or
    3. maximum is closer to 0.5 than minimum
    """
    a = 2*(y0-y1)+g0+g1
    b = -3*(y0-y1)-2*g0-g1
    p = np.array([a, b, g0, y0])
    r = np.roots(np.polyder(p))
    if not np.isreal(r).all():
        return None, None
    r = sorted(x.real for x in r)
    if p[0] > 0:
        maxim, minim = r
    else:
        minim, maxim = r
    if 0 < maxim < 1 and abs(minim-0.5) > abs(maxim-0.5):
        return None, None
    return minim, np.polyval(p, minim) 

Example 7

def get_mu_tensor(self):
    const_fact = self._dist_to_opt_avg**2 * self._h_min**2 / 2 / self._grad_var
    coef = tf.Variable([-1.0, 3.0, 0.0, 1.0], dtype=tf.float32, name="cubic_solver_coef")
    coef = tf.scatter_update(coef, tf.constant(2), -(3 + const_fact) )        
    roots = tf.py_func(np.roots, [coef], Tout=tf.complex64, stateful=False)
    
    # filter out the correct root
    root_idx = tf.logical_and(tf.logical_and(tf.greater(tf.real(roots), tf.constant(0.0) ),
      tf.less(tf.real(roots), tf.constant(1.0) ) ), tf.less(tf.abs(tf.imag(roots) ), 1e-5) )
    # in case there are two duplicated roots satisfying the above condition
    root = tf.reshape(tf.gather(tf.gather(roots, tf.where(root_idx) ), tf.constant(0) ), shape=[] )
    tf.assert_equal(tf.size(root), tf.constant(1) )

    dr = self._h_max / self._h_min
    mu = tf.maximum(tf.real(root)**2, ( (tf.sqrt(dr) - 1)/(tf.sqrt(dr) + 1) )**2)    
    return mu 

Example 8

def polyroots(p, realroots=False, condition=lambda r: True):
    """
    Returns the roots of a polynomial with coefficients given in p.
      p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
    INPUT:
    p - Rank-1 array-like object of polynomial coefficients.
    realroots - a boolean.  If true, only real roots will be returned  and the
        condition function can be written assuming all roots are real.
    condition - a boolean-valued function.  Only roots satisfying this will be
        returned.  If realroots==True, these conditions should assume the roots
        are real.
    OUTPUT:
    A list containing the roots of the polynomial.
    NOTE:  This uses np.isclose and np.roots"""
    roots = np.roots(p)
    if realroots:
        roots = [r.real for r in roots if isclose(r.imag, 0)]
    roots = [r for r in roots if condition(r)]

    duplicates = []
    for idx, (r1, r2) in enumerate(combinations(roots, 2)):
        if isclose(r1, r2):
            duplicates.append(idx)
    return [r for idx, r in enumerate(roots) if idx not in duplicates] 

Example 9

def polyroots(p, realroots=False, condition=lambda r: True):
    """
    Returns the roots of a polynomial with coefficients given in p.
      p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
    INPUT:
    p - Rank-1 array-like object of polynomial coefficients.
    realroots - a boolean.  If true, only real roots will be returned  and the
        condition function can be written assuming all roots are real.
    condition - a boolean-valued function.  Only roots satisfying this will be
        returned.  If realroots==True, these conditions should assume the roots
        are real.
    OUTPUT:
    A list containing the roots of the polynomial.
    NOTE:  This uses np.isclose and np.roots"""
    roots = np.roots(p)
    if realroots:
        roots = [r.real for r in roots if isclose(r.imag, 0)]
    roots = [r for r in roots if condition(r)]

    duplicates = []
    for idx, (r1, r2) in enumerate(combinations(roots, 2)):
        if isclose(r1, r2):
            duplicates.append(idx)
    return [r for idx, r in enumerate(roots) if idx not in duplicates] 

Example 10

def calc_cp(alpha,beta,k):
    gamma = []
    for i in range(k+1):
        num = factorial(alpha+k)*factorial(alpha+beta+k+i)
        denom = factorial(alpha+i)*factorial(k-i)*factorial(i)
        gamma.insert(i,num/denom)

    poly = []
    for i in range(k+1):
        if i == 0:
            poly.insert(i,gamma[i])
        else:
            prod = [1]
            j=1
            while j<=i:
                prod=conv(prod,[1,-1])
                j=j+1
            while len(poly)<len(prod):
                poly.insert(0,0)
            prod = [gamma[i]*t for t in prod]
            poly = [sum(pair) for pair in zip(poly,prod)]

    cp = numpy.roots(poly)
    return cp 

Example 11

def pitch_from_relden(relden, cf, sw):
    """
    This function calculates the pitch of cuboct of a given relative density, chamfer factor, and strut width.
    :param relden: float. Desired relative density
    :param cf: float. Chamfer factor of voxel
    :param sw: float. strut width of voxel
    :return: lattice pitch
    """
    chamheight = float(sw) / cf
    l_2 = sw / 2.0 + chamheight
    l_3 = l_2 + sw * np.cos(math.radians(45))  # horizontal position of points
    l_4 = np.sqrt(2) * (l_3 - sw / 2.0)
    tan_theta = ((l_3 - l_2) / ((l_4 / 2.0) - (np.sqrt(2) * chamheight / 2.0)))
    v1 = l_2 * (sw * sw + 4 * sw * (l_3 - sw / 2.0) + 2 * (l_3 - sw / 2.0) * (l_3 - sw / 2.0))
    h = (l_4 / 2.0) * tan_theta
    hs = chamheight * tan_theta * np.sqrt(2) / 2.0
    v2 = ((l_4 * l_4 * h) - (2 * (chamheight * chamheight * hs))) / 3.0
    v3 = 4 * sw * (0.5 * (l_3 - l_2) * (l_3 - l_2) + (l_3 - l_2) * chamheight)
    v4 = sw * sw * (l_3 - l_2)
    node_volume = v1 + v2 + v3 + v4

    c1 = relden
    c2 = (-6) * np.sqrt(2)*sw *sw
    c3 = -6*node_volume + 12*sw*sw*np.sqrt(2)*(l_2 + l_3)
    return max(np.roots([c1, 0, c2, c3])) 

Example 12

def tune_everything(x0squared, C, T, gmin, gmax):
  # 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
  lr_max = (1 + np.sqrt(mu) )**2 / max_curv

  alpha_star = lr_min
  mustar = mu

  return alpha_star, mustar 

Example 13

def test_roots(self):
        assert_array_equal(np.roots([1, 0, 0]), [0, 0]) 

Example 14

def GetRootsPolynomial(a, b, xE, yE, xC, yC, r):
    '''
    
    '''
    
    # Define some stuff
    r2 = r * r;
    a2 = a * a;
    b2 = b * b;
    a2b2 = a2 / b2;
    x0 = xE - xC;
    y0 = yE - yC;
    y2 = y0 * y0;
    x2 = x0 * x0;
    
    # Get the coefficients
    A = a2b2 - 1.;
    B = -2. * x0 * a2b2;
    C = r2 - y2 - a2 + a2b2 * x2;
    D = 4. * y2 * a2b2;
    c4 = A * A;
    c3 = 2. * A * B;
    c2 = 2. * A * C + B * B + D;
    c1 = 2. * B * C - 2. * D * x0;
    c0 = C * C - (b2 - x2) * D;
        
    # Get the real roots
    roots = [r.real + xC for r in np.roots([c4, c3, c2, c1, c0]) 
             if np.abs(r.imag) < tol]
    return roots 

Example 15

def filter_is_stable(a):
    """
    Check if filter coefficients of IIR filter are stable.
    
    Parameters
    ----------
    a: list or 1darray of number
        Denominator filter coefficients a.

    Returns
    -------
    is_stable: bool
        Filter is stable or not.  
    Notes
    ----
    Filter is stable if absolute value of all  roots is smaller than 1,
    see [1]_.
    
    References
    ----------
    .. [1] HYRY, "SciPy 'lfilter' returns only NaNs" StackOverflow,
       http://stackoverflow.com/a/8812737/1469195
    """
    assert a[0] == 1.0, (
        "a[0] should normally be zero, did you accidentally supply b?\n"
        "a: {:s}".format(str(a)))
    # from http://stackoverflow.com/a/8812737/1469195
    return np.all(np.abs(np.roots(a))<1) 

Example 16

def tune_everything(x0squared, C, T, gmin, gmax):
  # 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
  lr_max = (1 + np.sqrt(mu) )**2 / max_curv

  alpha_star = lr_min
  mustar = mu

  return alpha_star, mustar 

Example 17

def get_Z(self, T, p):
        """
        Returns the compressibility factor of a real gas.

        :param T: temperature (K)
        :type T: float
        :param p: pressure (Pa)
        :type p: float

        :return: compressibility factor
        :rtype: float
        """
        # a = self.get_a(T)
        # coeffs = [
        #         1,
        #         (-self.R*T - self.b*p + self.c*p + self.d*p)/(self.R*T),
        #         (-self.R*T*self.d*p + a*p - self.b*self.d*p**2 + self.c*self.d*p**2 + self.e*p**2)/(self.R**2*T**2),
        #         (-self.R*T*self.e*p**2 - a*self.b*p**2 + a*self.c*p**2 - self.b*self.e*p**3 + self.c*self.e*p**3)/(self.R**3*T**3)
        # ]
        coeffs = [1, self.get_A(T, p), self.get_B(T, p), self.get_C(T, p)]
        #print(coeffs)

        roots = np.roots(coeffs)
        real_roots = roots[np.isreal(roots)].real

        if len(real_roots) == 1:
            real_roots = real_roots[0]

        return real_roots

    # Partial derivative of Z wrt T at constant p 

Example 18

def get_Z(self, T, p):

        A = self.get_A(T, p)
        B = self.get_B(T, p)

        # Solve the cubic equation for compressibility factor z
        coeffs = [1, -1, A-B-B**2, -A*B]
        roots = np.roots(coeffs)

        real_roots = roots[np.isreal(roots)].real
        valid_roots = real_roots[real_roots > p*self.b/self.R/T]

        return valid_roots 

Example 19

def get_Z(self, T, p):
        #A = self.a * p / self.R**2 / T**2.5
        #B = self.b * p / self.R / T
        A = self.get_A(T, p)
        B = self.get_B(T, p)

        # Solve the cubic equation for compressibility factor z
        # Z^3 - Z^2 + (A-B-B**2)*Z - A*B = 0
        coeffs = [1, -1, A-B-B**2, -A*B]
        roots = np.roots(coeffs)

        real_roots = roots[np.isreal(roots)].real
        valid_roots = real_roots[real_roots > p*self.b/self.R/T]

        return valid_roots

    # dZ/dT at const. p 

Example 20

def lpc_to_lsf(all_lpc):
    if len(all_lpc.shape) < 2:
        all_lpc = all_lpc[None]
    order = all_lpc.shape[1] - 1
    all_lsf = np.zeros((len(all_lpc), order))
    for i in range(len(all_lpc)):
        lpc = all_lpc[i]
        lpc1 = np.append(lpc, 0)
        lpc2 = lpc1[::-1]
        sum_filt = lpc1 + lpc2
        diff_filt = lpc1 - lpc2

        if order % 2 != 0:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, 0, -1])
            deconv_sum = sum_filt
        else:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, -1])
            deconv_sum, _ = sg.deconvolve(sum_filt, [1, 1])

        roots_diff = np.roots(deconv_diff)
        roots_sum = np.roots(deconv_sum)
        angle_diff = np.angle(roots_diff[::2])
        angle_sum = np.angle(roots_sum[::2])
        lsf = np.sort(np.hstack((angle_diff, angle_sum)))
        if len(lsf) != 0:
            all_lsf[i] = lsf
    return np.squeeze(all_lsf) 

Example 21

def lpc_to_lsf(all_lpc):
    if len(all_lpc.shape) < 2:
        all_lpc = all_lpc[None]
    order = all_lpc.shape[1] - 1
    all_lsf = np.zeros((len(all_lpc), order))
    for i in range(len(all_lpc)):
        lpc = all_lpc[i]
        lpc1 = np.append(lpc, 0)
        lpc2 = lpc1[::-1]
        sum_filt = lpc1 + lpc2
        diff_filt = lpc1 - lpc2

        if order % 2 != 0:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, 0, -1])
            deconv_sum = sum_filt
        else:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, -1])
            deconv_sum, _ = sg.deconvolve(sum_filt, [1, 1])

        roots_diff = np.roots(deconv_diff)
        roots_sum = np.roots(deconv_sum)
        angle_diff = np.angle(roots_diff[::2])
        angle_sum = np.angle(roots_sum[::2])
        lsf = np.sort(np.hstack((angle_diff, angle_sum)))
        if len(lsf) != 0:
            all_lsf[i] = lsf
    return np.squeeze(all_lsf) 

Example 22

def test_roots(self):
        assert_array_equal(np.roots([1, 0, 0]), [0, 0]) 

Example 23

def fit_quartic(y0, y1, g0, g1):
    """Fit constrained quartic polynomial to function values and erivatives at x = 0,1.

    Returns position and function value of minimum or None if fit fails or has
    a maximum. Quartic polynomial is constrained such that it's 2nd derivative
    is zero at just one point. This ensures that it has just one local
    extremum.  No such or two such quartic polynomials always exist. From the
    two, the one with lower minimum is chosen.
    """
    def g(y0, y1, g0, g1, c):
        a = c+3*(y0-y1)+2*g0+g1
        b = -2*c-4*(y0-y1)-3*g0-g1
        return np.array([a, b, c, g0, y0])

    def quart_min(p):
        r = np.roots(np.polyder(p))
        is_real = np.isreal(r)
        if is_real.sum() == 1:
            minim = r[is_real][0].real
        else:
            minim = r[(r == max(-abs(r))) | r == -max(-abs(r))][0].real
        return minim, np.polyval(p, minim)

    D = -(g0+g1)**2-2*g0*g1+6*(y1-y0)*(g0+g1)-6*(y1-y0)**2  # discriminant of d^2y/dx^2=0
    if D < 1e-11:
        return None, None
    else:
        m = -5*g0-g1-6*y0+6*y1
        p1 = g(y0, y1, g0, g1, .5*(m+np.sqrt(2*D)))
        p2 = g(y0, y1, g0, g1, .5*(m-np.sqrt(2*D)))
        if p1[0] < 0 and p2[0] < 0:
            return None, None
        [minim1, minval1] = quart_min(p1)
        [minim2, minval2] = quart_min(p2)
        if minval1 < minval2:
            return minim1, minval1
        else:
            return minim2, minval2 

Example 24

def signalMinimum(img, fitParams=None, n_std=3):
    '''
    intersection between signal and background peak
    '''
    if fitParams is None:
        fitParams = FitHistogramPeaks(img).fitParams
    assert len(fitParams) > 1, 'need 2 peaks so get minimum signal'

    i = signalPeakIndex(fitParams)
    signal = fitParams[i]
    bg = getBackgroundPeak(fitParams)
    smn = signal[1] - n_std * signal[2]
    bmx = bg[1] + n_std * bg[2]
    if smn > bmx:
        return smn
    # peaks are overlapping
    # define signal min. as intersection between both Gaussians

    def solve(p1, p2):
        s1, m1, std1 = p1
        s2, m2, std2 = p2
        a = (1 / (2 * std1**2)) - (1 / (2 * std2**2))
        b = (m2 / (std2**2)) - (m1 / (std1**2))
        c = (m1**2 / (2 * std1**2)) - (m2**2 / (2 * std2**2)) - \
            np.log(((std2 * s1) / (std1 * s2)))
        return np.roots([a, b, c])
    i = solve(bg, signal)
    try:
        return i[np.logical_and(i > bg[1], i < signal[1])][0]
    except IndexError:
        # this error shouldn't occur... well
        return max(smn, bmx) 

Example 25

def getSignalMinimum(fitParams, n_std=3):
    assert len(fitParams) > 0, 'need min. 1 peak so get minimum signal'
    if len(fitParams) == 1:
        signal = fitParams[0]
        return signal[1] - n_std * signal[2]

    i = signalPeakIndex(fitParams)
    signal = fitParams[i]
    bg = fitParams[i - 1]
    #bg = getBackgroundPeak(fitParams)
    smn = signal[1] - n_std * signal[2]
    bmx = bg[1] + n_std * bg[2]
    if smn > bmx:
        return smn
    # peaks are overlapping
    # define signal min. as intersection between both Gaussians

    def solve(p1, p2):
        s1, m1, std1 = p1
        s2, m2, std2 = p2
        a = (1 / (2 * std1**2)) - (1 / (2 * std2**2))
        b = (m2 / (std2**2)) - (m1 / (std1**2))
        c = (m1**2 / (2 * std1**2)) - (m2**2 / (2 * std2**2)) - \
            np.log(((std2 * s1) / (std1 * s2)))
        return np.roots([a, b, c])

    i = solve(bg, signal)
    try:
        return i[np.logical_and(i > bg[1], i < signal[1])][0]
    except IndexError:
        # something didnt work out - fallback
        return smn 

Example 26

def tuneEverything(self, x0squared, c, t, gmin, gmax):
        # 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 27

def update_X(self, X, mu, k=20):
        U, S, VT = svdp(X, k=k)
        P = np.c_[np.ones((k, 1)), 1-S, 1./2./mu-S]
        sigma_star = np.zeros(k)
        for t in range(k):
            p = P[t, :]
            delta = p[1]**2 - 4 * p[0] * p[2]
            if delta <= 0:
                sigma_star[t] = 0.
            else:
                solution = np.roots(p)
                solution = sorted(solution, key=abs)
                solution = np.array(solution)
                if solution[0] * solution[1] <= 0:
                    sigma_star[t] = solution[1]
                elif solution[1] < 0:
                    sigma_star[t] = 0.
                else:
                    f = np.log(1 + solution[1]) + mu * (solution[1] - s[t])**2
                    if f > mu * s[t]**2:
                        sigma_star[t] = 0.
                    else:
                        sigma_star[t] = solution[1]

        sigma_star = np.diag(sigma_star)
        sigma_star = np.dot(np.dot(U, sigma_star), VT)
        return sigma_star 

Example 28

def update_X(X, mu, k=6):
    #U, S, VT = svdp(X, k=k)
    U, S, VT = svds(X, k=k, which='LM')
    P = np.c_[np.ones((k, 1)), 1-S, 1./2./mu-S]
    sigma_star = np.zeros(k)
    for t in range(k):
        p = P[t, :]
        delta = p[1]**2 - 4 * p[0] * p[2]
        if delta <= 0:
            sigma_star[t] = 0.
        else:
            solution = np.roots(p)
            solution = solution.tolist()
            solution.sort(key=abs)
            solution = np.array(solution)
            if solution[0] * solution[1] <= 0:
                sigma_star[t] = solution[1]
            elif solution[1] < 0:
                sigma_star[t] = 0.
            else:
                f = np.log(1 + solution[1]) + mu * (solution[1] - s[t])**2
                if f > mu *s[t]**2:
                    sigma_star[t] = 0.
                else:
                    sigma_star[t] = solution[1]

    sigma_star = sp.csr_matrix(np.diag(sigma_star))
    sigma_star = safe_sparse_dot(safe_sparse_dot(U, sigma_star), VT)
    sigma_star[abs(sigma_star)<1e-10] = 0
    return sp.lil_matrix(sigma_star) 

Example 29

def real_roots(coeffs):
    """Get real roots of a polynomial.

    Args:
        coeffs (List[Float]): List of polynomial coefficients.

    Returns:
        numpy.ndarray: The (sorted) real roots of the polynomial.
    """
    all_roots = np.roots(coeffs)
    filtered = all_roots[all_roots.imag == 0.0].real
    return np.sort(filtered) 

Example 30

def roots(a):

    return matlabarray(np.roots(np.asarray(a).ravel())) 

Example 31

def estimate_time_constant(fluor, p = 2, sn = None, lags = 5, fudge_factor = 1.):
    """    
    Estimate AR model parameters through the autocovariance function    
    Inputs
    ----------
    fluor        : nparray
        One dimensional array containing the fluorescence intensities with
        one entry per time-bin.
    p            : positive integer
        order of AR system  
    sn           : float
        noise standard deviation, estimated if not provided.
    lags         : positive integer
        number of additional lags where he autocovariance is computed
    fudge_factor : float (0< fudge_factor <= 1)
        shrinkage factor to reduce bias
        
    Return
    -----------
    g       : estimated coefficients of the AR process
    """    
    

    if sn is None:
        sn = GetSn(fluor)
        
    lags += p
    xc = axcov(fluor,lags)        
    xc = xc[:,np.newaxis]
    
    A = scipy.linalg.toeplitz(xc[lags+np.arange(lags)],xc[lags+np.arange(p)]) - sn**2*np.eye(lags,p)
    g = np.linalg.lstsq(A,xc[lags+1:])[0]
    gr = np.roots(np.concatenate([np.array([1]),-g.flatten()]))
    gr = (gr+gr.conjugate())/2.
    gr[gr>1] = 0.95 + np.random.normal(0,0.01,np.sum(gr>1))
    gr[gr<0] = 0.15 + np.random.normal(0,0.01,np.sum(gr<0))
    g = np.poly(fudge_factor*gr)
    g = -g[1:]    
        
    return g.flatten() 

Example 32

def _compute_minimum(self, a1, p1, dp1, a2, p2, dp2):
        cubic_mtx = np.zeros((4, 4))
        cubic_mtx[0, :] = [1., a1, a1 ** 2, a1 ** 3]
        cubic_mtx[1, :] = [1., a2, a2 ** 2, a2 ** 3]
        cubic_mtx[2, :] = [0., 1., 2 * a1, 3 * a1 ** 2]
        cubic_mtx[3, :] = [0., 1., 2 * a2, 3 * a2 ** 2]
        c0, c1, c2, c3 = np.linalg.solve(cubic_mtx, [p1, p2, dp1, dp2])

        d0 = c1
        d1 = 2 * c2
        d2 = 3 * c3
        r1, r2 = np.roots([d2, d1, d0])

        a = None
        p = max(p1, p2)
        if (a1 <= r1 <= a2 or a2 <= r1 <= a1) and np.isreal(r1):
            px = self._phi(r1)
            if px < p:
                a = r1
                p = px
                dp = self._dphi(r1)
        if (a1 <= r2 <= a2 or a2 <= r2 <= a1) and np.isreal(r2):
            px = self._phi(r2)
            if px < p:
                a = r2
                p = px
                dp = self._dphi(r2)

        return a, p, dp 

Example 33

def roots_of_characteristic(self):
        """
        This function calculates z_0 and the 2m roots of the characteristic equation 
        associated with the Euler equation (1.7)
        
        Note:
        ------
        numpy.poly1d(roots, True) defines a polynomial using its roots that can be
        evaluated at any point. If x_1, x_2, ... , x_m are the roots then 
            p(x) = (x - x_1)(x - x_2)...(x - x_m)
        
        """
        m = self.m
        ? = self.?
        
        # Calculate the roots of the 2m-polynomial
        roots = np.roots(?)
        # sort the roots according to their length (in descending order) 
        roots_sorted = roots[np.argsort(abs(roots))[::-1]]    

        z_0 = ?.sum() / np.poly1d(roots, True)(1)
        z_1_to_m = roots_sorted[:m]     # we need only those outside the unit circle
        
        ? = 1 / z_1_to_m

        return z_1_to_m, z_0, ? 

Example 34

def test_roots(self):
        assert_array_equal(np.roots([1, 0, 0]), [0, 0]) 

Example 35

def tune_everything(x0squared, C, T, gmin, gmax):
  # 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
  lr_max = (1 + np.sqrt(mu) )**2 / max_curv

  alpha_star = lr_min
  mustar = mu

  return alpha_star, mustar 

Example 36

def tfps(beta=0.6, ca=1., de2=0.000545, theta=0., kk=1.):
   """
      beta: Total plasma beta,
      ca: Alfven speed based on mean field, 
      de2: me/mi, 
      theta: Angle of propagation in degrees
      kk: Wavenumber of interest in units of kdi
      
      Output is frequencies of the roots and the phase speeds w/k
      The roots are w[0]:Fast/Whistler, w[1]:Alfven/KAW, w[2]: Slow/Cyclotron
   """
   tht = theta*pi/180.
   ct = np.cos(tht)
   st = np.sin(tht)
   tt = st/ct
   cs=np.sqrt(beta/2.)*ca
   di =1.
   caksq=np.cos(tht)**2*ca**2
   cmsq=ca**2 + cs**2
   
   pcs=np.zeros(4)
   D = 1 + kk**2*de2
   # Find out the root of the quadratic dispersion relation
   pcs[0] = 1.
   pcs[1] = -(ca**2/D + cs**2 + caksq*(1+kk**2*di**2/D)/D)*kk**2
   pcs[2] = caksq*kk**4*(ca**2/D + cs**2 + cs**2*(1+kk**2*di**2/D))/D
   pcs[3] = -(cs**2*kk**6*caksq**2)/D**2
   w2 = np.roots(pcs); w = np.sqrt(w2)
   speeds= w/kk
   return w,speeds 

Example 37

def tfps(beta=0.6, ca=1., de2=0.000545, theta=0., kk=1.):
   """
      beta: Total plasma beta,
      ca: Alfven speed based on mean field, 
      de2: me/mi, 
      theta: Angle of propagation as fraction of pi/2), 
      kk: Wavenumber of interest in units of kdi
      
      Output is frequencies of the roots and the phase speeds w/k
      The roots are w[0]:Fast/Whistler, w[1]:Alfven/KAW, w[2]: Slow/Cyclotron
   """
#  theta=float(raw_input("Angle of propagation as fraction of pi/2: ")) 
#  kk=float(raw_input("What k value? "))
#  ca = float(raw_input("Alfven Speed: "))
#  beta=float(raw_input("Total plasma beta?: "))
#  de2= float(raw_input("me/mi : "))

   ct = cos(theta*pi/2)
   st = sin(theta*pi/2)
   tt = st/ct
   cs=sqrt(beta/2.)*ca
   di =1.
   caksq=cos(theta*pi/2.)**2*ca**2
   cmsq=ca**2 + cs**2
   
   pcs=np.zeros(4)
   D = 1 + kk**2*de2
   # Find out the root of the quadratic dispersion relation
   pcs[0] = 1.
   pcs[1] = -(ca**2/D + cs**2 + caksq*(1+kk**2*di**2/D)/D)*kk**2
   pcs[2] = caksq*kk**4*(ca**2/D + cs**2 + cs**2*(1+kk**2*di**2/D))/D
   pcs[3] = -(cs**2*kk**6*caksq**2)/D**2
   w2 = np.roots(pcs); w = sqrt(w2)
   speeds= w/kk
   return w,speeds 

Example 38

def get_mu(self):
    coef = [-1.0, 3.0, 0.0, 1.0]
    coef[2] = -(3 + self._dist_to_opt**2 * self._h_min**2 / 2 / self._grad_var)
    roots = np.roots(coef)
    root = roots[np.logical_and(np.logical_and(np.real(roots) > 0.0, 
      np.real(roots) < 1.0), np.imag(roots) < 1e-5) ]
    assert root.size == 1
    dr = self._h_max / self._h_min
    self._mu_t = max(np.real(root)[0]**2, ( (np.sqrt(dr) - 1) / (np.sqrt(dr) + 1) )**2 )
    return 

Example 39

def test_roots(self):
        assert_array_equal(np.roots([1, 0, 0]), [0, 0]) 

Example 40

def get_H(self, mass_upper):
        """Solve for H+, the concentration of hydrogen ions
        (the (pH) of seawater).

        :param mass_upper: Carbon mass in ocenas in GtC
        :type mass_upper: float
        :return: H
        :rtype: float
        """
        p0 = 1
        p1 = (self.k_1 - mass_upper * self.k_1 / self.Alk)
        p2 = (1 - 2 * mass_upper / self.Alk) * self.k_1 * self.k_2
        return max(np.roots([p0, p1, p2])) 

Example 41

def get_H(self, mass_upper):
        """Solve for H+, the concentration of hydrogen ions
        (the (pH) of seawater).

        :param mass_upper: Carbon mass in ocenas in GtC
        :type mass_upper: float
        :return: H
        :rtype: float
        """
        p0 = 1
        p1 = (self.k_1 - mass_upper * self.k_1 / self.Alk)
        p2 = (1 - 2 * mass_upper / self.Alk) * self.k_1 * self.k_2
        return max(np.roots([p0, p1, p2])) 

Example 42

def test_roots(self):
        assert_array_equal(np.roots([1, 0, 0]), [0, 0]) 

Example 43

def polyroots01(p):
    """Returns the real roots between 0 and 1 of the polynomial with
    coefficients given in p,
      p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
    p can also be a np.poly1d object.  See polyroots for more information."""
    return polyroots(p, realroots=True, condition=lambda tval: 0 <= tval <= 1) 

Example 44

def polyroots01(p):
    """Returns the real roots between 0 and 1 of the polynomial with
    coefficients given in p,
      p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
    p can also be a np.poly1d object.  See polyroots for more information."""
    return polyroots(p, realroots=True, condition=lambda tval: 0 <= tval <= 1) 

Example 45

def circleby2ptsradius(pt1Vec, pt2Vec, radius):
    
    if radius < la.norm(pt1Vec-pt2Vec)/2:
        centerVec1 = np.array(['NaN', 'NaN'])
        centerVec2 = np.array(['NaN', 'NaN'])
    else:
        a = pt1Vec[0]**2-pt2Vec[0]**2
        b = pt1Vec[1]**2-pt2Vec[1]**2
        c = -2*(pt1Vec[0]-pt2Vec[0])
        d = -2*(pt1Vec[1]-pt2Vec[1])
        e = a+b
        Coeff1 = 1+(d/c)**2
        Coeff2 = ((2*d*e/c**2) +(pt2Vec[0]*2*d/c) -2*pt2Vec[1])
        Coeff3 = ((e/c)**2+(2*pt2Vec[0]*e/c)+pt2Vec[0]**2+pt2Vec[1]**2-\
            radius**2)
        
        All_coeff = np.array([Coeff1, Coeff2, Coeff3])
        Eq_root = np.roots(All_coeff)
        
        # centersArray--center of the circle. It is a 2x2 matrix. First row
        # represents first possible center (x1, y1) and second row is the 
        # second possible center.
        centersArray = np.zeros((len(Eq_root),2))
        
        for i in list(range(len(Eq_root))):
            x = -(e+d*Eq_root[i])/c
            centersArray[i,0] = x
            centersArray[i,1] = Eq_root[i]
        
        # Check if values in centerArray are real or unreal numbers
        if la.norm(centersArray.imag) != 0:
            print('Circle with the specified radius does not pass through',
                'specified points', sep=" ")
            centerVec1 = np.array(['NaN', 'NaN'])
            centerVec2 = np.array(['NaN', 'NaN'])
        else:
            centerVec1 = centersArray[0, :]
            centerVec2 = centersArray[1, :]
        
    return centerVec1, centerVec2 

Example 46

def tuneEverything(self, x0squared, c, t, gmin, gmax):
    # 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 47

def check_mcs_step_septum(septum, proton, step_in):
    """
    From K2 scattering routine. It is basically treatment of edge effect inside the septum jaws
    @param septum: septum object
    @param proton: Proton object
    @param step_in: Random step chosen for this iteration
    """

    theta_0 = 13.6e-3 / (proton.mom_p) / np.sqrt(septum.rad_length) * np.sqrt(step_in) * (1. + 0.038 * np.log(step_in / septum.rad_length))



    x_temp = proton.x - septum.aper_u if proton.x <= (septum.aper_u + septum.thickness / 2.) else proton.x - septum.aper_u - septum.thickness

    a = 9 * theta_0 ** 2 / 3.
    b = -1 * proton.x_p ** 2
    c = -2 * proton.x_p * x_temp
    d = -x_temp ** 2
    # print 'a', a, 'b', b, 'c', c, 'd', d
    # test = np.linspace(-10, 10, 100)
    # test_y = a * test ** 3 + b * test ** 2 + test * c + d
    # plt.plot(test, test_y, 'lime')
    sol = np.roots([a, b, c, d])
    # print sol
    sol_re = [100000.]
    for i in sol:
        if np.imag(i) == 0:
            sol_re.append(float(np.real(i)))
    min_sol = min(sol_re)
    if min_sol < step_in:
        if min_sol < septum.rad_length * 0.1e-2:
            return septum.rad_length * 0.1e-2
        else:
            return min(sol_re)
    else:
        return step_in 

Example 48

def lpc_to_lsf(all_lpc):
    if len(all_lpc.shape) < 2:
        all_lpc = all_lpc[None]
    order = all_lpc.shape[1] - 1
    all_lsf = np.zeros((len(all_lpc), order))
    for i in range(len(all_lpc)):
        lpc = all_lpc[i]
        lpc1 = np.append(lpc, 0)
        lpc2 = lpc1[::-1]
        sum_filt = lpc1 + lpc2
        diff_filt = lpc1 - lpc2

        if order % 2 != 0:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, 0, -1])
            deconv_sum = sum_filt
        else:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, -1])
            deconv_sum, _ = sg.deconvolve(sum_filt, [1, 1])

        roots_diff = np.roots(deconv_diff)
        roots_sum = np.roots(deconv_sum)
        angle_diff = np.angle(roots_diff[::2])
        angle_sum = np.angle(roots_sum[::2])
        lsf = np.sort(np.hstack((angle_diff, angle_sum)))
        if len(lsf) != 0:
            all_lsf[i] = lsf
    return np.squeeze(all_lsf) 

Example 49

def lpc_to_lsf(all_lpc):
    if len(all_lpc.shape) < 2:
        all_lpc = all_lpc[None]
    order = all_lpc.shape[1] - 1
    all_lsf = np.zeros((len(all_lpc), order))
    for i in range(len(all_lpc)):
        lpc = all_lpc[i]
        lpc1 = np.append(lpc, 0)
        lpc2 = lpc1[::-1]
        sum_filt = lpc1 + lpc2
        diff_filt = lpc1 - lpc2

        if order % 2 != 0:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, 0, -1])
            deconv_sum = sum_filt
        else:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, -1])
            deconv_sum, _ = sg.deconvolve(sum_filt, [1, 1])

        roots_diff = np.roots(deconv_diff)
        roots_sum = np.roots(deconv_sum)
        angle_diff = np.angle(roots_diff[::2])
        angle_sum = np.angle(roots_sum[::2])
        lsf = np.sort(np.hstack((angle_diff, angle_sum)))
        if len(lsf) != 0:
            all_lsf[i] = lsf
    return np.squeeze(all_lsf) 

Example 50

def lpc_to_lsf(all_lpc):
    if len(all_lpc.shape) < 2:
        all_lpc = all_lpc[None]
    order = all_lpc.shape[1] - 1
    all_lsf = np.zeros((len(all_lpc), order))
    for i in range(len(all_lpc)):
        lpc = all_lpc[i]
        lpc1 = np.append(lpc, 0)
        lpc2 = lpc1[::-1]
        sum_filt = lpc1 + lpc2
        diff_filt = lpc1 - lpc2

        if order % 2 != 0:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, 0, -1])
            deconv_sum = sum_filt
        else:
            deconv_diff, _ = sg.deconvolve(diff_filt, [1, -1])
            deconv_sum, _ = sg.deconvolve(sum_filt, [1, 1])

        roots_diff = np.roots(deconv_diff)
        roots_sum = np.roots(deconv_sum)
        angle_diff = np.angle(roots_diff[::2])
        angle_sum = np.angle(roots_sum[::2])
        lsf = np.sort(np.hstack((angle_diff, angle_sum)))
        if len(lsf) != 0:
            all_lsf[i] = lsf
    return np.squeeze(all_lsf) 
点赞

发表评论

电子邮件地址不会被公开。 必填项已用*标注