Python numpy.roots() Examples

The following are 30 code examples of numpy.roots(). 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 File: power_curve.py    From antifier with 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 #2
Source File: psycrometry.py    From pychemqt with 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 #3
Source File: impz.py    From parametric_modeling with 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
Source File: k_analysis.py    From ocelot with 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 #5
Source File: lpc.py    From pyvocoder with 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 #6
Source File: yellowfin.py    From MobileNet with 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 #7
Source File: _stroud_1966.py    From quadpy with 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 #8
Source File: polytools.py    From svgpathtools with 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 #9
Source File: test_filter_design.py    From GraphicDesignPatternByPython with 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 #10
Source File: thermo_bulk.py    From pyiron with 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 #11
Source File: ShortTermFeatures.py    From pyAudioAnalysis with 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 #12
Source File: control_and_filter.py    From QuantEcon.lectures.code with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
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 #13
Source File: Math.py    From pyberny with 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 #14
Source File: _stroud_1964.py    From quadpy with GNU General Public License v3.0 5 votes vote down vote up
def _stroud_1964(variant_a, n):
    degree = 3

    if n == 2:
        # The roots sum up to 1
        r, s, t = mp.polyroots([1, -1, 0.25, -1.0 / 60.0])
        data = [(1.0 / (n * (n + 1)), rd(n + 1, [(r, 1), (s, 1), (t, 1)]))]
    else:
        assert n > 2

        # Stroud's book only gives numerical values for certain n; the article explains
        # it in more detail, namely: r is a root of a polynomial of degree 3.
        rts = numpy.sort(
            numpy.roots([n + 1, -3, 3 / (n + 2), -1 / ((n + 2) * (n + 3))])
        )

        # all roots are real-valued
        if n > 8:
            assert not variant_a, "Choose variant b for n >= 9."

        r = rts[0] if variant_a else rts[1]

        # s and t are zeros of a polynomial of degree 2
        s, t = numpy.sort(
            numpy.roots(
                [
                    1,
                    -(1 - (n - 1) * r),
                    n / (2 * (n + 2)) - (n - 1) * r + n * (n - 1) / 2 * r ** 2,
                ]
            )
        )

        data = [(1 / (n * (n + 1)), rd(n + 1, [(r, n - 1), (s, 1), (t, 1)]))]

    points, weights = untangle(data)

    name = "Stroud 1964{}".format("a" if variant_a else "b")
    return TnScheme(name, n, weights, points, degree, source) 
Example #15
Source File: bandpass_filter.py    From spiketoolkit with MIT License 5 votes vote down vote up
def __init__(self, recording, freq_min=300, freq_max=6000, freq_wid=1000, filter_type='fft', order=3,
                 chunk_size=30000, cache_chunks=False):
        assert HAVE_BFR, "To use the BandpassFilterRecording, install scipy: \n\n pip install scipy\n\n"
        self._freq_min = freq_min
        self._freq_max = freq_max
        self._freq_wid = freq_wid
        self._type = filter_type
        self._order = order
        self._chunk_size = chunk_size

        if self._type == 'butter':
            fn = recording.get_sampling_frequency() / 2.
            band = np.array([self._freq_min, self._freq_max]) / fn

            self._b, self._a = ss.butter(self._order, band, btype='bandpass')

            if not np.all(np.abs(np.roots(self._a)) < 1):
                raise ValueError('Filter is not stable')
        FilterRecording.__init__(self, recording=recording, chunk_size=chunk_size, cache_chunks=cache_chunks)
        self.copy_channel_properties(recording)

        self.is_filtered = True
        self._kwargs = {'recording': recording.make_serialized_dict(), 'freq_min': freq_min, 'freq_max': freq_max,
                        'freq_wid': freq_wid, 'filter_type': filter_type, 'order': order,
                        'chunk_size': chunk_size, 'cache_chunks': cache_chunks} 
Example #16
Source File: main.py    From DARENet with MIT License 5 votes vote down vote up
def get_p_given_budget(budget):
    assert budget >= FLOPs[0], "budget invalid"
    coeff = []
    for s in range(STAGES):
        coeff.append(FLOPs[STAGES - 1 - s] - budget)
    roots = np.roots(coeff)
    q = roots.real[abs(roots.imag) < 1e-5][0]
    p = []
    for i in range(STAGES):
        p.append(q ** i)
    total = np.sum(p)
    p = [element * 1.0 / total for element in p]
    return p 
Example #17
Source File: fluid_flow_geometry.py    From ross with MIT License 5 votes vote down vote up
def calculate_eccentricity_ratio(modified_s):
    """Calculate the eccentricity ratio for a given load using the sommerfeld number.
    Suitable only for short bearings.
    Parameters
    ----------
    modified_s: float
        The modified sommerfeld number.
    Returns
    -------
    float
        The eccentricity ratio.
    Examples
    --------
    >>> modified_s = 6.329494061103412
    >>> calculate_eccentricity_ratio(modified_s) # doctest: +ELLIPSIS
    0.04999...
    """
    coefficients = [
        1,
        -4,
        (6 - (modified_s ** 2) * (16 - np.pi ** 2)),
        -(4 + (np.pi ** 2) * (modified_s ** 2)),
        1,
    ]
    roots = np.roots(coefficients)
    for i in range(0, len(roots)):
        if 0 <= roots[i] <= 1:
            return np.sqrt(roots[i].real)
    raise ValueError("Eccentricity ratio could not be calculated.") 
Example #18
Source File: questionnaire.py    From reportgen with MIT License 5 votes vote down vote up
def clean_ftime(ftime,cut_percent=0.25):
    '''
    ftime 是完成问卷的秒数
    思路:
    1、只考虑截断问卷完成时间较小的样本
    2、找到完成时间变化的拐点,即需要截断的时间点
    返回:r
    建议截断<r的样本
    '''
    t_min=int(ftime.min())
    t_cut=int(ftime.quantile(cut_percent))
    x=np.array(range(t_min,t_cut))
    y=np.array([len(ftime[ftime<=i]) for i in range(t_min,t_cut)])
    z1 = np.polyfit(x, y, 4) # 拟合得到的函数
    z2=np.polyder(z1,2) #求二阶导数
    r=np.roots(np.polyder(z2,1))
    r=int(r[0])
    return r



## ===========================================================
#
#
#                     数据分析和输出                          #
#
#
## ========================================================== 
Example #19
Source File: partitioners.py    From pyFTS with GNU General Public License v3.0 5 votes vote down vote up
def __init__(self, data, part, **kwargs):
        """"""
        super(PolynomialNonStationaryPartitioner, self).__init__(name=part.name, data=data, npart=part.partitions,
                                                                 func=part.membership_function, names=part.setnames,
                                                                 prefix=part.prefix, transformation=part.transformation,
                                                                 indexer=part.indexer, preprocess=False)

        self.sets = {}

        loc_params, wid_params = self.get_polynomial_perturbations(data, **kwargs)

        if self.ordered_sets is None and self.setnames is not None:
            self.ordered_sets = part.setnames
        else:
            self.ordered_sets = stationary_fs.set_ordered(part.sets)

        for ct, key in enumerate(self.ordered_sets):
            set = part.sets[key]
            loc_roots = np.roots(loc_params[ct])[0]
            wid_roots = np.roots(wid_params[ct])[0]
            tmp = common.FuzzySet(set.name, set.mf, set.parameters,
                           location=perturbation.polynomial,
                           location_params=loc_params[ct],
                           location_roots=loc_roots, #**kwargs)
                           width=perturbation.polynomial,
                           width_params=wid_params[ct],
                           width_roots=wid_roots, **kwargs)

            self.sets[set.name] = tmp 
Example #20
Source File: test_polynomial.py    From ImageFusion with MIT License 5 votes vote down vote up
def test_roots(self):
        assert_array_equal(np.roots([1, 0, 0]), [0, 0]) 
Example #21
Source File: test_polynomial.py    From predictive-maintenance-using-machine-learning with 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 #22
Source File: yellowfin_test.py    From fine-lm with 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 #23
Source File: yellowfin.py    From pytorch-planet-amazon with Apache License 2.0 5 votes vote down vote up
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 #24
Source File: yellowfin_test.py    From YellowFin_Pytorch with 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 #25
Source File: notch_filter.py    From spiketoolkit with MIT License 5 votes vote down vote up
def __init__(self, recording, freq=3000, q=30, chunk_size=30000, cache_chunks=False):
        assert HAVE_NFR, "To use the NotchFilterRecording, install scipy: \n\n pip install scipy\n\n"
        self._freq = freq
        self._q = q
        fn = 0.5 * float(recording.get_sampling_frequency())
        self._b, self._a = ss.iirnotch(self._freq / fn, self._q)

        if not np.all(np.abs(np.roots(self._a)) < 1):
            raise ValueError('Filter is not stable')
        FilterRecording.__init__(self, recording=recording, chunk_size=chunk_size, cache_chunks=cache_chunks)
        self.copy_channel_properties(recording)
        self.is_filtered = self._recording.is_filtered

        self._kwargs = {'recording': recording.make_serialized_dict(), 'freq': freq,
                        'q': q, 'chunk_size': chunk_size, 'cache_chunks': cache_chunks} 
Example #26
Source File: quadrocoptertrajectory.py    From crossgap_il_rl with GNU General Public License v2.0 5 votes vote down vote up
def get_min_max_acc(self, t1, t2):
        """Return the extrema of the acceleration trajectory between t1 and t2."""
        if self._accPeakTimes[0] is None:
            #uninitialised: calculate the roots of the polynomial
            if self._a:
                #solve a quadratic
                det = self._b*self._b - 2*self._g*self._a
                if det<0:
                    #no real roots
                    self._accPeakTimes[0] = 0
                    self._accPeakTimes[1] = 0
                else:
                    self._accPeakTimes[0] = (-self._b + np.sqrt(det))/self._a
                    self._accPeakTimes[1] = (-self._b - np.sqrt(det))/self._a
            else:
                #_g + _b*t == 0:
                if self._b:
                    self._accPeakTimes[0] = -self._g/self._b
                    self._accPeakTimes[1] = 0
                else:
                    self._accPeakTimes[0] = 0
                    self._accPeakTimes[1] = 0

        #Evaluate the acceleration at the boundaries of the period:
        aMinOut = min(self.get_acceleration(t1), self.get_acceleration(t2))
        aMaxOut = max(self.get_acceleration(t1), self.get_acceleration(t2))

        #Evaluate at the maximum/minimum times:
        for i in [0,1]:
            if self._accPeakTimes[i] <= t1: continue
            if self._accPeakTimes[i] >= t2: continue
            
            aMinOut = min(aMinOut, self.get_acceleration(self._accPeakTimes[i]))
            aMaxOut = max(aMaxOut, self.get_acceleration(self._accPeakTimes[i]))
        return (aMinOut, aMaxOut) 
Example #27
Source File: gauss_method.py    From orbitdeterminator with MIT License 5 votes vote down vote up
def read_args():
    parser = argparse.ArgumentParser()
    parser.add_argument('-f', '--file_path', type=str, help="path to IOD-formatted data file", default='../example_data/SATOBS-ML-19200716.txt')
    parser.add_argument('-o', '--obs_array', help="list of lines in file to be read", type=str, default=None)
    parser.add_argument('-b', '--body_name', type=str, help="observed object/body name", default=None)
    parser.add_argument('-r', '--root_index', nargs='*', help="user selection for multiple roots of Gauss polynomial (see docs for more information)", default=None)
    parser.add_argument('-i', '--iterations', type=int, help="number of iterations of Gauss method refinement", default=0)
    parser.add_argument('-p', '--plot', default=True, type=lambda x: (str(x).lower() == 'true'))
    return parser.parse_args() 
Example #28
Source File: common.py    From GraphicDesignPatternByPython with 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 #29
Source File: test_polynomial.py    From GraphicDesignPatternByPython with MIT License 5 votes vote down vote up
def test_roots(self):
        assert_array_equal(np.roots([1, 0, 0]), [0, 0]) 
Example #30
Source File: inversion.py    From oggm with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def _inversion_poly(a3, a0):
    """Solve for degree 5 polynomial with coefficients a5=1, a3, a0."""
    sols = np.roots([1., 0., a3, 0., 0., a0])
    test = (np.isreal(sols)*np.greater(sols, [0]*len(sols)))
    return sols[test][0].real