From d54e7351aeda804de6e7b5963bb2b9fa5c76d027 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Fri, 25 Sep 2015 16:29:09 +0100 Subject: ENH: Allow specification of terms to fit in legfit The argument `deg` is enhanced to allow an array_like argument to past which specifies which terms to include in the fit. The returned coef array is exapnded to have entries of 0 for all coefficients which were not included in the fit. --- numpy/polynomial/legendre.py | 47 ++++++++++++++++++++++++++++++++++++-------- 1 file changed, 39 insertions(+), 8 deletions(-) (limited to 'numpy') diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py index c91cb72ec..2035ba6e9 100644 --- a/numpy/polynomial/legendre.py +++ b/numpy/polynomial/legendre.py @@ -1418,8 +1418,14 @@ def legfit(x, y, deg, rcond=None, full=False, w=None): y-coordinates of the sample points. Several data sets of sample points sharing the same x-coordinates can be fitted at once by passing in a 2D-array that contains one dataset per column. - deg : int - Degree of the fitting polynomial + deg : int or array_like + Degree of the fitting polynomial. If `deg` is a single integer + all terms up to and including the `deg`'th term are included. + `deg` may alternatively be a list or array specifying which + terms in the Legendre expansion to include in the fit. + + .. versionchanged:: 1.11.0 + `deg` may be a list specifying which terms to fit rcond : float, optional Relative condition number of the fit. Singular values smaller than this relative to the largest singular value will be ignored. The @@ -1440,9 +1446,11 @@ def legfit(x, y, deg, rcond=None, full=False, w=None): Returns ------- coef : ndarray, shape (M,) or (M, K) - Legendre coefficients ordered from low to high. If `y` was 2-D, - the coefficients for the data in column k of `y` are in column - `k`. + Legendre coefficients ordered from low to high. If `y` was + 2-D, the coefficients for the data in column k of `y` are in + column `k`. If `deg` is specified as a list, coefficients for + terms not included in the fit are set equal to zero in the + returned `coef`. [residuals, rank, singular_values, rcond] : list These values are only returned if `full` = True @@ -1511,12 +1519,14 @@ def legfit(x, y, deg, rcond=None, full=False, w=None): -------- """ - order = int(deg) + 1 x = np.asarray(x) + 0.0 y = np.asarray(y) + 0.0 + deg = np.asarray([deg,], dtype=int).flatten() # check arguments. - if deg < 0: + if deg.size < 1: + raise TypeError("expected deg to be one or more integers") + if deg.min() < 0: raise ValueError("expected deg >= 0") if x.ndim != 1: raise TypeError("expected 1D vector for x") @@ -1527,8 +1537,20 @@ def legfit(x, y, deg, rcond=None, full=False, w=None): if len(x) != len(y): raise TypeError("expected x and y to have same length") + if deg.size == 1: + restricted_fit = False + lmax = deg[0] + order = lmax + 1 + else: + restricted_fit = True + lmax = deg.max() + order = deg.size + # set up the least squares matrices in transposed form - lhs = legvander(x, deg).T + van = legvander(x, lmax) + if restricted_fit: + van = van[:, deg] + lhs = van.T rhs = y.T if w is not None: w = np.asarray(w) + 0.0 @@ -1556,6 +1578,15 @@ def legfit(x, y, deg, rcond=None, full=False, w=None): c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond) c = (c.T/scl).T + # Expand c to include non-fitted coefficients which are set to zero + if restricted_fit: + if c.ndim == 2: + cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype) + else: + cc = np.zeros(lmax+1, dtype=c.dtype) + cc[deg] = c + c = cc + # warn on rank reduction if rank != order and not full: msg = "The fit may be poorly conditioned" -- cgit v1.2.1 From 40951167a51392205427c59a9d8425c36e89df08 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Fri, 25 Sep 2015 23:01:57 +0100 Subject: TST: Add tests for legfit with deg specified as list --- numpy/polynomial/tests/test_legendre.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) (limited to 'numpy') diff --git a/numpy/polynomial/tests/test_legendre.py b/numpy/polynomial/tests/test_legendre.py index 8ac1feb58..9c259d14c 100644 --- a/numpy/polynomial/tests/test_legendre.py +++ b/numpy/polynomial/tests/test_legendre.py @@ -388,6 +388,9 @@ class TestFitting(TestCase): def f(x): return x*(x - 1)*(x - 2) + def f2(x): + return x**4 + x**2 + 1 + # Test exceptions assert_raises(ValueError, leg.legfit, [1], [1], -1) assert_raises(TypeError, leg.legfit, [[1]], [1], 0) @@ -397,6 +400,9 @@ class TestFitting(TestCase): assert_raises(TypeError, leg.legfit, [1], [1, 2], 0) assert_raises(TypeError, leg.legfit, [1], [1], 0, w=[[1]]) assert_raises(TypeError, leg.legfit, [1], [1], 0, w=[1, 1]) + assert_raises(ValueError, leg.legfit, [1], [1], [-1,]) + assert_raises(ValueError, leg.legfit, [1], [1], [2, -1, 6]) + assert_raises(TypeError, leg.legfit, [1], [1], []) # Test fit x = np.linspace(0, 2) @@ -405,13 +411,25 @@ class TestFitting(TestCase): coef3 = leg.legfit(x, y, 3) assert_equal(len(coef3), 4) assert_almost_equal(leg.legval(x, coef3), y) + coef3 = leg.legfit(x, y, [0, 1, 2, 3]) + assert_equal(len(coef3), 4) + assert_almost_equal(leg.legval(x, coef3), y) # coef4 = leg.legfit(x, y, 4) assert_equal(len(coef4), 5) assert_almost_equal(leg.legval(x, coef4), y) + coef4 = leg.legfit(x, y, [0, 1, 2, 3, 4]) + assert_equal(len(coef4), 5) + assert_almost_equal(leg.legval(x, coef4), y) + # check things still work if deg is not in strict increasing + coef4 = leg.legfit(x, y, [2, 3, 4, 1, 0]) + assert_equal(len(coef4), 5) + assert_almost_equal(leg.legval(x, coef4), y) # coef2d = leg.legfit(x, np.array([y, y]).T, 3) assert_almost_equal(coef2d, np.array([coef3, coef3]).T) + coef2d = leg.legfit(x, np.array([y, y]).T, [0, 1, 2, 3]) + assert_almost_equal(coef2d, np.array([coef3, coef3]).T) # test weighting w = np.zeros_like(x) yw = y.copy() @@ -419,13 +437,26 @@ class TestFitting(TestCase): y[0::2] = 0 wcoef3 = leg.legfit(x, yw, 3, w=w) assert_almost_equal(wcoef3, coef3) + wcoef3 = leg.legfit(x, yw, [0, 1, 2, 3], w=w) + assert_almost_equal(wcoef3, coef3) # wcoef2d = leg.legfit(x, np.array([yw, yw]).T, 3, w=w) assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T) + wcoef2d = leg.legfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w) + assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T) # test scaling with complex values x points whose square # is zero when summed. x = [1, 1j, -1, -1j] assert_almost_equal(leg.legfit(x, x, 1), [0, 1]) + assert_almost_equal(leg.legfit(x, x, [0, 1]), [0, 1]) + # test fitting only even Legendre polynomials + x = np.linspace(-1, 1) + y = f2(x) + coef1 = leg.legfit(x, y, 4) + assert_almost_equal(leg.legval(x, coef1), y) + coef2 = leg.legfit(x, y, [0, 2, 4]) + assert_almost_equal(leg.legval(x, coef2), y) + assert_almost_equal(coef1, coef2) class TestCompanion(TestCase): -- cgit v1.2.1 From 942f294c06b0285ea3cf2bf223a63700a1ed50f5 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Mon, 12 Oct 2015 20:04:36 +0100 Subject: ENH: Allow specification of terms to fit in chebfit The argument `deg` is enhanced to allow an array_like argument to past which specifies which terms to include in the fit. The returned coef array is exapnded to have entries of 0 for all coefficients which were not included in the fit. --- numpy/polynomial/chebyshev.py | 39 ++++++++++++++++++++++++++++++++++----- 1 file changed, 34 insertions(+), 5 deletions(-) (limited to 'numpy') diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py index de5cbb734..9db613b78 100644 --- a/numpy/polynomial/chebyshev.py +++ b/numpy/polynomial/chebyshev.py @@ -1617,8 +1617,14 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None): y-coordinates of the sample points. Several data sets of sample points sharing the same x-coordinates can be fitted at once by passing in a 2D-array that contains one dataset per column. - deg : int - Degree of the fitting series + deg : int or array_like + Degree of the fitting series. If `deg` is a single integer + all terms up to and including the `deg`'th term are included. + `deg` may alternatively be a list or array specifying which + terms in the Legendre expansion to include in the fit. + + .. versionchanged:: 1.11.0 + `deg` may be a list specifying which terms to fit rcond : float, optional Relative condition number of the fit. Singular values smaller than this relative to the largest singular value will be ignored. The @@ -1710,12 +1716,14 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None): -------- """ - order = int(deg) + 1 x = np.asarray(x) + 0.0 y = np.asarray(y) + 0.0 + deg = np.asarray([deg,], dtype=int).flatten() # check arguments. - if deg < 0: + if deg.size < 1: + raise TypeError("expected deg to be one or more integers") + if deg.min() < 0: raise ValueError("expected deg >= 0") if x.ndim != 1: raise TypeError("expected 1D vector for x") @@ -1726,8 +1734,20 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None): if len(x) != len(y): raise TypeError("expected x and y to have same length") + if deg.size == 1: + restricted_fit = False + lmax = deg[0] + order = lmax + 1 + else: + restricted_fit = True + lmax = deg.max() + order = deg.size + # set up the least squares matrices in transposed form - lhs = chebvander(x, deg).T + van = chebvander(x, lmax) + if restricted_fit: + van = van[:, deg] + lhs = van.T rhs = y.T if w is not None: w = np.asarray(w) + 0.0 @@ -1755,6 +1775,15 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None): c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond) c = (c.T/scl).T + # Expand c to include non-fitted coefficients which are set to zero + if restricted_fit: + if c.ndim == 2: + cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype) + else: + cc = np.zeros(lmax+1, dtype=c.dtype) + cc[deg] = c + c = cc + # warn on rank reduction if rank != order and not full: msg = "The fit may be poorly conditioned" -- cgit v1.2.1 From 1a9fb061bb4f217f335616c65abd36644c2f2ac7 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Mon, 12 Oct 2015 20:23:56 +0100 Subject: TST: Add tests for chebfit with deg specified as list --- numpy/polynomial/tests/test_chebyshev.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) (limited to 'numpy') diff --git a/numpy/polynomial/tests/test_chebyshev.py b/numpy/polynomial/tests/test_chebyshev.py index a596905f6..8d992c4f0 100644 --- a/numpy/polynomial/tests/test_chebyshev.py +++ b/numpy/polynomial/tests/test_chebyshev.py @@ -399,6 +399,9 @@ class TestFitting(TestCase): def f(x): return x*(x - 1)*(x - 2) + def f2(x): + return x**4 + x**2 + 1 + # Test exceptions assert_raises(ValueError, cheb.chebfit, [1], [1], -1) assert_raises(TypeError, cheb.chebfit, [[1]], [1], 0) @@ -408,6 +411,9 @@ class TestFitting(TestCase): assert_raises(TypeError, cheb.chebfit, [1], [1, 2], 0) assert_raises(TypeError, cheb.chebfit, [1], [1], 0, w=[[1]]) assert_raises(TypeError, cheb.chebfit, [1], [1], 0, w=[1, 1]) + assert_raises(ValueError, cheb.chebfit, [1], [1], [-1,]) + assert_raises(ValueError, cheb.chebfit, [1], [1], [2, -1, 6]) + assert_raises(TypeError, cheb.chebfit, [1], [1], []) # Test fit x = np.linspace(0, 2) @@ -416,13 +422,25 @@ class TestFitting(TestCase): coef3 = cheb.chebfit(x, y, 3) assert_equal(len(coef3), 4) assert_almost_equal(cheb.chebval(x, coef3), y) + coef3 = cheb.chebfit(x, y, [0, 1, 2, 3]) + assert_equal(len(coef3), 4) + assert_almost_equal(cheb.chebval(x, coef3), y) # coef4 = cheb.chebfit(x, y, 4) assert_equal(len(coef4), 5) assert_almost_equal(cheb.chebval(x, coef4), y) + coef4 = cheb.chebfit(x, y, [0, 1, 2, 3, 4]) + assert_equal(len(coef4), 5) + assert_almost_equal(cheb.chebval(x, coef4), y) + # check things still work if deg is not in strict increasing + coef4 = cheb.chebfit(x, y, [2, 3, 4, 1, 0]) + assert_equal(len(coef4), 5) + assert_almost_equal(cheb.chebval(x, coef4), y) # coef2d = cheb.chebfit(x, np.array([y, y]).T, 3) assert_almost_equal(coef2d, np.array([coef3, coef3]).T) + coef2d = cheb.chebfit(x, np.array([y, y]).T, [0, 1, 2, 3]) + assert_almost_equal(coef2d, np.array([coef3, coef3]).T) # test weighting w = np.zeros_like(x) yw = y.copy() @@ -430,13 +448,26 @@ class TestFitting(TestCase): y[0::2] = 0 wcoef3 = cheb.chebfit(x, yw, 3, w=w) assert_almost_equal(wcoef3, coef3) + wcoef3 = cheb.chebfit(x, yw, [0, 1, 2, 3], w=w) + assert_almost_equal(wcoef3, coef3) # wcoef2d = cheb.chebfit(x, np.array([yw, yw]).T, 3, w=w) assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T) + wcoef2d = cheb.chebfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w) + assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T) # test scaling with complex values x points whose square # is zero when summed. x = [1, 1j, -1, -1j] assert_almost_equal(cheb.chebfit(x, x, 1), [0, 1]) + assert_almost_equal(cheb.chebfit(x, x, [0, 1]), [0, 1]) + # test fitting only even polynomials + x = np.linspace(-1, 1) + y = f2(x) + coef1 = cheb.chebfit(x, y, 4) + assert_almost_equal(cheb.chebval(x, coef1), y) + coef2 = cheb.chebfit(x, y, [0, 2, 4]) + assert_almost_equal(cheb.chebval(x, coef2), y) + assert_almost_equal(coef1, coef2) class TestCompanion(TestCase): -- cgit v1.2.1 From acc294dcd7d3998f8f2b82cad8f3ed6db48c1f00 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Mon, 12 Oct 2015 21:03:09 +0100 Subject: ENH: Allow specification of terms to fit in hermfit The argument `deg` is enhanced to allow an array_like argument to past which specifies which terms to include in the fit. The returned coef array is exapnded to have entries of 0 for all coefficients which were not included in the fit. --- numpy/polynomial/hermite.py | 39 ++++++++++++++++++++++++++++++++++----- 1 file changed, 34 insertions(+), 5 deletions(-) (limited to 'numpy') diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py index 2ce1d97a8..5d4b357fe 100644 --- a/numpy/polynomial/hermite.py +++ b/numpy/polynomial/hermite.py @@ -1388,8 +1388,14 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None): y-coordinates of the sample points. Several data sets of sample points sharing the same x-coordinates can be fitted at once by passing in a 2D-array that contains one dataset per column. - deg : int - Degree of the fitting polynomial + deg : int or array_like + Degree of the fitting polynomial. If `deg` is a single integer + all terms up to and including the `deg`'th term are included. + `deg` may alternatively be a list or array specifying which + terms in the Legendre expansion to include in the fit. + + .. versionchanged:: 1.11.0 + `deg` may be a list specifying which terms to fit rcond : float, optional Relative condition number of the fit. Singular values smaller than this relative to the largest singular value will be ignored. The @@ -1486,12 +1492,14 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None): array([ 0.97902637, 1.99849131, 3.00006 ]) """ - order = int(deg) + 1 x = np.asarray(x) + 0.0 y = np.asarray(y) + 0.0 + deg = np.asarray([deg,], dtype=int).flatten() # check arguments. - if deg < 0: + if deg.size < 1: + raise TypeError("expected deg to be one or more integers") + if deg.min() < 0: raise ValueError("expected deg >= 0") if x.ndim != 1: raise TypeError("expected 1D vector for x") @@ -1502,8 +1510,20 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None): if len(x) != len(y): raise TypeError("expected x and y to have same length") + if deg.size == 1: + restricted_fit = False + lmax = deg[0] + order = lmax + 1 + else: + restricted_fit = True + lmax = deg.max() + order = deg.size + # set up the least squares matrices in transposed form - lhs = hermvander(x, deg).T + van = hermvander(x, lmax) + if restricted_fit: + van = van[:, deg] + lhs = van.T rhs = y.T if w is not None: w = np.asarray(w) + 0.0 @@ -1531,6 +1551,15 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None): c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond) c = (c.T/scl).T + # Expand c to include non-fitted coefficients which are set to zero + if restricted_fit: + if c.ndim == 2: + cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype) + else: + cc = np.zeros(lmax+1, dtype=c.dtype) + cc[deg] = c + c = cc + # warn on rank reduction if rank != order and not full: msg = "The fit may be poorly conditioned" -- cgit v1.2.1 From b8180bbf97505a544324c90f407cdf2c5913c612 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Mon, 12 Oct 2015 21:10:22 +0100 Subject: TST: Add tests for hermfit with deg specified as list --- numpy/polynomial/tests/test_hermite.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) (limited to 'numpy') diff --git a/numpy/polynomial/tests/test_hermite.py b/numpy/polynomial/tests/test_hermite.py index e67625a88..04da72b26 100644 --- a/numpy/polynomial/tests/test_hermite.py +++ b/numpy/polynomial/tests/test_hermite.py @@ -387,6 +387,9 @@ class TestFitting(TestCase): def f(x): return x*(x - 1)*(x - 2) + def f2(x): + return x**4 + x**2 + 1 + # Test exceptions assert_raises(ValueError, herm.hermfit, [1], [1], -1) assert_raises(TypeError, herm.hermfit, [[1]], [1], 0) @@ -396,6 +399,9 @@ class TestFitting(TestCase): assert_raises(TypeError, herm.hermfit, [1], [1, 2], 0) assert_raises(TypeError, herm.hermfit, [1], [1], 0, w=[[1]]) assert_raises(TypeError, herm.hermfit, [1], [1], 0, w=[1, 1]) + assert_raises(ValueError, herm.hermfit, [1], [1], [-1,]) + assert_raises(ValueError, herm.hermfit, [1], [1], [2, -1, 6]) + assert_raises(TypeError, herm.hermfit, [1], [1], []) # Test fit x = np.linspace(0, 2) @@ -404,13 +410,25 @@ class TestFitting(TestCase): coef3 = herm.hermfit(x, y, 3) assert_equal(len(coef3), 4) assert_almost_equal(herm.hermval(x, coef3), y) + coef3 = herm.hermfit(x, y, [0, 1, 2, 3]) + assert_equal(len(coef3), 4) + assert_almost_equal(herm.hermval(x, coef3), y) # coef4 = herm.hermfit(x, y, 4) assert_equal(len(coef4), 5) assert_almost_equal(herm.hermval(x, coef4), y) + coef4 = herm.hermfit(x, y, [0, 1, 2, 3, 4]) + assert_equal(len(coef4), 5) + assert_almost_equal(herm.hermval(x, coef4), y) + # check things still work if deg is not in strict increasing + coef4 = herm.hermfit(x, y, [2, 3, 4, 1, 0]) + assert_equal(len(coef4), 5) + assert_almost_equal(herm.hermval(x, coef4), y) # coef2d = herm.hermfit(x, np.array([y, y]).T, 3) assert_almost_equal(coef2d, np.array([coef3, coef3]).T) + coef2d = herm.hermfit(x, np.array([y, y]).T, [0, 1, 2, 3]) + assert_almost_equal(coef2d, np.array([coef3, coef3]).T) # test weighting w = np.zeros_like(x) yw = y.copy() @@ -418,13 +436,26 @@ class TestFitting(TestCase): y[0::2] = 0 wcoef3 = herm.hermfit(x, yw, 3, w=w) assert_almost_equal(wcoef3, coef3) + wcoef3 = herm.hermfit(x, yw, [0, 1, 2, 3], w=w) + assert_almost_equal(wcoef3, coef3) # wcoef2d = herm.hermfit(x, np.array([yw, yw]).T, 3, w=w) assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T) + wcoef2d = herm.hermfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w) + assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T) # test scaling with complex values x points whose square # is zero when summed. x = [1, 1j, -1, -1j] assert_almost_equal(herm.hermfit(x, x, 1), [0, .5]) + assert_almost_equal(herm.hermfit(x, x, [0, 1]), [0, .5]) + # test fitting only even Legendre polynomials + x = np.linspace(-1, 1) + y = f2(x) + coef1 = herm.hermfit(x, y, 4) + assert_almost_equal(herm.hermval(x, coef1), y) + coef2 = herm.hermfit(x, y, [0, 2, 4]) + assert_almost_equal(herm.hermval(x, coef2), y) + assert_almost_equal(coef1, coef2) class TestCompanion(TestCase): -- cgit v1.2.1 From 8dcf03358a97bf3a577025c03438d381c8b4cb72 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Thu, 5 Nov 2015 15:03:27 +0000 Subject: ENH: Allow specification of terms to fit in lagfit The argument `deg` is enhanced to allow an array_like argument to past which specifies which terms to include in the fit. The returned coef array is exapnded to have entries of 0 for all coefficients which were not included in the fit. --- numpy/polynomial/laguerre.py | 39 ++++++++++++++++++++++++++++++++++----- 1 file changed, 34 insertions(+), 5 deletions(-) (limited to 'numpy') diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py index fffe9e6b6..280e28159 100644 --- a/numpy/polynomial/laguerre.py +++ b/numpy/polynomial/laguerre.py @@ -1387,8 +1387,14 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None): y-coordinates of the sample points. Several data sets of sample points sharing the same x-coordinates can be fitted at once by passing in a 2D-array that contains one dataset per column. - deg : int - Degree of the fitting polynomial + deg : int or array_like + Degree of the fitting polynomial. If `deg` is a single integer + all terms up to and including the `deg`'th term are included. + `deg` may alternatively be a list or array specifying which + terms in the Legendre expansion to include in the fit. + + .. versionchanged:: 1.11.0 + `deg` may be a list specifying which terms to fit rcond : float, optional Relative condition number of the fit. Singular values smaller than this relative to the largest singular value will be ignored. The @@ -1485,12 +1491,14 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None): array([ 0.96971004, 2.00193749, 3.00288744]) """ - order = int(deg) + 1 x = np.asarray(x) + 0.0 y = np.asarray(y) + 0.0 + deg = np.asarray([deg,], dtype=int).flatten() # check arguments. - if deg < 0: + if deg.size < 1: + raise TypeError("expected deg to be one or more integers") + if deg.min() < 0: raise ValueError("expected deg >= 0") if x.ndim != 1: raise TypeError("expected 1D vector for x") @@ -1501,8 +1509,20 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None): if len(x) != len(y): raise TypeError("expected x and y to have same length") + if deg.size == 1: + restricted_fit = False + lmax = deg[0] + order = lmax + 1 + else: + restricted_fit = True + lmax = deg.max() + order = deg.size + # set up the least squares matrices in transposed form - lhs = lagvander(x, deg).T + van = lagvander(x, lmax) + if restricted_fit: + van = van[:, deg] + lhs = van.T rhs = y.T if w is not None: w = np.asarray(w) + 0.0 @@ -1530,6 +1550,15 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None): c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond) c = (c.T/scl).T + # Expand c to include non-fitted coefficients which are set to zero + if restricted_fit: + if c.ndim == 2: + cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype) + else: + cc = np.zeros(lmax+1, dtype=c.dtype) + cc[deg] = c + c = cc + # warn on rank reduction if rank != order and not full: msg = "The fit may be poorly conditioned" -- cgit v1.2.1 From 84e0b6ec1b686b908eb2d8bba38337a13f02cb82 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Thu, 5 Nov 2015 15:40:25 +0000 Subject: TST: Add tests for lagfit with deg specified as list --- numpy/polynomial/tests/test_laguerre.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) (limited to 'numpy') diff --git a/numpy/polynomial/tests/test_laguerre.py b/numpy/polynomial/tests/test_laguerre.py index 1dc57a960..c25476088 100644 --- a/numpy/polynomial/tests/test_laguerre.py +++ b/numpy/polynomial/tests/test_laguerre.py @@ -393,6 +393,9 @@ class TestFitting(TestCase): assert_raises(TypeError, lag.lagfit, [1], [1, 2], 0) assert_raises(TypeError, lag.lagfit, [1], [1], 0, w=[[1]]) assert_raises(TypeError, lag.lagfit, [1], [1], 0, w=[1, 1]) + assert_raises(ValueError, lag.lagfit, [1], [1], [-1,]) + assert_raises(ValueError, lag.lagfit, [1], [1], [2, -1, 6]) + assert_raises(TypeError, lag.lagfit, [1], [1], []) # Test fit x = np.linspace(0, 2) @@ -401,13 +404,21 @@ class TestFitting(TestCase): coef3 = lag.lagfit(x, y, 3) assert_equal(len(coef3), 4) assert_almost_equal(lag.lagval(x, coef3), y) + coef3 = lag.lagfit(x, y, [0, 1, 2, 3]) + assert_equal(len(coef3), 4) + assert_almost_equal(lag.lagval(x, coef3), y) # coef4 = lag.lagfit(x, y, 4) assert_equal(len(coef4), 5) assert_almost_equal(lag.lagval(x, coef4), y) + coef4 = lag.lagfit(x, y, [0, 1, 2, 3, 4]) + assert_equal(len(coef4), 5) + assert_almost_equal(lag.lagval(x, coef4), y) # coef2d = lag.lagfit(x, np.array([y, y]).T, 3) assert_almost_equal(coef2d, np.array([coef3, coef3]).T) + coef2d = lag.lagfit(x, np.array([y, y]).T, [0, 1, 2, 3]) + assert_almost_equal(coef2d, np.array([coef3, coef3]).T) # test weighting w = np.zeros_like(x) yw = y.copy() @@ -415,13 +426,18 @@ class TestFitting(TestCase): y[0::2] = 0 wcoef3 = lag.lagfit(x, yw, 3, w=w) assert_almost_equal(wcoef3, coef3) + wcoef3 = lag.lagfit(x, yw, [0, 1, 2, 3], w=w) + assert_almost_equal(wcoef3, coef3) # wcoef2d = lag.lagfit(x, np.array([yw, yw]).T, 3, w=w) assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T) + wcoef2d = lag.lagfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w) + assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T) # test scaling with complex values x points whose square # is zero when summed. x = [1, 1j, -1, -1j] assert_almost_equal(lag.lagfit(x, x, 1), [1, -1]) + assert_almost_equal(lag.lagfit(x, x, [0, 1]), [1, -1]) class TestCompanion(TestCase): -- cgit v1.2.1 From b904ef19b55796ea6d0e43d00a551fc841833b78 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Wed, 2 Dec 2015 20:27:25 +0000 Subject: ENH: Allow specification of terms to fit in polyfit The argument `deg` is enhanced to allow an array_like argument to past which specifies which terms to include in the fit. The returned coef array is exapnded to have entries of 0 for all coefficients which were not included in the fit. --- numpy/polynomial/polynomial.py | 39 ++++++++++++++++++++++++++++++++++----- 1 file changed, 34 insertions(+), 5 deletions(-) (limited to 'numpy') diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py index 60e339a1d..7c922c11b 100644 --- a/numpy/polynomial/polynomial.py +++ b/numpy/polynomial/polynomial.py @@ -1217,8 +1217,14 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None): sharing the same x-coordinates can be (independently) fit with one call to `polyfit` by passing in for `y` a 2-D array that contains one data set per column. - deg : int - Degree of the polynomial(s) to be fit. + deg : int or array_like + Degree of the fitting polynomial. If `deg` is a single integer + all terms up to and including the `deg`'th term are included. + `deg` may alternatively be a list or array specifying which + terms in the Legendre expansion to include in the fit. + + .. versionchanged:: 1.11.0 + `deg` may be a list specifying which terms to fit rcond : float, optional Relative condition number of the fit. Singular values smaller than `rcond`, relative to the largest singular value, will be @@ -1332,12 +1338,14 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None): 0.50443316, 0.28853036]), 1.1324274851176597e-014] """ - order = int(deg) + 1 x = np.asarray(x) + 0.0 y = np.asarray(y) + 0.0 + deg = np.asarray([deg,], dtype=int).flatten() # check arguments. - if deg < 0: + if deg.size < 1: + raise TypeError("expected deg to be one or more integers") + if deg.min() < 0: raise ValueError("expected deg >= 0") if x.ndim != 1: raise TypeError("expected 1D vector for x") @@ -1348,8 +1356,20 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None): if len(x) != len(y): raise TypeError("expected x and y to have same length") + if deg.size == 1: + restricted_fit = False + lmax = deg[0] + order = lmax + 1 + else: + restricted_fit = True + lmax = deg.max() + order = deg.size + # set up the least squares matrices in transposed form - lhs = polyvander(x, deg).T + van = polyvander(x, lmax) + if restricted_fit: + van = van[:, deg] + lhs = van.T rhs = y.T if w is not None: w = np.asarray(w) + 0.0 @@ -1377,6 +1397,15 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None): c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond) c = (c.T/scl).T + # Expand c to include non-fitted coefficients which are set to zero + if restricted_fit: + if c.ndim == 2: + cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype) + else: + cc = np.zeros(lmax+1, dtype=c.dtype) + cc[deg] = c + c = cc + # warn on rank reduction if rank != order and not full: msg = "The fit may be poorly conditioned" -- cgit v1.2.1 From fd5d1a4893d1a4b04d6df94ff83e09a6e4f12df6 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Thu, 3 Dec 2015 14:11:59 +0000 Subject: TST: Add tests for polyfit with deg specified as list --- numpy/polynomial/tests/test_polynomial.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) (limited to 'numpy') diff --git a/numpy/polynomial/tests/test_polynomial.py b/numpy/polynomial/tests/test_polynomial.py index c806a8497..00a52ebce 100644 --- a/numpy/polynomial/tests/test_polynomial.py +++ b/numpy/polynomial/tests/test_polynomial.py @@ -419,6 +419,9 @@ class TestMisc(TestCase): def f(x): return x*(x - 1)*(x - 2) + def f2(x): + return x**4 + x**2 + 1 + # Test exceptions assert_raises(ValueError, poly.polyfit, [1], [1], -1) assert_raises(TypeError, poly.polyfit, [[1]], [1], 0) @@ -428,6 +431,9 @@ class TestMisc(TestCase): assert_raises(TypeError, poly.polyfit, [1], [1, 2], 0) assert_raises(TypeError, poly.polyfit, [1], [1], 0, w=[[1]]) assert_raises(TypeError, poly.polyfit, [1], [1], 0, w=[1, 1]) + assert_raises(ValueError, poly.polyfit, [1], [1], [-1,]) + assert_raises(ValueError, poly.polyfit, [1], [1], [2, -1, 6]) + assert_raises(TypeError, poly.polyfit, [1], [1], []) # Test fit x = np.linspace(0, 2) @@ -436,13 +442,21 @@ class TestMisc(TestCase): coef3 = poly.polyfit(x, y, 3) assert_equal(len(coef3), 4) assert_almost_equal(poly.polyval(x, coef3), y) + coef3 = poly.polyfit(x, y, [0, 1, 2, 3]) + assert_equal(len(coef3), 4) + assert_almost_equal(poly.polyval(x, coef3), y) # coef4 = poly.polyfit(x, y, 4) assert_equal(len(coef4), 5) assert_almost_equal(poly.polyval(x, coef4), y) + coef4 = poly.polyfit(x, y, [0, 1, 2, 3, 4]) + assert_equal(len(coef4), 5) + assert_almost_equal(poly.polyval(x, coef4), y) # coef2d = poly.polyfit(x, np.array([y, y]).T, 3) assert_almost_equal(coef2d, np.array([coef3, coef3]).T) + coef2d = poly.polyfit(x, np.array([y, y]).T, [0, 1, 2, 3]) + assert_almost_equal(coef2d, np.array([coef3, coef3]).T) # test weighting w = np.zeros_like(x) yw = y.copy() @@ -450,13 +464,26 @@ class TestMisc(TestCase): yw[0::2] = 0 wcoef3 = poly.polyfit(x, yw, 3, w=w) assert_almost_equal(wcoef3, coef3) + wcoef3 = poly.polyfit(x, yw, [0, 1, 2, 3], w=w) + assert_almost_equal(wcoef3, coef3) # wcoef2d = poly.polyfit(x, np.array([yw, yw]).T, 3, w=w) assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T) + wcoef2d = poly.polyfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w) + assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T) # test scaling with complex values x points whose square # is zero when summed. x = [1, 1j, -1, -1j] assert_almost_equal(poly.polyfit(x, x, 1), [0, 1]) + assert_almost_equal(poly.polyfit(x, x, [0, 1]), [0, 1]) + # test fitting only even Polyendre polynomials + x = np.linspace(-1, 1) + y = f2(x) + coef1 = poly.polyfit(x, y, 4) + assert_almost_equal(poly.polyval(x, coef1), y) + coef2 = poly.polyfit(x, y, [0, 2, 4]) + assert_almost_equal(poly.polyval(x, coef2), y) + assert_almost_equal(coef1, coef2) def test_polytrim(self): coef = [2, -1, 1, 0] -- cgit v1.2.1 From 6411ec505ed7b36c9768accd338a96f5a64eba93 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Thu, 3 Dec 2015 15:30:16 +0000 Subject: ENH: Allow specification of terms to fit in hermefit The argument `deg` is enhanced to allow an array_like argument to past which specifies which terms to include in the fit. The returned coef array is exapnded to have entries of 0 for all coefficients which were not included in the fit. --- numpy/polynomial/hermite_e.py | 39 ++++++++++++++++++++++++++++++++++----- 1 file changed, 34 insertions(+), 5 deletions(-) (limited to 'numpy') diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py index 394f256a5..da441af83 100644 --- a/numpy/polynomial/hermite_e.py +++ b/numpy/polynomial/hermite_e.py @@ -1385,8 +1385,14 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None): y-coordinates of the sample points. Several data sets of sample points sharing the same x-coordinates can be fitted at once by passing in a 2D-array that contains one dataset per column. - deg : int - Degree of the fitting polynomial + deg : int or array_like + Degree of the fitting polynomial. If `deg` is a single integer + all terms up to and including the `deg`'th term are included. + `deg` may alternatively be a list or array specifying which + terms in the Legendre expansion to include in the fit. + + .. versionchanged:: 1.11.0 + `deg` may be a list specifying which terms to fit rcond : float, optional Relative condition number of the fit. Singular values smaller than this relative to the largest singular value will be ignored. The @@ -1483,12 +1489,14 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None): array([ 1.01690445, 1.99951418, 2.99948696]) """ - order = int(deg) + 1 x = np.asarray(x) + 0.0 y = np.asarray(y) + 0.0 + deg = np.asarray([deg,], dtype=int).flatten() # check arguments. - if deg < 0: + if deg.size < 1: + raise TypeError("expected deg to be one or more integers") + if deg.min() < 0: raise ValueError("expected deg >= 0") if x.ndim != 1: raise TypeError("expected 1D vector for x") @@ -1499,8 +1507,20 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None): if len(x) != len(y): raise TypeError("expected x and y to have same length") + if deg.size == 1: + restricted_fit = False + lmax = deg[0] + order = lmax + 1 + else: + restricted_fit = True + lmax = deg.max() + order = deg.size + # set up the least squares matrices in transposed form - lhs = hermevander(x, deg).T + van = hermevander(x, lmax) + if restricted_fit: + van = van[:, deg] + lhs = van.T rhs = y.T if w is not None: w = np.asarray(w) + 0.0 @@ -1528,6 +1548,15 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None): c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond) c = (c.T/scl).T + # Expand c to include non-fitted coefficients which are set to zero + if restricted_fit: + if c.ndim == 2: + cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype) + else: + cc = np.zeros(lmax+1, dtype=c.dtype) + cc[deg] = c + c = cc + # warn on rank reduction if rank != order and not full: msg = "The fit may be poorly conditioned" -- cgit v1.2.1 From c65093c96746fd5513a648522e105d7df8e8c912 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Thu, 3 Dec 2015 15:45:38 +0000 Subject: TST: Add tests for hermefit with deg specified as list --- numpy/polynomial/tests/test_hermite_e.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) (limited to 'numpy') diff --git a/numpy/polynomial/tests/test_hermite_e.py b/numpy/polynomial/tests/test_hermite_e.py index f8601a828..1162502dc 100644 --- a/numpy/polynomial/tests/test_hermite_e.py +++ b/numpy/polynomial/tests/test_hermite_e.py @@ -388,6 +388,9 @@ class TestFitting(TestCase): def f(x): return x*(x - 1)*(x - 2) + def f2(x): + return x**4 + x**2 + 1 + # Test exceptions assert_raises(ValueError, herme.hermefit, [1], [1], -1) assert_raises(TypeError, herme.hermefit, [[1]], [1], 0) @@ -397,6 +400,9 @@ class TestFitting(TestCase): assert_raises(TypeError, herme.hermefit, [1], [1, 2], 0) assert_raises(TypeError, herme.hermefit, [1], [1], 0, w=[[1]]) assert_raises(TypeError, herme.hermefit, [1], [1], 0, w=[1, 1]) + assert_raises(ValueError, herme.hermefit, [1], [1], [-1,]) + assert_raises(ValueError, herme.hermefit, [1], [1], [2, -1, 6]) + assert_raises(TypeError, herme.hermefit, [1], [1], []) # Test fit x = np.linspace(0, 2) @@ -405,13 +411,25 @@ class TestFitting(TestCase): coef3 = herme.hermefit(x, y, 3) assert_equal(len(coef3), 4) assert_almost_equal(herme.hermeval(x, coef3), y) + coef3 = herme.hermefit(x, y, [0, 1, 2, 3]) + assert_equal(len(coef3), 4) + assert_almost_equal(herme.hermeval(x, coef3), y) # coef4 = herme.hermefit(x, y, 4) assert_equal(len(coef4), 5) assert_almost_equal(herme.hermeval(x, coef4), y) + coef4 = herme.hermefit(x, y, [0, 1, 2, 3, 4]) + assert_equal(len(coef4), 5) + assert_almost_equal(herme.hermeval(x, coef4), y) + # check things still work if deg is not in strict increasing + coef4 = herme.hermefit(x, y, [2, 3, 4, 1, 0]) + assert_equal(len(coef4), 5) + assert_almost_equal(herme.hermeval(x, coef4), y) # coef2d = herme.hermefit(x, np.array([y, y]).T, 3) assert_almost_equal(coef2d, np.array([coef3, coef3]).T) + coef2d = herme.hermefit(x, np.array([y, y]).T, [0, 1, 2, 3]) + assert_almost_equal(coef2d, np.array([coef3, coef3]).T) # test weighting w = np.zeros_like(x) yw = y.copy() @@ -419,13 +437,26 @@ class TestFitting(TestCase): y[0::2] = 0 wcoef3 = herme.hermefit(x, yw, 3, w=w) assert_almost_equal(wcoef3, coef3) + wcoef3 = herme.hermefit(x, yw, [0, 1, 2, 3], w=w) + assert_almost_equal(wcoef3, coef3) # wcoef2d = herme.hermefit(x, np.array([yw, yw]).T, 3, w=w) assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T) + wcoef2d = herme.hermefit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w) + assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T) # test scaling with complex values x points whose square # is zero when summed. x = [1, 1j, -1, -1j] assert_almost_equal(herme.hermefit(x, x, 1), [0, 1]) + assert_almost_equal(herme.hermefit(x, x, [0, 1]), [0, 1]) + # test fitting only even Legendre polynomials + x = np.linspace(-1, 1) + y = f2(x) + coef1 = herme.hermefit(x, y, 4) + assert_almost_equal(herme.hermeval(x, coef1), y) + coef2 = herme.hermefit(x, y, [0, 2, 4]) + assert_almost_equal(herme.hermeval(x, coef2), y) + assert_almost_equal(coef1, coef2) class TestCompanion(TestCase): -- cgit v1.2.1 From 5f7b1af4e652c2ab631634d9778a3a4015e41ced Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Thu, 3 Dec 2015 14:41:00 +0000 Subject: ENH: Allow specification of terms to fit in fit method --- numpy/polynomial/_polybase.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) (limited to 'numpy') diff --git a/numpy/polynomial/_polybase.py b/numpy/polynomial/_polybase.py index 41bdf0484..37eb59f03 100644 --- a/numpy/polynomial/_polybase.py +++ b/numpy/polynomial/_polybase.py @@ -739,8 +739,14 @@ class ABCPolyBase(object): y-coordinates of the sample points. Several data sets of sample points sharing the same x-coordinates can be fitted at once by passing in a 2D-array that contains one dataset per column. - deg : int - Degree of the fitting polynomial. + deg : int or array_like + Degree of the fitting polynomial. If `deg` is a single integer + all terms up to and including the `deg`'th term are included. + `deg` may alternatively be a list or array specifying which + terms in the Legendre expansion to include in the fit. + + .. versionchanged:: 1.11.0 + `deg` may be a list specifying which terms to fit domain : {None, [beg, end], []}, optional Domain to use for the returned series. If ``None``, then a minimal domain that covers the points `x` is chosen. If -- cgit v1.2.1 From 4dd71a3ab7bd019e36998e7c5f98ec2345539f18 Mon Sep 17 00:00:00 2001 From: Jonathan Underwood Date: Thu, 3 Dec 2015 14:49:31 +0000 Subject: TST: Add tests for check_fit with deg specified as list --- numpy/polynomial/tests/test_classes.py | 9 +++++++++ 1 file changed, 9 insertions(+) (limited to 'numpy') diff --git a/numpy/polynomial/tests/test_classes.py b/numpy/polynomial/tests/test_classes.py index cd5a54687..a7cf7209c 100644 --- a/numpy/polynomial/tests/test_classes.py +++ b/numpy/polynomial/tests/test_classes.py @@ -173,11 +173,18 @@ def check_fit(Poly): assert_almost_equal(p(x), y) assert_almost_equal(p.domain, d) assert_almost_equal(p.window, w) + p = Poly.fit(x, y, [0, 1, 2, 3], domain=d, window=w) + assert_almost_equal(p(x), y) + assert_almost_equal(p.domain, d) + assert_almost_equal(p.window, w) # check with class domain default p = Poly.fit(x, y, 3, []) assert_equal(p.domain, Poly.domain) assert_equal(p.window, Poly.window) + p = Poly.fit(x, y, [0, 1, 2, 3], []) + assert_equal(p.domain, Poly.domain) + assert_equal(p.window, Poly.window) # check that fit accepts weights. w = np.zeros_like(x) @@ -185,7 +192,9 @@ def check_fit(Poly): w[::2] = 1 p1 = Poly.fit(x[::2], z[::2], 3) p2 = Poly.fit(x, z, 3, w=w) + p3 = Poly.fit(x, z, [0, 1, 2, 3], w=w) assert_almost_equal(p1(x), p2(x)) + assert_almost_equal(p2(x), p3(x)) def check_equal(Poly): -- cgit v1.2.1