Python numpy.poly1d() Examples

def remove_linear_BG_XAS_preedge(
    xmcd_data, scanparams, process_parameters=None, process_number=-1
):
    """Should remove a linear bg based on the preedge average"""
    preedge_spectrum = get_preedge_spectrum(process_parameters, xmcd_data)

    preedge_poly = np.poly1d(
        np.polyfit(preedge_spectrum["Energy"], preedge_spectrum["XAS"], 1)
    )

    xas_bg = preedge_poly(xmcd_data["Energy"])

    for xas in ["XAS+", "XAS-", "XAS"]:
        xmcd_data[xas] -= xas_bg

    return (xmcd_data, {"xas_bg_poly_coeffs": " ".join(map(str, preedge_poly.coeffs))}) 

def test_poly1d_math(self):
        # here we use some simple coeffs to make calculations easier
        p = np.poly1d([1., 2, 4])
        q = np.poly1d([4., 2, 1])
        assert_equal(p/q, (np.poly1d([0.25]), np.poly1d([1.5, 3.75])))
        assert_equal(p.integ(), np.poly1d([1/3, 1., 4., 0.]))
        assert_equal(p.integ(1), np.poly1d([1/3, 1., 4., 0.]))

        p = np.poly1d([1., 2, 3])
        q = np.poly1d([3., 2, 1])
        assert_equal(p * q, np.poly1d([3., 8., 14., 8., 3.]))
        assert_equal(p + q, np.poly1d([4., 4., 4.]))
        assert_equal(p - q, np.poly1d([-2., 0., 2.]))
        assert_equal(p ** 4, np.poly1d([1., 8., 36., 104., 214., 312., 324., 216., 81.]))
        assert_equal(p(q), np.poly1d([9., 12., 16., 8., 6.]))
        assert_equal(q(p), np.poly1d([3., 12., 32., 40., 34.]))
        assert_equal(p.deriv(), np.poly1d([2., 2.]))
        assert_equal(p.deriv(2), np.poly1d([2.]))
        assert_equal(np.polydiv(np.poly1d([1, 0, -1]), np.poly1d([1, 1])),
                     (np.poly1d([1., -1.]), np.poly1d([0.]))) 

def test_poly1d_str_and_repr(self):
        p = np.poly1d([1., 2, 3])
        assert_equal(repr(p), 'poly1d([1., 2., 3.])')
        assert_equal(str(p),
                     '   2\n'
                     '1 x + 2 x + 3')

        q = np.poly1d([3., 2, 1])
        assert_equal(repr(q), 'poly1d([3., 2., 1.])')
        assert_equal(str(q),
                     '   2\n'
                     '3 x + 2 x + 1')

        r = np.poly1d([1.89999 + 2j, -3j, -5.12345678, 2 + 1j])
        assert_equal(str(r),
                     '            3      2\n'
                     '(1.9 + 2j) x - 3j x - 5.123 x + (2 + 1j)')

        assert_equal(str(np.poly1d([-3, -2, -1])),
                     '    2\n'
                     '-3 x - 2 x - 1') 

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)) 

def __init__(self, roots, weights=None, hn=1.0, kn=1.0, wfunc=None, limits=None, monic=0,eval_func=None):
        np.poly1d.__init__(self, roots, r=1)
        equiv_weights = [weights[k] / wfunc(roots[k]) for k in range(len(roots))]
        self.__dict__['weights'] = np.array(list(zip(roots,weights,equiv_weights)))
        self.__dict__['weight_func'] = wfunc
        self.__dict__['limits'] = limits
        mu = sqrt(hn)
        if monic:
            evf = eval_func
            if evf:
                eval_func = lambda x: evf(x)/kn
            mu = mu / abs(kn)
            kn = 1.0
        self.__dict__['normcoef'] = mu
        self.__dict__['coeffs'] *= kn

        # Note: eval_func will be discarded on arithmetic
        self.__dict__['_eval_func'] = eval_func 

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 

def __call__(self, v):
        if self._eval_func and not isinstance(v, np.poly1d):
            return self._eval_func(v)
        else:
            return np.poly1d.__call__(self, v) 

def test_sweep_poly_cubic(self):
        p = np.poly1d([2.0, 1.0, 0.0, -2.0])
        t = np.linspace(0, 2.0, 10000)
        phase = waveforms._sweep_poly_phase(t, p)
        tf, f = compute_frequency(t, phase)
        expected = p(tf)
        abserr = np.max(np.abs(f - expected))
        assert_(abserr < 1e-6) 

def test_sweep_poly_quad2(self):
        p = np.poly1d([1.0, 0.0, -2.0])
        t = np.linspace(0, 3.0, 10000)
        phase = waveforms._sweep_poly_phase(t, p)
        tf, f = compute_frequency(t, phase)
        expected = p(tf)
        abserr = np.max(np.abs(f - expected))
        assert_(abserr < 1e-6) 

def test_sweep_poly_const(self):
        p = np.poly1d(2.0)
        t = np.linspace(0, 3.0, 10000)
        phase = waveforms._sweep_poly_phase(t, p)
        tf, f = compute_frequency(t, phase)
        expected = p(tf)
        abserr = np.max(np.abs(f - expected))
        assert_(abserr < 1e-6) 

def test_sweep_poly_quad1(self):
        p = np.poly1d([1.0, 0.0, 1.0])
        t = np.linspace(0, 3.0, 10000)
        phase = waveforms._sweep_poly_phase(t, p)
        tf, f = compute_frequency(t, phase)
        expected = p(tf)
        abserr = np.max(np.abs(f - expected))
        assert_(abserr < 1e-6) 

def test_str_leading_zeros(self):
        p = np.poly1d([4, 3, 2, 1])
        p[3] = 0
        assert_equal(str(p),
                     "   2\n"
                     "3 x + 2 x + 1")

        p = np.poly1d([1, 2])
        p[0] = 0
        p[1] = 0
        assert_equal(str(p), " \n0") 

def test_sweep_poly_cubic2(self):
        """Use an array of coefficients instead of a poly1d."""
        p = np.array([2.0, 1.0, 0.0, -2.0])
        t = np.linspace(0, 2.0, 10000)
        phase = waveforms._sweep_poly_phase(t, p)
        tf, f = compute_frequency(t, phase)
        expected = np.poly1d(p)(tf)
        abserr = np.max(np.abs(f - expected))
        assert_(abserr < 1e-6) 

def check_poly(self, func, cls, param_ranges=[], x_range=[], nn=10,
                   nparam=10, nx=10, rtol=1e-8):
        np.random.seed(1234)

        dataset = []
        for n in np.arange(nn):
            params = [a + (b-a)*np.random.rand(nparam) for a,b in param_ranges]
            params = np.asarray(params).T
            if not param_ranges:
                params = [0]
            for p in params:
                if param_ranges:
                    p = (n,) + tuple(p)
                else:
                    p = (n,)
                x = x_range[0] + (x_range[1] - x_range[0])*np.random.rand(nx)
                x[0] = x_range[0]  # always include domain start point
                x[1] = x_range[1]  # always include domain end point
                poly = np.poly1d(cls(*p))
                z = np.c_[np.tile(p, (nx,1)), x, poly(x)]
                dataset.append(z)

        dataset = np.concatenate(dataset, axis=0)

        def polyfunc(*p):
            p = (p[0].astype(int),) + p[1:]
            return func(*p)

        olderr = np.seterr(all='raise')
        try:
            ds = FuncData(polyfunc, dataset, list(range(len(param_ranges)+2)), -1,
                          rtol=rtol)
            ds.check()
        finally:
            np.seterr(**olderr) 

def test_poly_eq(self):
        # Ticket #554
        x = np.poly1d([1, 2, 3])
        y = np.poly1d([3, 4])
        assert_(x != y)
        assert_(x == x) 

def test_poly_div(self):
        # Ticket #553
        u = np.poly1d([1, 2, 3])
        v = np.poly1d([1, 2, 3, 4, 5])
        q, r = np.polydiv(u, v)
        assert_equal(q*v + r, u) 

def test_poly1d_nan_roots(self):
        # Ticket #396
        p = np.poly1d([np.nan, np.nan, 1], r=0)
        assert_raises(np.linalg.LinAlgError, getattr, p, "r") 

def test_poly1d(self):
        # Ticket #28
        assert_equal(np.poly1d([1]) - np.poly1d([1, 0]),
                     np.poly1d([-1, 1])) 

def test_poly_eq(self):
        p = np.poly1d([1, 2, 3])
        p2 = np.poly1d([1, 2, 4])
        assert_equal(p == None, False)
        assert_equal(p != None, True)
        assert_equal(p == p, True)
        assert_equal(p == p2, False)
        assert_equal(p != p2, True) 

def test_integ_coeffs(self):
        p = np.poly1d([3, 2, 1])
        p2 = p.integ(3, k=[9, 7, 6])
        assert_(
            (p2.coeffs == [1/4./5., 1/3./4., 1/2./3., 9/1./2., 7, 6]).all()) 

def test_complex(self):
        p = np.poly1d([3j, 2j, 1j])
        p2 = p.integ()
        assert_((p2.coeffs == [1j, 1j, 1j, 0]).all())
        p2 = p.deriv()
        assert_((p2.coeffs == [6j, 2j]).all()) 

def test_objects(self):
        from decimal import Decimal
        p = np.poly1d([Decimal('4.0'), Decimal('3.0'), Decimal('2.0')])
        p2 = p * Decimal('1.333333333333333')
        assert_(p2[1] == Decimal("3.9999999999999990"))
        p2 = p.deriv()
        assert_(p2[1] == Decimal('8.0'))
        p2 = p.integ()
        assert_(p2[3] == Decimal("1.333333333333333333333333333"))
        assert_(p2[2] == Decimal('1.5'))
        assert_(np.issubdtype(p2.coeffs.dtype, np.object_))
        p = np.poly([Decimal(1), Decimal(2)])
        assert_equal(np.poly([Decimal(1), Decimal(2)]),
                     [1, Decimal(-3), Decimal(2)]) 

def test_str_leading_zeros(self):
        p = np.poly1d([4, 3, 2, 1])
        p[3] = 0
        assert_equal(str(p),
                     "   2\n"
                     "3 x + 2 x + 1")

        p = np.poly1d([1, 2])
        p[0] = 0
        p[1] = 0
        assert_equal(str(p), " \n0") 

def fit_y(self, X, Y, x1, x2):
        len(X) != 0
        # if X only include one point, the function will get line y=Y[0]
        if np.sum(X == X[0]) == len(X):
            return Y[0], Y[0]
        p = np.poly1d(np.polyfit(X, Y, 1))
        return p(x1), p(x2) 

def fit_y(self, X, Y, x1, x2):
        len(X) != 0
        # if X only include one point, the function will get line y=Y[0]
        if np.sum(X == X[0]) == len(X):
            return Y[0], Y[0]
        p = np.poly1d(np.polyfit(X, Y, 1))
        return p(x1), p(x2) 

def test_polyder_return_type(self):
        # Ticket #1249
        assert_(isinstance(np.polyder(np.poly1d([1]), 0), np.poly1d))
        assert_(isinstance(np.polyder([1], 0), np.ndarray))
        assert_(isinstance(np.polyder(np.poly1d([1]), 1), np.poly1d))
        assert_(isinstance(np.polyder([1], 1), np.ndarray)) 

def test_poly_eq(self, level=rlevel):
        # Ticket #554
        x = np.poly1d([1, 2, 3])
        y = np.poly1d([3, 4])
        assert_(x != y)
        assert_(x == x) 

def test_poly_div(self, level=rlevel):
        # Ticket #553
        u = np.poly1d([1, 2, 3])
        v = np.poly1d([1, 2, 3, 4, 5])
        q, r = np.polydiv(u, v)
        assert_equal(q*v + r, u) 

def test_poly1d(self, level=rlevel):
        # Ticket #28
        assert_equal(np.poly1d([1]) - np.poly1d([1, 0]),
                     np.poly1d([-1, 1])) 

def test_integ_coeffs(self):
        p = np.poly1d([3, 2, 1])
        p2 = p.integ(3, k=[9, 7, 6])
        assert_(
            (p2.coeffs == [1/4./5., 1/3./4., 1/2./3., 9/1./2., 7, 6]).all())