Python numpy.roots() Examples

The following are 30 code examples for showing how to use numpy.roots(). 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: pyberny   Author: jhrmnn   File: Math.py    License: Mozilla Public License 2.0 6 votes vote down vote up
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 2
Project: ocelot   Author: ocelot-collab   File: k_analysis.py    License: GNU General Public License v3.0 6 votes vote down vote up
def data_analysis(e_ph, flux, method="least"):

    if method == "least":
        coeffs = np.polyfit(x=e_ph, y=flux, deg=11)
        polynom = np.poly1d(coeffs)


        x = np.linspace(e_ph[0], e_ph[-1], num=100)
        pd = np.polyder(polynom, m=1)
        indx = np.argmax(np.abs(pd(x)))
        eph_c = x[indx]

        pd2 = np.polyder(polynom, m=2)
        p2_roots = np.roots(pd2)
        p2_roots = p2_roots[p2_roots[:].imag == 0]
        p2_roots = np.real(p2_roots)
        Eph_fin = find_nearest(p2_roots,eph_c)
        return Eph_fin, polynom

    elif method == "new method":
        pass

        #plt.plot(Etotal, total, "ro")
        #plt.plot(x, polynom(x)) 
Example 3
Project: parametric_modeling   Author: awesomebytes   File: impz.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def impz(b,a):
    """Pseudo implementation of the impz method of MATLAB"""
#% Compute time vector
# M = 0;  NN = [];
# if isempty(N)
#   % if not specified, determine the length
#   if isTF
#     N = impzlength(b,a,.00005);
#   else
#     N  = impzlength(b,.00005);
#   end
    p = np.roots(a)
    N = stableNmarginal_length(p, 0.00005, 0)
    N = len(b) * len(b) * len(b) # MATLAB AUTOFINDS THE SIZE HERE... 
    #TODO: Implement some way of finding the autosieze of this... I used a couple of examples... matlab gave 43 as length we give 64
    x = zeros(N)
    x[0] = 1
    h = lfilter(b,a, x)
    return h 
Example 4
Project: svgpathtools   Author: mathandy   File: polytools.py    License: MIT License 6 votes vote down vote up
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 5
Project: pyiron   Author: pyiron   File: thermo_bulk.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def get_minimum_energy_path(self, pressure=None):
        """

        Args:
            pressure:

        Returns:

        """
        if pressure is not None:
            raise NotImplemented()
        v_min_lst = []
        for c in self._coeff.T:
            v_min = np.roots(np.polyder(c, 1))
            p_der2 = np.polyder(c, 2)
            p_val2 = np.polyval(p_der2, v_min)
            v_m_lst = v_min[p_val2 > 0]
            if len(v_m_lst) > 0:
                v_min_lst.append(v_m_lst[0])
            else:
                v_min_lst.append(np.nan)
        return np.array(v_min_lst) 
Example 6
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 7
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_filter_design.py    License: MIT License 6 votes vote down vote up
def test_output_order(self):
        zc, zr = _cplxreal(np.roots(array([1, 0, 0, 1])))
        assert_allclose(np.append(zc, zr), [1/2 + 1j*sin(pi/3), -1])

        eps = spacing(1)

        a = [0+1j, 0-1j, eps + 1j, eps - 1j, -eps + 1j, -eps - 1j,
             1, 4, 2, 3, 0, 0,
             2+3j, 2-3j,
             1-eps + 1j, 1+2j, 1-2j, 1+eps - 1j,  # sorts out of order
             3+1j, 3+1j, 3+1j, 3-1j, 3-1j, 3-1j,
             2-3j, 2+3j]
        zc, zr = _cplxreal(a)
        assert_allclose(zc, [1j, 1j, 1j, 1+1j, 1+2j, 2+3j, 2+3j, 3+1j, 3+1j,
                             3+1j])
        assert_allclose(zr, [0, 0, 1, 2, 3, 4])

        z = array([1-eps + 1j, 1+2j, 1-2j, 1+eps - 1j, 1+eps+3j, 1-2*eps-3j,
                   0+1j, 0-1j, 2+4j, 2-4j, 2+3j, 2-3j, 3+7j, 3-7j, 4-eps+1j,
                   4+eps-2j, 4-1j, 4-eps+2j])

        zc, zr = _cplxreal(z)
        assert_allclose(zc, [1j, 1+1j, 1+2j, 1+3j, 2+3j, 2+4j, 3+7j, 4+1j,
                             4+2j])
        assert_equal(zr, []) 
Example 8
Project: pychemqt   Author: jjgomera   File: psycrometry.py    License: GNU General Public License v3.0 6 votes vote down vote up
def _v(self, P, T, phi_w):
        """Calculate the real volume of a humid air using the virial equation
        of state

        Parameters
        ----------
        T : float
            Temperature, [K]
        phi_w : float
            Molar fraction of water, [-]

        Returns
        -------
        """
        vir = self._virialMixture(T, phi_w)
        Bm = vir["Bm"]
        Cm = vir["Cm"]

        vm = roots([1, -R*T/P, -R*T*Bm/P, -R*T*Cm/P])
        if vm[0].imag == 0.0:
            v = vm[0].real
        else:
            v = vm[2].real

        return v 
Example 9
Project: pyvocoder   Author: shamidreza   File: lpc.py    License: GNU General Public License v2.0 6 votes vote down vote up
def poly2lsf(a):
    a = a / a[0]        
    A = np.r_[a, 0.0]
    B = A[::-1]
    P = A - B  
    Q = A + B  
    
    P = deconvolve(P, np.array([1.0, -1.0]))[0]
    Q = deconvolve(Q, np.array([1.0, 1.0]))[0]
    
    roots_P = np.roots(P)
    roots_Q = np.roots(Q)
    
    angles_P = np.angle(roots_P[::2])
    angles_Q = np.angle(roots_Q[::2])
    angles_P[angles_P < 0.0] += np.pi
    angles_Q[angles_Q < 0.0] += np.pi
    lsf = np.sort(np.r_[angles_P, angles_Q])
    return lsf 
Example 10
Project: MobileNet   Author: Zehaos   File: yellowfin.py    License: Apache License 2.0 6 votes vote down vote up
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 11
Project: antifier   Author: john-38787364   File: power_curve.py    License: MIT License 6 votes vote down vote up
def get_speed(power,Cx,f,W,slope,headwind,elevation):
    # slope in percent
    # headwind in m/s at 10 m high
    # elevation in meters
    air_pressure = 1 - 0.000104 * elevation
    Cx = Cx*air_pressure
    G = 9.81
    headwind = (0.1**0.143) * headwind
    roots = np.roots([Cx, 2*Cx*headwind, Cx*headwind**2 + W*G*(slope/100.0+f), -power])
    roots = np.real(roots[np.imag(roots) == 0])
    roots = roots[roots>0]
    speed = np.min(roots)
    if speed + headwind < 0:
        roots = np.roots([-Cx, -2*Cx*headwind, -Cx*headwind**2 + W*G*(slope/100.0+f), -power])
        roots = np.real(roots[np.imag(roots) == 0])
        roots = roots[roots>0]
        if len(roots) > 0:
            speed = np.min(roots)  
    return speed 
Example 12
Project: quadpy   Author: nschloe   File: _stroud_1966.py    License: GNU General Public License v3.0 6 votes vote down vote up
def stroud_1966_2(n):
    degree = 3
    # r is a smallest real-valued root of a polynomial of degree 3
    rts = numpy.roots(
        [2 * (n - 2) * (n + 1) * (n + 3), -(5 * n ** 2 + 5 * n - 18), 4 * n, -1]
    )
    r = numpy.min([r.real for r in rts if abs(r.imag) < 1.0e-15])

    s = 1 - n * r
    t = 0.5

    B = (n - 2) / (1 - 2 * n * r ** 2 - 2 * (1 - n * r) ** 2) / (n + 1) / (n + 2)
    C = 2 * (1 / (n + 1) - B) / n

    data = [(B, rd(n + 1, [(r, n), (s, 1)])), (C, rd(n + 1, [(t, 2)]))]

    points, weights = untangle(data)
    return TnScheme("Stroud 1966-II", n, weights, points, degree, source) 
Example 13
Project: pyAudioAnalysis   Author: tyiannak   File: ShortTermFeatures.py    License: Apache License 2.0 6 votes vote down vote up
def phormants(x, sampling_rate):
    N = len(x)
    w = np.hamming(N)

    # Apply window and high pass filter.
    x1 = x * w
    x1 = lfilter([1], [1., 0.63], x1)

    # Get LPC.
    ncoeff = 2 + sampling_rate / 1000
    A, e, k = lpc(x1, ncoeff)
    # A, e, k = lpc(x1, 8)

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

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

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

    return frqs 
Example 14
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 15
Project: pyberny   Author: jhrmnn   File: Math.py    License: Mozilla Public License 2.0 5 votes vote down vote up
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)

    # discriminant of d^2y/dx^2=0
    D = -((g0 + g1) ** 2) - 2 * g0 * g1 + 6 * (y1 - y0) * (g0 + g1) - 6 * (y1 - y0) ** 2
    if D < 1e-11:
        return None, None
    else:
        m = -5 * g0 - g1 - 6 * y0 + 6 * y1
        p1 = g(y0, y1, g0, g1, 0.5 * (m + np.sqrt(2 * D)))
        p2 = g(y0, y1, g0, g1, 0.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 16
Project: python-control   Author: python-control   File: margins.py    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def _polysqr(pol):
    """return a polynomial squared"""
    return np.polymul(pol, pol)

# Took the framework for the old function by
# Sawyer B. Fuller <minster@caltech.edu>, removed a lot of the innards
# and replaced with analytical polynomial functions for LTI systems.
#
# idea for the frequency data solution copied/adapted from
# https://github.com/alchemyst/Skogestad-Python/blob/master/BODE.py
# Rene van Paassen <rene.vanpaassen@gmail.com>
#
# RvP, July 8, 2014, corrected to exclude phase=0 crossing for the gain
#                    margin polynomial
# RvP, July 8, 2015, augmented to calculate all phase/gain crossings with
#                    frd data. Correct to return smallest phase
#                    margin, smallest gain margin and their frequencies
# RvP, Jun 10, 2017, modified the inclusion of roots found for phase
#                    crossing to include all >= 0, made subsequent calc
#                    insensitive to div by 0
#                    also changed the selection of which crossings to
#                    return on basis of "A note on the Gain and Phase
#                    Margin Concepts" Journal of Control and Systems
#                    Engineering, Yazdan Bavafi-Toosi, Dec 2015, vol 3
#                    issue 1, pp 51-59, closer to Matlab behavior, but
#                    not completely identical in edge cases, which don't
#                    cross but touch gain=1 
Example 17
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 18
Project: recruit   Author: Frank-qlu   File: test_polynomial.py    License: Apache License 2.0 5 votes vote down vote up
def test_roots(self):
        assert_array_equal(np.roots([1, 0, 0]), [0, 0]) 
Example 19
Project: ibllib   Author: int-brain-lab   File: metrics.py    License: MIT License 5 votes vote down vote up
def fp_est(ts, rp=0.002):
    '''
    An estimate of the fraction of false positive spikes in a unit
    (see Hill et al. (2011) J Neurosci 31: 8699-8705).

    Parameters
    ----------
    ts : ndarray_like
        The timestamps (in s) of the spikes.
    rp : float
        The refractory period (in s).

    Returns
    -------
    fp : float
        An estimate of the fraction of false positives.

    Examples
    --------
    1) Compute false positive estimate for unit 1.
        >>> ts = units_b['times']['1']
        >>> fp = bb.metrics.fp_est(ts)
    '''

    # Get number of spikes, number of isi violations, and time from first to final spike.
    n_spks = len(ts)
    n_isi_viol = len(np.where(np.diff(ts < rp)[0]))
    t = ts[-1] - ts[0]

    # `fp` is min of roots of solved quadratic equation.
    c = (t * n_isi_viol) / (2 * rp * n_spks**2)  # 3rd term in quadratic
    fp = np.min(np.abs(np.roots([-1, 1, c])))  # solve quadratic
    return fp 
Example 20
Project: YellowFin_MXNet   Author: rzhu3   File: yellowfin.py    License: Apache License 2.0 5 votes vote down vote up
def single_step_mu_lr(self, C, D, h_min, h_max):
    coef = np.array([-1.0, 3.0, 0.0, 1.0])
    coef[2] = -(3 + D ** 2 * h_min ** 2 / 2 / C)
    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 = h_max / h_min
    mu_t = max(np.real(root)[0] ** 2, ((np.sqrt(dr) - 1) / (np.sqrt(dr) + 1)) ** 2)
    lr_t = (1.0 - math.sqrt(mu_t)) ** 2 / h_min
    return mu_t, lr_t 
Example 21
Project: YellowFin_MXNet   Author: rzhu3   File: yellowfin_test.py    License: Apache License 2.0 5 votes vote down vote up
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 22
Project: lambda-packs   Author: ryfeus   File: common.py    License: MIT License 5 votes vote down vote up
def intersect_trust_region(x, s, Delta):
    """Find the intersection of a line with the boundary of a trust region.
    
    This function solves the quadratic equation with respect to t
    ||(x + s*t)||**2 = Delta**2.
    
    Returns
    -------
    t_neg, t_pos : tuple of float
        Negative and positive roots.
    
    Raises
    ------
    ValueError
        If `s` is zero or `x` is not within the trust region.
    """
    a = np.dot(s, s)
    if a == 0:
        raise ValueError("`s` is zero.")

    b = np.dot(x, s)

    c = np.dot(x, x) - Delta**2
    if c > 0:
        raise ValueError("`x` is not within the trust region.")

    d = np.sqrt(b*b - a*c)  # Root from one fourth of the discriminant.

    # Computations below avoid loss of significance, see "Numerical Recipes".
    q = -(b + copysign(d, b))
    t1 = q / a
    t2 = c / q

    if t1 < t2:
        return t1, t2
    else:
        return t2, t1 
Example 23
Project: lambda-packs   Author: ryfeus   File: test_polynomial.py    License: MIT License 5 votes vote down vote up
def test_roots(self):
        assert_array_equal(np.roots([1, 0, 0]), [0, 0]) 
Example 24
Project: auto-alt-text-lambda-api   Author: abhisuri97   File: test_polynomial.py    License: MIT License 5 votes vote down vote up
def test_roots(self):
        assert_array_equal(np.roots([1, 0, 0]), [0, 0]) 
Example 25
Project: tensor2tensor   Author: tensorflow   File: yellowfin_test.py    License: Apache License 2.0 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 26
Project: CEASIOMpy   Author: cfsengineering   File: func_dynamic.py    License: Apache License 2.0 5 votes vote down vote up
def direc_mode_characteristic(roll,spiral,dr1,dr2):
    """ TODO:

    roll,spiral,dr1,dr2 roots of characteristic equation

    Args:

    Returns:
        roll_timecst,
        spiral_timecst,
        spiral_t2: time to double amplitude
        dr_freq  [rad/s]
        dr_damp
        dr_damp_freq (): product of freq and damping ratio[rad/s]

    """

    # Roll
    roll_timecst = - 1/ roll

    # Spiral
    spiral_timecst = -1/spiral
    spiral_t2 = ln(2) / spiral

    # Dutch Roll
    dr_freq = (np.sqrt(dr1*dr2)).real # [rad/s] Undamped natural frequency omega, of concise equations
    dr_damp = - dr1.real /dr_freq # [-] Damping of concise equations
    dr_damp_freq = dr_damp*dr_freq # [rad/s] time constante  T = 1/(XiOmega)

    return(roll_timecst, spiral_timecst, spiral_t2, dr_freq, dr_damp, dr_damp_freq) 
Example 27
Project: vnpy_crypto   Author: birforce   File: test_polynomial.py    License: MIT License 5 votes vote down vote up
def test_roots(self):
        assert_array_equal(np.roots([1, 0, 0]), [0, 0]) 
Example 28
Project: vnpy_crypto   Author: birforce   File: sarimax.py    License: MIT License 5 votes vote down vote up
def arroots(self):
        """
        (array) Roots of the reduced form autoregressive lag polynomial
        """
        return np.roots(self.polynomial_reduced_ar)**-1 
Example 29
Project: vnpy_crypto   Author: birforce   File: sarimax.py    License: MIT License 5 votes vote down vote up
def maroots(self):
        """
        (array) Roots of the reduced form moving average lag polynomial
        """
        return np.roots(self.polynomial_reduced_ma)**-1 
Example 30
Project: vnpy_crypto   Author: birforce   File: sarimax.py    License: MIT License 5 votes vote down vote up
def arfreq(self):
        """
        (array) Frequency of the roots of the reduced form autoregressive
        lag polynomial
        """
        z = self.arroots
        if not z.size:
            return
        return np.arctan2(z.imag, z.real) / (2 * np.pi)