Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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