Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_minpack_py.py: 9%

1import warnings 

2from . import _minpack 

3 

4import numpy as np 

5from numpy import (atleast_1d, triu, shape, transpose, zeros, prod, greater, 

6 asarray, inf, 

7 finfo, inexact, issubdtype, dtype) 

8from scipy import linalg 

9from scipy.linalg import svd, cholesky, solve_triangular, LinAlgError 

10from scipy._lib._util import _asarray_validated, _lazywhere, _contains_nan 

11from scipy._lib._util import getfullargspec_no_self as _getfullargspec 

12from ._optimize import OptimizeResult, _check_unknown_options, OptimizeWarning 

13from ._lsq import least_squares 

14# from ._lsq.common import make_strictly_feasible 

15from ._lsq.least_squares import prepare_bounds 

16from scipy.optimize._minimize import Bounds 

17 

18error = _minpack.error 

19 

20__all__ = ['fsolve', 'leastsq', 'fixed_point', 'curve_fit'] 

21 

22 

23def _check_func(checker, argname, thefunc, x0, args, numinputs, 

24 output_shape=None): 

25 res = atleast_1d(thefunc(*((x0[:numinputs],) + args))) 

26 if (output_shape is not None) and (shape(res) != output_shape): 

27 if (output_shape[0] != 1): 

28 if len(output_shape) > 1: 

29 if output_shape[1] == 1: 

30 return shape(res) 

31 msg = "{}: there is a mismatch between the input and output " \ 

32 "shape of the '{}' argument".format(checker, argname) 

33 func_name = getattr(thefunc, '__name__', None) 

34 if func_name: 

35 msg += " '%s'." % func_name 

36 else: 

37 msg += "." 

38 msg += f'Shape should be {output_shape} but it is {shape(res)}.' 

39 raise TypeError(msg) 

40 if issubdtype(res.dtype, inexact): 

41 dt = res.dtype 

42 else: 

43 dt = dtype(float) 

44 return shape(res), dt 

45 

46 

47def fsolve(func, x0, args=(), fprime=None, full_output=0, 

48 col_deriv=0, xtol=1.49012e-8, maxfev=0, band=None, 

49 epsfcn=None, factor=100, diag=None): 

50 """ 

51 Find the roots of a function. 

52 

53 Return the roots of the (non-linear) equations defined by 

54 ``func(x) = 0`` given a starting estimate. 

55 

56 Parameters 

57 ---------- 

58 func : callable ``f(x, *args)`` 

59 A function that takes at least one (possibly vector) argument, 

60 and returns a value of the same length. 

61 x0 : ndarray 

62 The starting estimate for the roots of ``func(x) = 0``. 

63 args : tuple, optional 

64 Any extra arguments to `func`. 

65 fprime : callable ``f(x, *args)``, optional 

66 A function to compute the Jacobian of `func` with derivatives 

67 across the rows. By default, the Jacobian will be estimated. 

68 full_output : bool, optional 

69 If True, return optional outputs. 

70 col_deriv : bool, optional 

71 Specify whether the Jacobian function computes derivatives down 

72 the columns (faster, because there is no transpose operation). 

73 xtol : float, optional 

74 The calculation will terminate if the relative error between two 

75 consecutive iterates is at most `xtol`. 

76 maxfev : int, optional 

77 The maximum number of calls to the function. If zero, then 

78 ``100*(N+1)`` is the maximum where N is the number of elements 

79 in `x0`. 

80 band : tuple, optional 

81 If set to a two-sequence containing the number of sub- and 

82 super-diagonals within the band of the Jacobi matrix, the 

83 Jacobi matrix is considered banded (only for ``fprime=None``). 

84 epsfcn : float, optional 

85 A suitable step length for the forward-difference 

86 approximation of the Jacobian (for ``fprime=None``). If 

87 `epsfcn` is less than the machine precision, it is assumed 

88 that the relative errors in the functions are of the order of 

89 the machine precision. 

90 factor : float, optional 

91 A parameter determining the initial step bound 

92 (``factor * || diag * x||``). Should be in the interval 

93 ``(0.1, 100)``. 

94 diag : sequence, optional 

95 N positive entries that serve as a scale factors for the 

96 variables. 

97 

98 Returns 

99 ------- 

100 x : ndarray 

101 The solution (or the result of the last iteration for 

102 an unsuccessful call). 

103 infodict : dict 

104 A dictionary of optional outputs with the keys: 

105 

106 ``nfev`` 

107 number of function calls 

108 ``njev`` 

109 number of Jacobian calls 

110 ``fvec`` 

111 function evaluated at the output 

112 ``fjac`` 

113 the orthogonal matrix, q, produced by the QR 

114 factorization of the final approximate Jacobian 

115 matrix, stored column wise 

116 ``r`` 

117 upper triangular matrix produced by QR factorization 

118 of the same matrix 

119 ``qtf`` 

120 the vector ``(transpose(q) * fvec)`` 

121 

122 ier : int 

123 An integer flag. Set to 1 if a solution was found, otherwise refer 

124 to `mesg` for more information. 

125 mesg : str 

126 If no solution is found, `mesg` details the cause of failure. 

127 

128 See Also 

129 -------- 

130 root : Interface to root finding algorithms for multivariate 

131 functions. See the ``method='hybr'`` in particular. 

132 

133 Notes 

134 ----- 

135 ``fsolve`` is a wrapper around MINPACK's hybrd and hybrj algorithms. 

136 

137 Examples 

138 -------- 

139 Find a solution to the system of equations: 

140 ``x0*cos(x1) = 4, x1*x0 - x1 = 5``. 

141 

142 >>> import numpy as np 

143 >>> from scipy.optimize import fsolve 

144 >>> def func(x): 

145 ... return [x[0] * np.cos(x[1]) - 4, 

146 ... x[1] * x[0] - x[1] - 5] 

147 >>> root = fsolve(func, [1, 1]) 

148 >>> root 

149 array([6.50409711, 0.90841421]) 

150 >>> np.isclose(func(root), [0.0, 0.0]) # func(root) should be almost 0.0. 

151 array([ True, True]) 

152 

153 """ 

154 options = {'col_deriv': col_deriv, 

155 'xtol': xtol, 

156 'maxfev': maxfev, 

157 'band': band, 

158 'eps': epsfcn, 

159 'factor': factor, 

160 'diag': diag} 

161 

162 res = _root_hybr(func, x0, args, jac=fprime, **options) 

163 if full_output: 

164 x = res['x'] 

165 info = {k: res.get(k) 

166 for k in ('nfev', 'njev', 'fjac', 'r', 'qtf') if k in res} 

167 info['fvec'] = res['fun'] 

168 return x, info, res['status'], res['message'] 

169 else: 

170 status = res['status'] 

171 msg = res['message'] 

172 if status == 0: 

173 raise TypeError(msg) 

174 elif status == 1: 

175 pass 

176 elif status in [2, 3, 4, 5]: 

177 warnings.warn(msg, RuntimeWarning) 

178 else: 

179 raise TypeError(msg) 

180 return res['x'] 

181 

182 

183def _root_hybr(func, x0, args=(), jac=None, 

184 col_deriv=0, xtol=1.49012e-08, maxfev=0, band=None, eps=None, 

185 factor=100, diag=None, **unknown_options): 

186 """ 

187 Find the roots of a multivariate function using MINPACK's hybrd and 

188 hybrj routines (modified Powell method). 

189 

190 Options 

191 ------- 

192 col_deriv : bool 

193 Specify whether the Jacobian function computes derivatives down 

194 the columns (faster, because there is no transpose operation). 

195 xtol : float 

196 The calculation will terminate if the relative error between two 

197 consecutive iterates is at most `xtol`. 

198 maxfev : int 

199 The maximum number of calls to the function. If zero, then 

200 ``100*(N+1)`` is the maximum where N is the number of elements 

201 in `x0`. 

202 band : tuple 

203 If set to a two-sequence containing the number of sub- and 

204 super-diagonals within the band of the Jacobi matrix, the 

205 Jacobi matrix is considered banded (only for ``fprime=None``). 

206 eps : float 

207 A suitable step length for the forward-difference 

208 approximation of the Jacobian (for ``fprime=None``). If 

209 `eps` is less than the machine precision, it is assumed 

210 that the relative errors in the functions are of the order of 

211 the machine precision. 

212 factor : float 

213 A parameter determining the initial step bound 

214 (``factor * || diag * x||``). Should be in the interval 

215 ``(0.1, 100)``. 

216 diag : sequence 

217 N positive entries that serve as a scale factors for the 

218 variables. 

219 

220 """ 

221 _check_unknown_options(unknown_options) 

222 epsfcn = eps 

223 

224 x0 = asarray(x0).flatten() 

225 n = len(x0) 

226 if not isinstance(args, tuple): 

227 args = (args,) 

228 shape, dtype = _check_func('fsolve', 'func', func, x0, args, n, (n,)) 

229 if epsfcn is None: 

230 epsfcn = finfo(dtype).eps 

231 Dfun = jac 

232 if Dfun is None: 

233 if band is None: 

234 ml, mu = -10, -10 

235 else: 

236 ml, mu = band[:2] 

237 if maxfev == 0: 

238 maxfev = 200 * (n + 1) 

239 retval = _minpack._hybrd(func, x0, args, 1, xtol, maxfev, 

240 ml, mu, epsfcn, factor, diag) 

241 else: 

242 _check_func('fsolve', 'fprime', Dfun, x0, args, n, (n, n)) 

243 if (maxfev == 0): 

244 maxfev = 100 * (n + 1) 

245 retval = _minpack._hybrj(func, Dfun, x0, args, 1, 

246 col_deriv, xtol, maxfev, factor, diag) 

247 

248 x, status = retval[0], retval[-1] 

249 

250 errors = {0: "Improper input parameters were entered.", 

251 1: "The solution converged.", 

252 2: "The number of calls to function has " 

253 "reached maxfev = %d." % maxfev, 

254 3: "xtol=%f is too small, no further improvement " 

255 "in the approximate\n solution " 

256 "is possible." % xtol, 

257 4: "The iteration is not making good progress, as measured " 

258 "by the \n improvement from the last five " 

259 "Jacobian evaluations.", 

260 5: "The iteration is not making good progress, " 

261 "as measured by the \n improvement from the last " 

262 "ten iterations.", 

263 'unknown': "An error occurred."} 

264 

265 info = retval[1] 

266 info['fun'] = info.pop('fvec') 

267 sol = OptimizeResult(x=x, success=(status == 1), status=status) 

268 sol.update(info) 

269 try: 

270 sol['message'] = errors[status] 

271 except KeyError: 

272 sol['message'] = errors['unknown'] 

273 

274 return sol 

275 

276 

277LEASTSQ_SUCCESS = [1, 2, 3, 4] 

278LEASTSQ_FAILURE = [5, 6, 7, 8] 

279 

280 

281def leastsq(func, x0, args=(), Dfun=None, full_output=False, 

282 col_deriv=False, ftol=1.49012e-8, xtol=1.49012e-8, 

283 gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None): 

284 """ 

285 Minimize the sum of squares of a set of equations. 

286 

287 :: 

288 

289 x = arg min(sum(func(y)**2,axis=0)) 

290 y 

291 

292 Parameters 

293 ---------- 

294 func : callable 

295 Should take at least one (possibly length ``N`` vector) argument and 

296 returns ``M`` floating point numbers. It must not return NaNs or 

297 fitting might fail. ``M`` must be greater than or equal to ``N``. 

298 x0 : ndarray 

299 The starting estimate for the minimization. 

300 args : tuple, optional 

301 Any extra arguments to func are placed in this tuple. 

302 Dfun : callable, optional 

303 A function or method to compute the Jacobian of func with derivatives 

304 across the rows. If this is None, the Jacobian will be estimated. 

305 full_output : bool, optional 

306 If ``True``, return all optional outputs (not just `x` and `ier`). 

307 col_deriv : bool, optional 

308 If ``True``, specify that the Jacobian function computes derivatives 

309 down the columns (faster, because there is no transpose operation). 

310 ftol : float, optional 

311 Relative error desired in the sum of squares. 

312 xtol : float, optional 

313 Relative error desired in the approximate solution. 

314 gtol : float, optional 

315 Orthogonality desired between the function vector and the columns of 

316 the Jacobian. 

317 maxfev : int, optional 

318 The maximum number of calls to the function. If `Dfun` is provided, 

319 then the default `maxfev` is 100*(N+1) where N is the number of elements 

320 in x0, otherwise the default `maxfev` is 200*(N+1). 

321 epsfcn : float, optional 

322 A variable used in determining a suitable step length for the forward- 

323 difference approximation of the Jacobian (for Dfun=None). 

324 Normally the actual step length will be sqrt(epsfcn)*x 

325 If epsfcn is less than the machine precision, it is assumed that the 

326 relative errors are of the order of the machine precision. 

327 factor : float, optional 

328 A parameter determining the initial step bound 

329 (``factor * || diag * x||``). Should be in interval ``(0.1, 100)``. 

330 diag : sequence, optional 

331 N positive entries that serve as a scale factors for the variables. 

332 

333 Returns 

334 ------- 

335 x : ndarray 

336 The solution (or the result of the last iteration for an unsuccessful 

337 call). 

338 cov_x : ndarray 

339 The inverse of the Hessian. `fjac` and `ipvt` are used to construct an 

340 estimate of the Hessian. A value of None indicates a singular matrix, 

341 which means the curvature in parameters `x` is numerically flat. To 

342 obtain the covariance matrix of the parameters `x`, `cov_x` must be 

343 multiplied by the variance of the residuals -- see curve_fit. Only 

344 returned if `full_output` is ``True``. 

345 infodict : dict 

346 a dictionary of optional outputs with the keys: 

347 

348 ``nfev`` 

349 The number of function calls 

350 ``fvec`` 

351 The function evaluated at the output 

352 ``fjac`` 

353 A permutation of the R matrix of a QR 

354 factorization of the final approximate 

355 Jacobian matrix, stored column wise. 

356 Together with ipvt, the covariance of the 

357 estimate can be approximated. 

358 ``ipvt`` 

359 An integer array of length N which defines 

360 a permutation matrix, p, such that 

361 fjac*p = q*r, where r is upper triangular 

362 with diagonal elements of nonincreasing 

363 magnitude. Column j of p is column ipvt(j) 

364 of the identity matrix. 

365 ``qtf`` 

366 The vector (transpose(q) * fvec). 

367 

368 Only returned if `full_output` is ``True``. 

369 mesg : str 

370 A string message giving information about the cause of failure. 

371 Only returned if `full_output` is ``True``. 

372 ier : int 

373 An integer flag. If it is equal to 1, 2, 3 or 4, the solution was 

374 found. Otherwise, the solution was not found. In either case, the 

375 optional output variable 'mesg' gives more information. 

376 

377 See Also 

378 -------- 

379 least_squares : Newer interface to solve nonlinear least-squares problems 

380 with bounds on the variables. See ``method='lm'`` in particular. 

381 

382 Notes 

383 ----- 

384 "leastsq" is a wrapper around MINPACK's lmdif and lmder algorithms. 

385 

386 cov_x is a Jacobian approximation to the Hessian of the least squares 

387 objective function. 

388 This approximation assumes that the objective function is based on the 

389 difference between some observed target data (ydata) and a (non-linear) 

390 function of the parameters `f(xdata, params)` :: 

391 

392 func(params) = ydata - f(xdata, params) 

393 

394 so that the objective function is :: 

395 

396 min sum((ydata - f(xdata, params))**2, axis=0) 

397 params 

398 

399 The solution, `x`, is always a 1-D array, regardless of the shape of `x0`, 

400 or whether `x0` is a scalar. 

401 

402 Examples 

403 -------- 

404 >>> from scipy.optimize import leastsq 

405 >>> def func(x): 

406 ... return 2*(x-3)**2+1 

407 >>> leastsq(func, 0) 

408 (array([2.99999999]), 1) 

409 

410 """ 

411 x0 = asarray(x0).flatten() 

412 n = len(x0) 

413 if not isinstance(args, tuple): 

414 args = (args,) 

415 shape, dtype = _check_func('leastsq', 'func', func, x0, args, n) 

416 m = shape[0] 

417 

418 if n > m: 

419 raise TypeError(f"Improper input: func input vector length N={n} must" 

420 f" not exceed func output vector length M={m}") 

421 

422 if epsfcn is None: 

423 epsfcn = finfo(dtype).eps 

424 

425 if Dfun is None: 

426 if maxfev == 0: 

427 maxfev = 200*(n + 1) 

428 retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol, 

429 gtol, maxfev, epsfcn, factor, diag) 

430 else: 

431 if col_deriv: 

432 _check_func('leastsq', 'Dfun', Dfun, x0, args, n, (n, m)) 

433 else: 

434 _check_func('leastsq', 'Dfun', Dfun, x0, args, n, (m, n)) 

435 if maxfev == 0: 

436 maxfev = 100 * (n + 1) 

437 retval = _minpack._lmder(func, Dfun, x0, args, full_output, 

438 col_deriv, ftol, xtol, gtol, maxfev, 

439 factor, diag) 

440 

441 errors = {0: ["Improper input parameters.", TypeError], 

442 1: ["Both actual and predicted relative reductions " 

443 "in the sum of squares\n are at most %f" % ftol, None], 

444 2: ["The relative error between two consecutive " 

445 "iterates is at most %f" % xtol, None], 

446 3: ["Both actual and predicted relative reductions in " 

447 "the sum of squares\n are at most {:f} and the " 

448 "relative error between two consecutive " 

449 "iterates is at \n most {:f}".format(ftol, xtol), None], 

450 4: ["The cosine of the angle between func(x) and any " 

451 "column of the\n Jacobian is at most %f in " 

452 "absolute value" % gtol, None], 

453 5: ["Number of calls to function has reached " 

454 "maxfev = %d." % maxfev, ValueError], 

455 6: ["ftol=%f is too small, no further reduction " 

456 "in the sum of squares\n is possible." % ftol, 

457 ValueError], 

458 7: ["xtol=%f is too small, no further improvement in " 

459 "the approximate\n solution is possible." % xtol, 

460 ValueError], 

461 8: ["gtol=%f is too small, func(x) is orthogonal to the " 

462 "columns of\n the Jacobian to machine " 

463 "precision." % gtol, ValueError]} 

464 

465 # The FORTRAN return value (possible return values are >= 0 and <= 8) 

466 info = retval[-1] 

467 

468 if full_output: 

469 cov_x = None 

470 if info in LEASTSQ_SUCCESS: 

471 # This was 

472 # perm = take(eye(n), retval[1]['ipvt'] - 1, 0) 

473 # r = triu(transpose(retval[1]['fjac'])[:n, :]) 

474 # R = dot(r, perm) 

475 # cov_x = inv(dot(transpose(R), R)) 

476 # but the explicit dot product was not necessary and sometimes 

477 # the result was not symmetric positive definite. See gh-4555. 

478 perm = retval[1]['ipvt'] - 1 

479 n = len(perm) 

480 r = triu(transpose(retval[1]['fjac'])[:n, :]) 

481 inv_triu = linalg.get_lapack_funcs('trtri', (r,)) 

482 try: 

483 # inverse of permuted matrix is a permuation of matrix inverse 

484 invR, trtri_info = inv_triu(r) # default: upper, non-unit diag 

485 if trtri_info != 0: # explicit comparison for readability 

486 raise LinAlgError(f'trtri returned info {trtri_info}') 

487 invR[perm] = invR.copy() 

488 cov_x = invR @ invR.T 

489 except (LinAlgError, ValueError): 

490 pass 

491 return (retval[0], cov_x) + retval[1:-1] + (errors[info][0], info) 

492 else: 

493 if info in LEASTSQ_FAILURE: 

494 warnings.warn(errors[info][0], RuntimeWarning) 

495 elif info == 0: 

496 raise errors[info][1](errors[info][0]) 

497 return retval[0], info 

498 

499 

500def _lightweight_memoizer(f): 

501 # very shallow memoization - only remember the first set of parameters 

502 # and corresponding function value to address gh-13670 

503 def _memoized_func(params): 

504 if np.all(_memoized_func.last_params == params): 

505 return _memoized_func.last_val 

506 

507 val = f(params) 

508 

509 if _memoized_func.last_params is None: 

510 _memoized_func.last_params = np.copy(params) 

511 _memoized_func.last_val = val 

512 

513 return val 

514 

515 _memoized_func.last_params = None 

516 _memoized_func.last_val = None 

517 return _memoized_func 

518 

519 

520def _wrap_func(func, xdata, ydata, transform): 

521 if transform is None: 

522 def func_wrapped(params): 

523 return func(xdata, *params) - ydata 

524 elif transform.ndim == 1: 

525 def func_wrapped(params): 

526 return transform * (func(xdata, *params) - ydata) 

527 else: 

528 # Chisq = (y - yd)^T C^{-1} (y-yd) 

529 # transform = L such that C = L L^T 

530 # C^{-1} = L^{-T} L^{-1} 

531 # Chisq = (y - yd)^T L^{-T} L^{-1} (y-yd) 

532 # Define (y-yd)' = L^{-1} (y-yd) 

533 # by solving 

534 # L (y-yd)' = (y-yd) 

535 # and minimize (y-yd)'^T (y-yd)' 

536 def func_wrapped(params): 

537 return solve_triangular(transform, func(xdata, *params) - ydata, lower=True) 

538 return func_wrapped 

539 

540 

541def _wrap_jac(jac, xdata, transform): 

542 if transform is None: 

543 def jac_wrapped(params): 

544 return jac(xdata, *params) 

545 elif transform.ndim == 1: 

546 def jac_wrapped(params): 

547 return transform[:, np.newaxis] * np.asarray(jac(xdata, *params)) 

548 else: 

549 def jac_wrapped(params): 

550 return solve_triangular(transform, np.asarray(jac(xdata, *params)), lower=True) 

551 return jac_wrapped 

552 

553 

554def _initialize_feasible(lb, ub): 

555 p0 = np.ones_like(lb) 

556 lb_finite = np.isfinite(lb) 

557 ub_finite = np.isfinite(ub) 

558 

559 mask = lb_finite & ub_finite 

560 p0[mask] = 0.5 * (lb[mask] + ub[mask]) 

561 

562 mask = lb_finite & ~ub_finite 

563 p0[mask] = lb[mask] + 1 

564 

565 mask = ~lb_finite & ub_finite 

566 p0[mask] = ub[mask] - 1 

567 

568 return p0 

569 

570 

571def curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, 

572 check_finite=None, bounds=(-np.inf, np.inf), method=None, 

573 jac=None, *, full_output=False, nan_policy=None, 

574 **kwargs): 

575 """ 

576 Use non-linear least squares to fit a function, f, to data. 

577 

578 Assumes ``ydata = f(xdata, *params) + eps``. 

579 

580 Parameters 

581 ---------- 

582 f : callable 

583 The model function, f(x, ...). It must take the independent 

584 variable as the first argument and the parameters to fit as 

585 separate remaining arguments. 

586 xdata : array_like 

587 The independent variable where the data is measured. 

588 Should usually be an M-length sequence or an (k,M)-shaped array for 

589 functions with k predictors, and each element should be float 

590 convertible if it is an array like object. 

591 ydata : array_like 

592 The dependent data, a length M array - nominally ``f(xdata, ...)``. 

593 p0 : array_like, optional 

594 Initial guess for the parameters (length N). If None, then the 

595 initial values will all be 1 (if the number of parameters for the 

596 function can be determined using introspection, otherwise a 

597 ValueError is raised). 

598 sigma : None or M-length sequence or MxM array, optional 

599 Determines the uncertainty in `ydata`. If we define residuals as 

600 ``r = ydata - f(xdata, *popt)``, then the interpretation of `sigma` 

601 depends on its number of dimensions: 

602 

603 - A 1-D `sigma` should contain values of standard deviations of 

604 errors in `ydata`. In this case, the optimized function is 

605 ``chisq = sum((r / sigma) ** 2)``. 

606 

607 - A 2-D `sigma` should contain the covariance matrix of 

608 errors in `ydata`. In this case, the optimized function is 

609 ``chisq = r.T @ inv(sigma) @ r``. 

610 

611 .. versionadded:: 0.19 

612 

613 None (default) is equivalent of 1-D `sigma` filled with ones. 

614 absolute_sigma : bool, optional 

615 If True, `sigma` is used in an absolute sense and the estimated parameter 

616 covariance `pcov` reflects these absolute values. 

617 

618 If False (default), only the relative magnitudes of the `sigma` values matter. 

619 The returned parameter covariance matrix `pcov` is based on scaling 

620 `sigma` by a constant factor. This constant is set by demanding that the 

621 reduced `chisq` for the optimal parameters `popt` when using the 

622 *scaled* `sigma` equals unity. In other words, `sigma` is scaled to 

623 match the sample variance of the residuals after the fit. Default is False. 

624 Mathematically, 

625 ``pcov(absolute_sigma=False) = pcov(absolute_sigma=True) * chisq(popt)/(M-N)`` 

626 check_finite : bool, optional 

627 If True, check that the input arrays do not contain nans of infs, 

628 and raise a ValueError if they do. Setting this parameter to 

629 False may silently produce nonsensical results if the input arrays 

630 do contain nans. Default is True if `nan_policy` is not specified 

631 explicitly and False otherwise. 

632 bounds : 2-tuple of array_like or `Bounds`, optional 

633 Lower and upper bounds on parameters. Defaults to no bounds. 

634 There are two ways to specify the bounds: 

635 

636 - Instance of `Bounds` class. 

637 

638 - 2-tuple of array_like: Each element of the tuple must be either 

639 an array with the length equal to the number of parameters, or a 

640 scalar (in which case the bound is taken to be the same for all 

641 parameters). Use ``np.inf`` with an appropriate sign to disable 

642 bounds on all or some parameters. 

643 

644 method : {'lm', 'trf', 'dogbox'}, optional 

645 Method to use for optimization. See `least_squares` for more details. 

646 Default is 'lm' for unconstrained problems and 'trf' if `bounds` are 

647 provided. The method 'lm' won't work when the number of observations 

648 is less than the number of variables, use 'trf' or 'dogbox' in this 

649 case. 

650 

651 .. versionadded:: 0.17 

652 jac : callable, string or None, optional 

653 Function with signature ``jac(x, ...)`` which computes the Jacobian 

654 matrix of the model function with respect to parameters as a dense 

655 array_like structure. It will be scaled according to provided `sigma`. 

656 If None (default), the Jacobian will be estimated numerically. 

657 String keywords for 'trf' and 'dogbox' methods can be used to select 

658 a finite difference scheme, see `least_squares`. 

659 

660 .. versionadded:: 0.18 

661 full_output : boolean, optional 

662 If True, this function returns additioal information: `infodict`, 

663 `mesg`, and `ier`. 

664 

665 .. versionadded:: 1.9 

666 nan_policy : {'raise', 'omit', None}, optional 

667 Defines how to handle when input contains nan. 

668 The following options are available (default is None): 

669 

670 * 'raise': throws an error 

671 * 'omit': performs the calculations ignoring nan values 

672 * None: no special handling of NaNs is performed 

673 (except what is done by check_finite); the behavior when NaNs 

674 are present is implementation-dependent and may change. 

675 

676 Note that if this value is specified explicitly (not None), 

677 `check_finite` will be set as False. 

678 

679 .. versionadded:: 1.11 

680 **kwargs 

681 Keyword arguments passed to `leastsq` for ``method='lm'`` or 

682 `least_squares` otherwise. 

683 

684 Returns 

685 ------- 

686 popt : array 

687 Optimal values for the parameters so that the sum of the squared 

688 residuals of ``f(xdata, *popt) - ydata`` is minimized. 

689 pcov : 2-D array 

690 The estimated approximate covariance of popt. The diagonals provide 

691 the variance of the parameter estimate. To compute one standard 

692 deviation errors on the parameters, use 

693 ``perr = np.sqrt(np.diag(pcov))``. Note that the relationship between 

694 `cov` and parameter error estimates is derived based on a linear 

695 approximation to the model function around the optimum [1]. 

696 When this approximation becomes inaccurate, `cov` may not provide an 

697 accurate measure of uncertainty. 

698 

699 How the `sigma` parameter affects the estimated covariance 

700 depends on `absolute_sigma` argument, as described above. 

701 

702 If the Jacobian matrix at the solution doesn't have a full rank, then 

703 'lm' method returns a matrix filled with ``np.inf``, on the other hand 

704 'trf' and 'dogbox' methods use Moore-Penrose pseudoinverse to compute 

705 the covariance matrix. Covariance matrices with large condition numbers 

706 (e.g. computed with `numpy.linalg.cond`) may indicate that results are 

707 unreliable. 

708 infodict : dict (returned only if `full_output` is True) 

709 a dictionary of optional outputs with the keys: 

710 

711 ``nfev`` 

712 The number of function calls. Methods 'trf' and 'dogbox' do not 

713 count function calls for numerical Jacobian approximation, 

714 as opposed to 'lm' method. 

715 ``fvec`` 

716 The residual values evaluated at the solution, for a 1-D `sigma` 

717 this is ``(f(x, *popt) - ydata)/sigma``. 

718 ``fjac`` 

719 A permutation of the R matrix of a QR 

720 factorization of the final approximate 

721 Jacobian matrix, stored column wise. 

722 Together with ipvt, the covariance of the 

723 estimate can be approximated. 

724 Method 'lm' only provides this information. 

725 ``ipvt`` 

726 An integer array of length N which defines 

727 a permutation matrix, p, such that 

728 fjac*p = q*r, where r is upper triangular 

729 with diagonal elements of nonincreasing 

730 magnitude. Column j of p is column ipvt(j) 

731 of the identity matrix. 

732 Method 'lm' only provides this information. 

733 ``qtf`` 

734 The vector (transpose(q) * fvec). 

735 Method 'lm' only provides this information. 

736 

737 .. versionadded:: 1.9 

738 mesg : str (returned only if `full_output` is True) 

739 A string message giving information about the solution. 

740 

741 .. versionadded:: 1.9 

742 ier : int (returnned only if `full_output` is True) 

743 An integer flag. If it is equal to 1, 2, 3 or 4, the solution was 

744 found. Otherwise, the solution was not found. In either case, the 

745 optional output variable `mesg` gives more information. 

746 

747 .. versionadded:: 1.9 

748 

749 Raises 

750 ------ 

751 ValueError 

752 if either `ydata` or `xdata` contain NaNs, or if incompatible options 

753 are used. 

754 

755 RuntimeError 

756 if the least-squares minimization fails. 

757 

758 OptimizeWarning 

759 if covariance of the parameters can not be estimated. 

760 

761 See Also 

762 -------- 

763 least_squares : Minimize the sum of squares of nonlinear functions. 

764 scipy.stats.linregress : Calculate a linear least squares regression for 

765 two sets of measurements. 

766 

767 Notes 

768 ----- 

769 Users should ensure that inputs `xdata`, `ydata`, and the output of `f` 

770 are ``float64``, or else the optimization may return incorrect results. 

771 

772 With ``method='lm'``, the algorithm uses the Levenberg-Marquardt algorithm 

773 through `leastsq`. Note that this algorithm can only deal with 

774 unconstrained problems. 

775 

776 Box constraints can be handled by methods 'trf' and 'dogbox'. Refer to 

777 the docstring of `least_squares` for more information. 

778 

779 References 

780 ---------- 

781 [1] K. Vugrin et al. Confidence region estimation techniques for nonlinear 

782 regression in groundwater flow: Three case studies. Water Resources 

783 Research, Vol. 43, W03423, :doi:`10.1029/2005WR004804` 

784 

785 Examples 

786 -------- 

787 >>> import numpy as np 

788 >>> import matplotlib.pyplot as plt 

789 >>> from scipy.optimize import curve_fit 

790 

791 >>> def func(x, a, b, c): 

792 ... return a * np.exp(-b * x) + c 

793 

794 Define the data to be fit with some noise: 

795 

796 >>> xdata = np.linspace(0, 4, 50) 

797 >>> y = func(xdata, 2.5, 1.3, 0.5) 

798 >>> rng = np.random.default_rng() 

799 >>> y_noise = 0.2 * rng.normal(size=xdata.size) 

800 >>> ydata = y + y_noise 

801 >>> plt.plot(xdata, ydata, 'b-', label='data') 

802 

803 Fit for the parameters a, b, c of the function `func`: 

804 

805 >>> popt, pcov = curve_fit(func, xdata, ydata) 

806 >>> popt 

807 array([2.56274217, 1.37268521, 0.47427475]) 

808 >>> plt.plot(xdata, func(xdata, *popt), 'r-', 

809 ... label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt)) 

810 

811 Constrain the optimization to the region of ``0 <= a <= 3``, 

812 ``0 <= b <= 1`` and ``0 <= c <= 0.5``: 

813 

814 >>> popt, pcov = curve_fit(func, xdata, ydata, bounds=(0, [3., 1., 0.5])) 

815 >>> popt 

816 array([2.43736712, 1. , 0.34463856]) 

817 >>> plt.plot(xdata, func(xdata, *popt), 'g--', 

818 ... label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt)) 

819 

820 >>> plt.xlabel('x') 

821 >>> plt.ylabel('y') 

822 >>> plt.legend() 

823 >>> plt.show() 

824 

825 For reliable results, the model `func` should not be overparametrized; 

826 redundant parameters can cause unreliable covariance matrices and, in some 

827 cases, poorer quality fits. As a quick check of whether the model may be 

828 overparameterized, calculate the condition number of the covariance matrix: 

829 

830 >>> np.linalg.cond(pcov) 

831 34.571092161547405 # may vary 

832 

833 The value is small, so it does not raise much concern. If, however, we were 

834 to add a fourth parameter ``d`` to `func` with the same effect as ``a``: 

835 

836 >>> def func(x, a, b, c, d): 

837 ... return a * d * np.exp(-b * x) + c # a and d are redundant 

838 >>> popt, pcov = curve_fit(func, xdata, ydata) 

839 >>> np.linalg.cond(pcov) 

840 1.13250718925596e+32 # may vary 

841 

842 Such a large value is cause for concern. The diagonal elements of the 

843 covariance matrix, which is related to uncertainty of the fit, gives more 

844 information: 

845 

846 >>> np.diag(pcov) 

847 array([1.48814742e+29, 3.78596560e-02, 5.39253738e-03, 2.76417220e+28]) # may vary 

848 

849 Note that the first and last terms are much larger than the other elements, 

850 suggesting that the optimal values of these parameters are ambiguous and 

851 that only one of these parameters is needed in the model. 

852 

853 """ # noqa 

854 if p0 is None: 

855 # determine number of parameters by inspecting the function 

856 sig = _getfullargspec(f) 

857 args = sig.args 

858 if len(args) < 2: 

859 raise ValueError("Unable to determine number of fit parameters.") 

860 n = len(args) - 1 

861 else: 

862 p0 = np.atleast_1d(p0) 

863 n = p0.size 

864 

865 if isinstance(bounds, Bounds): 

866 lb, ub = bounds.lb, bounds.ub 

867 else: 

868 lb, ub = prepare_bounds(bounds, n) 

869 if p0 is None: 

870 p0 = _initialize_feasible(lb, ub) 

871 

872 bounded_problem = np.any((lb > -np.inf) | (ub < np.inf)) 

873 if method is None: 

874 if bounded_problem: 

875 method = 'trf' 

876 else: 

877 method = 'lm' 

878 

879 if method == 'lm' and bounded_problem: 

880 raise ValueError("Method 'lm' only works for unconstrained problems. " 

881 "Use 'trf' or 'dogbox' instead.") 

882 

883 if check_finite is None: 

884 check_finite = True if nan_policy is None else False 

885 

886 # optimization may produce garbage for float32 inputs, cast them to float64 

887 if check_finite: 

888 ydata = np.asarray_chkfinite(ydata, float) 

889 else: 

890 ydata = np.asarray(ydata, float) 

891 

892 if isinstance(xdata, (list, tuple, np.ndarray)): 

893 # `xdata` is passed straight to the user-defined `f`, so allow 

894 # non-array_like `xdata`. 

895 if check_finite: 

896 xdata = np.asarray_chkfinite(xdata, float) 

897 else: 

898 xdata = np.asarray(xdata, float) 

899 

900 if ydata.size == 0: 

901 raise ValueError("`ydata` must not be empty!") 

902 

903 # nan handling is needed only if check_finite is False because if True, 

904 # the x-y data are already checked, and they don't contain nans. 

905 if not check_finite and nan_policy is not None: 

906 if nan_policy == "propagate": 

907 raise ValueError("`nan_policy='propagate'` is not supported " 

908 "by this function.") 

909 

910 policies = [None, 'raise', 'omit'] 

911 x_contains_nan, nan_policy = _contains_nan(xdata, nan_policy, 

912 policies=policies) 

913 y_contains_nan, nan_policy = _contains_nan(ydata, nan_policy, 

914 policies=policies) 

915 

916 if (x_contains_nan or y_contains_nan) and nan_policy == 'omit': 

917 # ignore NaNs for N dimensional arrays 

918 has_nan = np.isnan(xdata) 

919 has_nan = has_nan.any(axis=tuple(range(has_nan.ndim-1))) 

920 has_nan |= np.isnan(ydata) 

921 

922 xdata = xdata[..., ~has_nan] 

923 ydata = ydata[~has_nan] 

924 

925 # Determine type of sigma 

926 if sigma is not None: 

927 sigma = np.asarray(sigma) 

928 

929 # if 1-D, sigma are errors, define transform = 1/sigma 

930 if sigma.shape == (ydata.size, ): 

931 transform = 1.0 / sigma 

932 # if 2-D, sigma is the covariance matrix, 

933 # define transform = L such that L L^T = C 

934 elif sigma.shape == (ydata.size, ydata.size): 

935 try: 

936 # scipy.linalg.cholesky requires lower=True to return L L^T = A 

937 transform = cholesky(sigma, lower=True) 

938 except LinAlgError as e: 

939 raise ValueError("`sigma` must be positive definite.") from e 

940 else: 

941 raise ValueError("`sigma` has incorrect shape.") 

942 else: 

943 transform = None 

944 

945 func = _lightweight_memoizer(_wrap_func(f, xdata, ydata, transform)) 

946 

947 if callable(jac): 

948 jac = _lightweight_memoizer(_wrap_jac(jac, xdata, transform)) 

949 elif jac is None and method != 'lm': 

950 jac = '2-point' 

951 

952 if 'args' in kwargs: 

953 # The specification for the model function `f` does not support 

954 # additional arguments. Refer to the `curve_fit` docstring for 

955 # acceptable call signatures of `f`. 

956 raise ValueError("'args' is not a supported keyword argument.") 

957 

958 if method == 'lm': 

959 # if ydata.size == 1, this might be used for broadcast. 

960 if ydata.size != 1 and n > ydata.size: 

961 raise TypeError(f"The number of func parameters={n} must not" 

962 f" exceed the number of data points={ydata.size}") 

963 res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs) 

964 popt, pcov, infodict, errmsg, ier = res 

965 ysize = len(infodict['fvec']) 

966 cost = np.sum(infodict['fvec'] ** 2) 

967 if ier not in [1, 2, 3, 4]: 

968 raise RuntimeError("Optimal parameters not found: " + errmsg) 

969 else: 

970 # Rename maxfev (leastsq) to max_nfev (least_squares), if specified. 

971 if 'max_nfev' not in kwargs: 

972 kwargs['max_nfev'] = kwargs.pop('maxfev', None) 

973 

974 res = least_squares(func, p0, jac=jac, bounds=bounds, method=method, 

975 **kwargs) 

976 

977 if not res.success: 

978 raise RuntimeError("Optimal parameters not found: " + res.message) 

979 

980 infodict = dict(nfev=res.nfev, fvec=res.fun)