1""" 2Utility classes and functions for the polynomial modules. 3 4This module provides: error and warning objects; a polynomial base class; 5and some routines used in both the `polynomial` and `chebyshev` modules. 6 7Error objects 8------------- 9 10.. autosummary:: 11 :toctree: generated/ 12 13 PolyError base class for this sub-package's errors. 14 PolyDomainError raised when domains are mismatched. 15 16Warning objects 17--------------- 18 19.. autosummary:: 20 :toctree: generated/ 21 22 RankWarning raised in least-squares fit for rank-deficient matrix. 23 24Base class 25---------- 26 27.. autosummary:: 28 :toctree: generated/ 29 30 PolyBase Obsolete base class for the polynomial classes. Do not use. 31 32Functions 33--------- 34 35.. autosummary:: 36 :toctree: generated/ 37 38 as_series convert list of array_likes into 1-D arrays of common type. 39 trimseq remove trailing zeros. 40 trimcoef remove small trailing coefficients. 41 getdomain return the domain appropriate for a given set of abscissae. 42 mapdomain maps points between domains. 43 mapparms parameters of the linear map between domains. 44 45""" 46import operator 47import functools 48import warnings 49 50import numpy as np 51 52__all__ = [ 53 'RankWarning', 'PolyError', 'PolyDomainError', 'as_series', 'trimseq', 54 'trimcoef', 'getdomain', 'mapdomain', 'mapparms', 'PolyBase'] 55 56# 57# Warnings and Exceptions 58# 59 60class RankWarning(UserWarning): 61 """Issued by chebfit when the design matrix is rank deficient.""" 62 pass 63 64class PolyError(Exception): 65 """Base class for errors in this module.""" 66 pass 67 68class PolyDomainError(PolyError): 69 """Issued by the generic Poly class when two domains don't match. 70 71 This is raised when an binary operation is passed Poly objects with 72 different domains. 73 74 """ 75 pass 76 77# 78# Base class for all polynomial types 79# 80 81class PolyBase: 82 """ 83 Base class for all polynomial types. 84 85 Deprecated in numpy 1.9.0, use the abstract 86 ABCPolyBase class instead. Note that the latter 87 requires a number of virtual functions to be 88 implemented. 89 90 """ 91 pass 92 93# 94# Helper functions to convert inputs to 1-D arrays 95# 96def trimseq(seq): 97 """Remove small Poly series coefficients. 98 99 Parameters 100 ---------- 101 seq : sequence 102 Sequence of Poly series coefficients. This routine fails for 103 empty sequences. 104 105 Returns 106 ------- 107 series : sequence 108 Subsequence with trailing zeros removed. If the resulting sequence 109 would be empty, return the first element. The returned sequence may 110 or may not be a view. 111 112 Notes 113 ----- 114 Do not lose the type info if the sequence contains unknown objects. 115 116 """ 117 if len(seq) == 0: 118 return seq 119 else: 120 for i in range(len(seq) - 1, -1, -1): 121 if seq[i] != 0: 122 break 123 return seq[:i+1] 124 125 126def as_series(alist, trim=True): 127 """ 128 Return argument as a list of 1-d arrays. 129 130 The returned list contains array(s) of dtype double, complex double, or 131 object. A 1-d argument of shape ``(N,)`` is parsed into ``N`` arrays of 132 size one; a 2-d argument of shape ``(M,N)`` is parsed into ``M`` arrays 133 of size ``N`` (i.e., is "parsed by row"); and a higher dimensional array 134 raises a Value Error if it is not first reshaped into either a 1-d or 2-d 135 array. 136 137 Parameters 138 ---------- 139 alist : array_like 140 A 1- or 2-d array_like 141 trim : boolean, optional 142 When True, trailing zeros are removed from the inputs. 143 When False, the inputs are passed through intact. 144 145 Returns 146 ------- 147 [a1, a2,...] : list of 1-D arrays 148 A copy of the input data as a list of 1-d arrays. 149 150 Raises 151 ------ 152 ValueError 153 Raised when `as_series` cannot convert its input to 1-d arrays, or at 154 least one of the resulting arrays is empty. 155 156 Examples 157 -------- 158 >>> from numpy.polynomial import polyutils as pu 159 >>> a = np.arange(4) 160 >>> pu.as_series(a) 161 [array([0.]), array([1.]), array([2.]), array([3.])] 162 >>> b = np.arange(6).reshape((2,3)) 163 >>> pu.as_series(b) 164 [array([0., 1., 2.]), array([3., 4., 5.])] 165 166 >>> pu.as_series((1, np.arange(3), np.arange(2, dtype=np.float16))) 167 [array([1.]), array([0., 1., 2.]), array([0., 1.])] 168 169 >>> pu.as_series([2, [1.1, 0.]]) 170 [array([2.]), array([1.1])] 171 172 >>> pu.as_series([2, [1.1, 0.]], trim=False) 173 [array([2.]), array([1.1, 0. ])] 174 175 """ 176 arrays = [np.array(a, ndmin=1, copy=False) for a in alist] 177 if min([a.size for a in arrays]) == 0: 178 raise ValueError("Coefficient array is empty") 179 if any(a.ndim != 1 for a in arrays): 180 raise ValueError("Coefficient array is not 1-d") 181 if trim: 182 arrays = [trimseq(a) for a in arrays] 183 184 if any(a.dtype == np.dtype(object) for a in arrays): 185 ret = [] 186 for a in arrays: 187 if a.dtype != np.dtype(object): 188 tmp = np.empty(len(a), dtype=np.dtype(object)) 189 tmp[:] = a[:] 190 ret.append(tmp) 191 else: 192 ret.append(a.copy()) 193 else: 194 try: 195 dtype = np.common_type(*arrays) 196 except Exception as e: 197 raise ValueError("Coefficient arrays have no common type") from e 198 ret = [np.array(a, copy=True, dtype=dtype) for a in arrays] 199 return ret 200 201 202def trimcoef(c, tol=0): 203 """ 204 Remove "small" "trailing" coefficients from a polynomial. 205 206 "Small" means "small in absolute value" and is controlled by the 207 parameter `tol`; "trailing" means highest order coefficient(s), e.g., in 208 ``[0, 1, 1, 0, 0]`` (which represents ``0 + x + x**2 + 0*x**3 + 0*x**4``) 209 both the 3-rd and 4-th order coefficients would be "trimmed." 210 211 Parameters 212 ---------- 213 c : array_like 214 1-d array of coefficients, ordered from lowest order to highest. 215 tol : number, optional 216 Trailing (i.e., highest order) elements with absolute value less 217 than or equal to `tol` (default value is zero) are removed. 218 219 Returns 220 ------- 221 trimmed : ndarray 222 1-d array with trailing zeros removed. If the resulting series 223 would be empty, a series containing a single zero is returned. 224 225 Raises 226 ------ 227 ValueError 228 If `tol` < 0 229 230 See Also 231 -------- 232 trimseq 233 234 Examples 235 -------- 236 >>> from numpy.polynomial import polyutils as pu 237 >>> pu.trimcoef((0,0,3,0,5,0,0)) 238 array([0., 0., 3., 0., 5.]) 239 >>> pu.trimcoef((0,0,1e-3,0,1e-5,0,0),1e-3) # item == tol is trimmed 240 array([0.]) 241 >>> i = complex(0,1) # works for complex 242 >>> pu.trimcoef((3e-4,1e-3*(1-i),5e-4,2e-5*(1+i)), 1e-3) 243 array([0.0003+0.j , 0.001 -0.001j]) 244 245 """ 246 if tol < 0: 247 raise ValueError("tol must be non-negative") 248 249 [c] = as_series([c]) 250 [ind] = np.nonzero(np.abs(c) > tol) 251 if len(ind) == 0: 252 return c[:1]*0 253 else: 254 return c[:ind[-1] + 1].copy() 255 256def getdomain(x): 257 """ 258 Return a domain suitable for given abscissae. 259 260 Find a domain suitable for a polynomial or Chebyshev series 261 defined at the values supplied. 262 263 Parameters 264 ---------- 265 x : array_like 266 1-d array of abscissae whose domain will be determined. 267 268 Returns 269 ------- 270 domain : ndarray 271 1-d array containing two values. If the inputs are complex, then 272 the two returned points are the lower left and upper right corners 273 of the smallest rectangle (aligned with the axes) in the complex 274 plane containing the points `x`. If the inputs are real, then the 275 two points are the ends of the smallest interval containing the 276 points `x`. 277 278 See Also 279 -------- 280 mapparms, mapdomain 281 282 Examples 283 -------- 284 >>> from numpy.polynomial import polyutils as pu 285 >>> points = np.arange(4)**2 - 5; points 286 array([-5, -4, -1, 4]) 287 >>> pu.getdomain(points) 288 array([-5., 4.]) 289 >>> c = np.exp(complex(0,1)*np.pi*np.arange(12)/6) # unit circle 290 >>> pu.getdomain(c) 291 array([-1.-1.j, 1.+1.j]) 292 293 """ 294 [x] = as_series([x], trim=False) 295 if x.dtype.char in np.typecodes['Complex']: 296 rmin, rmax = x.real.min(), x.real.max() 297 imin, imax = x.imag.min(), x.imag.max() 298 return np.array((complex(rmin, imin), complex(rmax, imax))) 299 else: 300 return np.array((x.min(), x.max())) 301 302def mapparms(old, new): 303 """ 304 Linear map parameters between domains. 305 306 Return the parameters of the linear map ``offset + scale*x`` that maps 307 `old` to `new` such that ``old[i] -> new[i]``, ``i = 0, 1``. 308 309 Parameters 310 ---------- 311 old, new : array_like 312 Domains. Each domain must (successfully) convert to a 1-d array 313 containing precisely two values. 314 315 Returns 316 ------- 317 offset, scale : scalars 318 The map ``L(x) = offset + scale*x`` maps the first domain to the 319 second. 320 321 See Also 322 -------- 323 getdomain, mapdomain 324 325 Notes 326 ----- 327 Also works for complex numbers, and thus can be used to calculate the 328 parameters required to map any line in the complex plane to any other 329 line therein. 330 331 Examples 332 -------- 333 >>> from numpy.polynomial import polyutils as pu 334 >>> pu.mapparms((-1,1),(-1,1)) 335 (0.0, 1.0) 336 >>> pu.mapparms((1,-1),(-1,1)) 337 (-0.0, -1.0) 338 >>> i = complex(0,1) 339 >>> pu.mapparms((-i,-1),(1,i)) 340 ((1+1j), (1-0j)) 341 342 """ 343 oldlen = old[1] - old[0] 344 newlen = new[1] - new[0] 345 off = (old[1]*new[0] - old[0]*new[1])/oldlen 346 scl = newlen/oldlen 347 return off, scl 348 349def mapdomain(x, old, new): 350 """ 351 Apply linear map to input points. 352 353 The linear map ``offset + scale*x`` that maps the domain `old` to 354 the domain `new` is applied to the points `x`. 355 356 Parameters 357 ---------- 358 x : array_like 359 Points to be mapped. If `x` is a subtype of ndarray the subtype 360 will be preserved. 361 old, new : array_like 362 The two domains that determine the map. Each must (successfully) 363 convert to 1-d arrays containing precisely two values. 364 365 Returns 366 ------- 367 x_out : ndarray 368 Array of points of the same shape as `x`, after application of the 369 linear map between the two domains. 370 371 See Also 372 -------- 373 getdomain, mapparms 374 375 Notes 376 ----- 377 Effectively, this implements: 378 379 .. math :: 380 x\\_out = new[0] + m(x - old[0]) 381 382 where 383 384 .. math :: 385 m = \\frac{new[1]-new[0]}{old[1]-old[0]} 386 387 Examples 388 -------- 389 >>> from numpy.polynomial import polyutils as pu 390 >>> old_domain = (-1,1) 391 >>> new_domain = (0,2*np.pi) 392 >>> x = np.linspace(-1,1,6); x 393 array([-1. , -0.6, -0.2, 0.2, 0.6, 1. ]) 394 >>> x_out = pu.mapdomain(x, old_domain, new_domain); x_out 395 array([ 0. , 1.25663706, 2.51327412, 3.76991118, 5.02654825, # may vary 396 6.28318531]) 397 >>> x - pu.mapdomain(x_out, new_domain, old_domain) 398 array([0., 0., 0., 0., 0., 0.]) 399 400 Also works for complex numbers (and thus can be used to map any line in 401 the complex plane to any other line therein). 402 403 >>> i = complex(0,1) 404 >>> old = (-1 - i, 1 + i) 405 >>> new = (-1 + i, 1 - i) 406 >>> z = np.linspace(old[0], old[1], 6); z 407 array([-1. -1.j , -0.6-0.6j, -0.2-0.2j, 0.2+0.2j, 0.6+0.6j, 1. +1.j ]) 408 >>> new_z = pu.mapdomain(z, old, new); new_z 409 array([-1.0+1.j , -0.6+0.6j, -0.2+0.2j, 0.2-0.2j, 0.6-0.6j, 1.0-1.j ]) # may vary 410 411 """ 412 x = np.asanyarray(x) 413 off, scl = mapparms(old, new) 414 return off + scl*x 415 416 417def _nth_slice(i, ndim): 418 sl = [np.newaxis] * ndim 419 sl[i] = slice(None) 420 return tuple(sl) 421 422 423def _vander_nd(vander_fs, points, degrees): 424 r""" 425 A generalization of the Vandermonde matrix for N dimensions 426 427 The result is built by combining the results of 1d Vandermonde matrices, 428 429 .. math:: 430 W[i_0, \ldots, i_M, j_0, \ldots, j_N] = \prod_{k=0}^N{V_k(x_k)[i_0, \ldots, i_M, j_k]} 431 432 where 433 434 .. math:: 435 N &= \texttt{len(points)} = \texttt{len(degrees)} = \texttt{len(vander\_fs)} \\ 436 M &= \texttt{points[k].ndim} \\ 437 V_k &= \texttt{vander\_fs[k]} \\ 438 x_k &= \texttt{points[k]} \\ 439 0 \le j_k &\le \texttt{degrees[k]} 440 441 Expanding the one-dimensional :math:`V_k` functions gives: 442 443 .. math:: 444 W[i_0, \ldots, i_M, j_0, \ldots, j_N] = \prod_{k=0}^N{B_{k, j_k}(x_k[i_0, \ldots, i_M])} 445 446 where :math:`B_{k,m}` is the m'th basis of the polynomial construction used along 447 dimension :math:`k`. For a regular polynomial, :math:`B_{k, m}(x) = P_m(x) = x^m`. 448 449 Parameters 450 ---------- 451 vander_fs : Sequence[function(array_like, int) -> ndarray] 452 The 1d vander function to use for each axis, such as ``polyvander`` 453 points : Sequence[array_like] 454 Arrays of point coordinates, all of the same shape. The dtypes 455 will be converted to either float64 or complex128 depending on 456 whether any of the elements are complex. Scalars are converted to 457 1-D arrays. 458 This must be the same length as `vander_fs`. 459 degrees : Sequence[int] 460 The maximum degree (inclusive) to use for each axis. 461 This must be the same length as `vander_fs`. 462 463 Returns 464 ------- 465 vander_nd : ndarray 466 An array of shape ``points[0].shape + tuple(d + 1 for d in degrees)``. 467 """ 468 n_dims = len(vander_fs) 469 if n_dims != len(points): 470 raise ValueError( 471 f"Expected {n_dims} dimensions of sample points, got {len(points)}") 472 if n_dims != len(degrees): 473 raise ValueError( 474 f"Expected {n_dims} dimensions of degrees, got {len(degrees)}") 475 if n_dims == 0: 476 raise ValueError("Unable to guess a dtype or shape when no points are given") 477 478 # convert to the same shape and type 479 points = tuple(np.array(tuple(points), copy=False) + 0.0) 480 481 # produce the vandermonde matrix for each dimension, placing the last 482 # axis of each in an independent trailing axis of the output 483 vander_arrays = ( 484 vander_fs[i](points[i], degrees[i])[(...,) + _nth_slice(i, n_dims)] 485 for i in range(n_dims) 486 ) 487 488 # we checked this wasn't empty already, so no `initial` needed 489 return functools.reduce(operator.mul, vander_arrays) 490 491 492def _vander_nd_flat(vander_fs, points, degrees): 493 """ 494 Like `_vander_nd`, but flattens the last ``len(degrees)`` axes into a single axis 495 496 Used to implement the public ``<type>vander<n>d`` functions. 497 """ 498 v = _vander_nd(vander_fs, points, degrees) 499 return v.reshape(v.shape[:-len(degrees)] + (-1,)) 500 501 502def _fromroots(line_f, mul_f, roots): 503 """ 504 Helper function used to implement the ``<type>fromroots`` functions. 505 506 Parameters 507 ---------- 508 line_f : function(float, float) -> ndarray 509 The ``<type>line`` function, such as ``polyline`` 510 mul_f : function(array_like, array_like) -> ndarray 511 The ``<type>mul`` function, such as ``polymul`` 512 roots : 513 See the ``<type>fromroots`` functions for more detail 514 """ 515 if len(roots) == 0: 516 return np.ones(1) 517 else: 518 [roots] = as_series([roots], trim=False) 519 roots.sort() 520 p = [line_f(-r, 1) for r in roots] 521 n = len(p) 522 while n > 1: 523 m, r = divmod(n, 2) 524 tmp = [mul_f(p[i], p[i+m]) for i in range(m)] 525 if r: 526 tmp[0] = mul_f(tmp[0], p[-1]) 527 p = tmp 528 n = m 529 return p[0] 530 531 532def _valnd(val_f, c, *args): 533 """ 534 Helper function used to implement the ``<type>val<n>d`` functions. 535 536 Parameters 537 ---------- 538 val_f : function(array_like, array_like, tensor: bool) -> array_like 539 The ``<type>val`` function, such as ``polyval`` 540 c, args : 541 See the ``<type>val<n>d`` functions for more detail 542 """ 543 args = [np.asanyarray(a) for a in args] 544 shape0 = args[0].shape 545 if not all((a.shape == shape0 for a in args[1:])): 546 if len(args) == 3: 547 raise ValueError('x, y, z are incompatible') 548 elif len(args) == 2: 549 raise ValueError('x, y are incompatible') 550 else: 551 raise ValueError('ordinates are incompatible') 552 it = iter(args) 553 x0 = next(it) 554 555 # use tensor on only the first 556 c = val_f(x0, c) 557 for xi in it: 558 c = val_f(xi, c, tensor=False) 559 return c 560 561 562def _gridnd(val_f, c, *args): 563 """ 564 Helper function used to implement the ``<type>grid<n>d`` functions. 565 566 Parameters 567 ---------- 568 val_f : function(array_like, array_like, tensor: bool) -> array_like 569 The ``<type>val`` function, such as ``polyval`` 570 c, args : 571 See the ``<type>grid<n>d`` functions for more detail 572 """ 573 for xi in args: 574 c = val_f(xi, c) 575 return c 576 577 578def _div(mul_f, c1, c2): 579 """ 580 Helper function used to implement the ``<type>div`` functions. 581 582 Implementation uses repeated subtraction of c2 multiplied by the nth basis. 583 For some polynomial types, a more efficient approach may be possible. 584 585 Parameters 586 ---------- 587 mul_f : function(array_like, array_like) -> array_like 588 The ``<type>mul`` function, such as ``polymul`` 589 c1, c2 : 590 See the ``<type>div`` functions for more detail 591 """ 592 # c1, c2 are trimmed copies 593 [c1, c2] = as_series([c1, c2]) 594 if c2[-1] == 0: 595 raise ZeroDivisionError() 596 597 lc1 = len(c1) 598 lc2 = len(c2) 599 if lc1 < lc2: 600 return c1[:1]*0, c1 601 elif lc2 == 1: 602 return c1/c2[-1], c1[:1]*0 603 else: 604 quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype) 605 rem = c1 606 for i in range(lc1 - lc2, - 1, -1): 607 p = mul_f([0]*i + [1], c2) 608 q = rem[-1]/p[-1] 609 rem = rem[:-1] - q*p[:-1] 610 quo[i] = q 611 return quo, trimseq(rem) 612 613 614def _add(c1, c2): 615 """ Helper function used to implement the ``<type>add`` functions. """ 616 # c1, c2 are trimmed copies 617 [c1, c2] = as_series([c1, c2]) 618 if len(c1) > len(c2): 619 c1[:c2.size] += c2 620 ret = c1 621 else: 622 c2[:c1.size] += c1 623 ret = c2 624 return trimseq(ret) 625 626 627def _sub(c1, c2): 628 """ Helper function used to implement the ``<type>sub`` functions. """ 629 # c1, c2 are trimmed copies 630 [c1, c2] = as_series([c1, c2]) 631 if len(c1) > len(c2): 632 c1[:c2.size] -= c2 633 ret = c1 634 else: 635 c2 = -c2 636 c2[:c1.size] += c1 637 ret = c2 638 return trimseq(ret) 639 640 641def _fit(vander_f, x, y, deg, rcond=None, full=False, w=None): 642 """ 643 Helper function used to implement the ``<type>fit`` functions. 644 645 Parameters 646 ---------- 647 vander_f : function(array_like, int) -> ndarray 648 The 1d vander function, such as ``polyvander`` 649 c1, c2 : 650 See the ``<type>fit`` functions for more detail 651 """ 652 x = np.asarray(x) + 0.0 653 y = np.asarray(y) + 0.0 654 deg = np.asarray(deg) 655 656 # check arguments. 657 if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0: 658 raise TypeError("deg must be an int or non-empty 1-D array of int") 659 if deg.min() < 0: 660 raise ValueError("expected deg >= 0") 661 if x.ndim != 1: 662 raise TypeError("expected 1D vector for x") 663 if x.size == 0: 664 raise TypeError("expected non-empty vector for x") 665 if y.ndim < 1 or y.ndim > 2: 666 raise TypeError("expected 1D or 2D array for y") 667 if len(x) != len(y): 668 raise TypeError("expected x and y to have same length") 669 670 if deg.ndim == 0: 671 lmax = deg 672 order = lmax + 1 673 van = vander_f(x, lmax) 674 else: 675 deg = np.sort(deg) 676 lmax = deg[-1] 677 order = len(deg) 678 van = vander_f(x, lmax)[:, deg] 679 680 # set up the least squares matrices in transposed form 681 lhs = van.T 682 rhs = y.T 683 if w is not None: 684 w = np.asarray(w) + 0.0 685 if w.ndim != 1: 686 raise TypeError("expected 1D vector for w") 687 if len(x) != len(w): 688 raise TypeError("expected x and w to have same length") 689 # apply weights. Don't use inplace operations as they 690 # can cause problems with NA. 691 lhs = lhs * w 692 rhs = rhs * w 693 694 # set rcond 695 if rcond is None: 696 rcond = len(x)*np.finfo(x.dtype).eps 697 698 # Determine the norms of the design matrix columns. 699 if issubclass(lhs.dtype.type, np.complexfloating): 700 scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1)) 701 else: 702 scl = np.sqrt(np.square(lhs).sum(1)) 703 scl[scl == 0] = 1 704 705 # Solve the least squares problem. 706 c, resids, rank, s = np.linalg.lstsq(lhs.T/scl, rhs.T, rcond) 707 c = (c.T/scl).T 708 709 # Expand c to include non-fitted coefficients which are set to zero 710 if deg.ndim > 0: 711 if c.ndim == 2: 712 cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype) 713 else: 714 cc = np.zeros(lmax+1, dtype=c.dtype) 715 cc[deg] = c 716 c = cc 717 718 # warn on rank reduction 719 if rank != order and not full: 720 msg = "The fit may be poorly conditioned" 721 warnings.warn(msg, RankWarning, stacklevel=2) 722 723 if full: 724 return c, [resids, rank, s, rcond] 725 else: 726 return c 727 728 729def _pow(mul_f, c, pow, maxpower): 730 """ 731 Helper function used to implement the ``<type>pow`` functions. 732 733 Parameters 734 ---------- 735 vander_f : function(array_like, int) -> ndarray 736 The 1d vander function, such as ``polyvander`` 737 pow, maxpower : 738 See the ``<type>pow`` functions for more detail 739 mul_f : function(array_like, array_like) -> ndarray 740 The ``<type>mul`` function, such as ``polymul`` 741 """ 742 # c is a trimmed copy 743 [c] = as_series([c]) 744 power = int(pow) 745 if power != pow or power < 0: 746 raise ValueError("Power must be a non-negative integer.") 747 elif maxpower is not None and power > maxpower: 748 raise ValueError("Power is too large") 749 elif power == 0: 750 return np.array([1], dtype=c.dtype) 751 elif power == 1: 752 return c 753 else: 754 # This can be made more efficient by using powers of two 755 # in the usual way. 756 prd = c 757 for i in range(2, power + 1): 758 prd = mul_f(prd, c) 759 return prd 760 761 762def _deprecate_as_int(x, desc): 763 """ 764 Like `operator.index`, but emits a deprecation warning when passed a float 765 766 Parameters 767 ---------- 768 x : int-like, or float with integral value 769 Value to interpret as an integer 770 desc : str 771 description to include in any error message 772 773 Raises 774 ------ 775 TypeError : if x is a non-integral float or non-numeric 776 DeprecationWarning : if x is an integral float 777 """ 778 try: 779 return operator.index(x) 780 except TypeError as e: 781 # Numpy 1.17.0, 2019-03-11 782 try: 783 ix = int(x) 784 except TypeError: 785 pass 786 else: 787 if ix == x: 788 warnings.warn( 789 f"In future, this will raise TypeError, as {desc} will " 790 "need to be an integer not just an integral float.", 791 DeprecationWarning, 792 stacklevel=3 793 ) 794 return ix 795 796 raise TypeError(f"{desc} must be an integer") from e 797