1
2
3 """
4 Perform Levenberg-Marquardt least-squares minimization, based on MINPACK-1.
5
6 AUTHORS
7 The original version of this software, called LMFIT, was written in FORTRAN
8 as part of the MINPACK-1 package by XXX.
9
10 Craig Markwardt converted the FORTRAN code to IDL. The information for the
11 IDL version is:
12 Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
13 craigm@lheamail.gsfc.nasa.gov
14 UPDATED VERSIONs can be found on my WEB PAGE:
15 http://cow.physics.wisc.edu/~craigm/idl/idl.html
16
17 Mark Rivers created this Python version from Craig's IDL version.
18 Mark Rivers, University of Chicago
19 Building 434A, Argonne National Laboratory
20 9700 South Cass Avenue, Argonne, IL 60439
21 rivers@cars.uchicago.edu
22 Updated versions can be found at http://cars.uchicago.edu/software
23
24 Sergey Koposov converted the Mark's Python version from Numeric to numpy
25 Sergey Koposov, Max Planck Institute for Astronomy
26 Heidelberg, Germany, D-69117
27 koposov@mpia.de
28 Updated versions can be found at http://code.google.com/p/astrolibpy/source/browse/trunk/
29
30 DESCRIPTION
31
32 MPFIT uses the Levenberg-Marquardt technique to solve the
33 least-squares problem. In its typical use, MPFIT will be used to
34 fit a user-supplied function (the "model") to user-supplied data
35 points (the "data") by adjusting a set of parameters. MPFIT is
36 based upon MINPACK-1 (LMDIF.F) by More' and collaborators.
37
38 For example, a researcher may think that a set of observed data
39 points is best modelled with a Gaussian curve. A Gaussian curve is
40 parameterized by its mean, standard deviation and normalization.
41 MPFIT will, within certain constraints, find the set of parameters
42 which best fits the data. The fit is "best" in the least-squares
43 sense; that is, the sum of the weighted squared differences between
44 the model and data is minimized.
45
46 The Levenberg-Marquardt technique is a particular strategy for
47 iteratively searching for the best fit. This particular
48 implementation is drawn from MINPACK-1 (see NETLIB), and is much faster
49 and more accurate than the version provided in the Scientific Python package
50 in Scientific.Functions.LeastSquares.
51 This version allows upper and lower bounding constraints to be placed on each
52 parameter, or the parameter can be held fixed.
53
54 The user-supplied Python function should return an array of weighted
55 deviations between model and data. In a typical scientific problem
56 the residuals should be weighted so that each deviate has a
57 gaussian sigma of 1.0. If X represents values of the independent
58 variable, Y represents a measurement for each value of X, and ERR
59 represents the error in the measurements, then the deviates could
60 be calculated as follows:
61
62 DEVIATES = (Y - F(X)) / ERR
63
64 where F is the analytical function representing the model. You are
65 recommended to use the convenience functions MPFITFUN and
66 MPFITEXPR, which are driver functions that calculate the deviates
67 for you. If ERR are the 1-sigma uncertainties in Y, then
68
69 TOTAL( DEVIATES^2 )
70
71 will be the total chi-squared value. MPFIT will minimize the
72 chi-square value. The values of X, Y and ERR are passed through
73 MPFIT to the user-supplied function via the FUNCTKW keyword.
74
75 Simple constraints can be placed on parameter values by using the
76 PARINFO keyword to MPFIT. See below for a description of this
77 keyword.
78
79 MPFIT does not perform more general optimization tasks. See TNMIN
80 instead. MPFIT is customized, based on MINPACK-1, to the
81 least-squares minimization problem.
82
83
84 USER FUNCTION
85
86 The user must define a function which returns the appropriate
87 values as specified above. The function should return the weighted
88 deviations between the model and the data. It should also return a status
89 flag and an optional partial derivative array. For applications which
90 use finite-difference derivatives -- the default -- the user
91 function should be declared in the following way:
92
93 def myfunct(p, fjac=None, x=None, y=None, err=None)
94 # Parameter values are passed in "p"
95 # If fjac==None then partial derivatives should not be
96 # computed. It will always be None if MPFIT is called with default
97 # flag.
98 model = F(x, p)
99 # Non-negative status value means MPFIT should continue, negative means
100 # stop the calculation.
101 status = 0
102 return([status, (y-model)/err]
103
104 See below for applications with analytical derivatives.
105
106 The keyword parameters X, Y, and ERR in the example above are
107 suggestive but not required. Any parameters can be passed to
108 MYFUNCT by using the functkw keyword to MPFIT. Use MPFITFUN and
109 MPFITEXPR if you need ideas on how to do that. The function *must*
110 accept a parameter list, P.
111
112 In general there are no restrictions on the number of dimensions in
113 X, Y or ERR. However the deviates *must* be returned in a
114 one-dimensional Numeric array of type Float.
115
116 User functions may also indicate a fatal error condition using the
117 status return described above. If status is set to a number between
118 -15 and -1 then MPFIT will stop the calculation and return to the caller.
119
120
121 ANALYTIC DERIVATIVES
122
123 In the search for the best-fit solution, MPFIT by default
124 calculates derivatives numerically via a finite difference
125 approximation. The user-supplied function need not calculate the
126 derivatives explicitly. However, if you desire to compute them
127 analytically, then the AUTODERIVATIVE=0 keyword must be passed to MPFIT.
128 As a practical matter, it is often sufficient and even faster to allow
129 MPFIT to calculate the derivatives numerically, and so
130 AUTODERIVATIVE=0 is not necessary.
131
132 If AUTODERIVATIVE=0 is used then the user function must check the parameter
133 FJAC, and if FJAC!=None then return the partial derivative array in the
134 return list.
135 def myfunct(p, fjac=None, x=None, y=None, err=None)
136 # Parameter values are passed in "p"
137 # If FJAC!=None then partial derivatives must be comptuer.
138 # FJAC contains an array of len(p), where each entry
139 # is 1 if that parameter is free and 0 if it is fixed.
140 model = F(x, p)
141 Non-negative status value means MPFIT should continue, negative means
142 # stop the calculation.
143 status = 0
144 if (dojac):
145 pderiv = zeros([len(x), len(p)], Float)
146 for j in range(len(p)):
147 pderiv[:,j] = FGRAD(x, p, j)
148 else:
149 pderiv = None
150 return([status, (y-model)/err, pderiv]
151
152 where FGRAD(x, p, i) is a user function which must compute the
153 derivative of the model with respect to parameter P[i] at X. When
154 finite differencing is used for computing derivatives (ie, when
155 AUTODERIVATIVE=1), or when MPFIT needs only the errors but not the
156 derivatives the parameter FJAC=None.
157
158 Derivatives should be returned in the PDERIV array. PDERIV should be an m x
159 n array, where m is the number of data points and n is the number
160 of parameters. dp[i,j] is the derivative at the ith point with
161 respect to the jth parameter.
162
163 The derivatives with respect to fixed parameters are ignored; zero
164 is an appropriate value to insert for those derivatives. Upon
165 input to the user function, FJAC is set to a vector with the same
166 length as P, with a value of 1 for a parameter which is free, and a
167 value of zero for a parameter which is fixed (and hence no
168 derivative needs to be calculated).
169
170 If the data is higher than one dimensional, then the *last*
171 dimension should be the parameter dimension. Example: fitting a
172 50x50 image, "dp" should be 50x50xNPAR.
173
174
175 CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD
176
177 The behavior of MPFIT can be modified with respect to each
178 parameter to be fitted. A parameter value can be fixed; simple
179 boundary constraints can be imposed; limitations on the parameter
180 changes can be imposed; properties of the automatic derivative can
181 be modified; and parameters can be tied to one another.
182
183 These properties are governed by the PARINFO structure, which is
184 passed as a keyword parameter to MPFIT.
185
186 PARINFO should be a list of dictionaries, one list entry for each parameter.
187 Each parameter is associated with one element of the array, in
188 numerical order. The dictionary can have the following keys
189 (none are required, keys are case insensitive):
190
191 'value' - the starting parameter value (but see the START_PARAMS
192 parameter for more information).
193
194 'fixed' - a boolean value, whether the parameter is to be held
195 fixed or not. Fixed parameters are not varied by
196 MPFIT, but are passed on to MYFUNCT for evaluation.
197
198 'limited' - a two-element boolean array. If the first/second
199 element is set, then the parameter is bounded on the
200 lower/upper side. A parameter can be bounded on both
201 sides. Both LIMITED and LIMITS must be given
202 together.
203
204 'limits' - a two-element float array. Gives the
205 parameter limits on the lower and upper sides,
206 respectively. Zero, one or two of these values can be
207 set, depending on the values of LIMITED. Both LIMITED
208 and LIMITS must be given together.
209
210 'parname' - a string, giving the name of the parameter. The
211 fitting code of MPFIT does not use this tag in any
212 way. However, the default iterfunct will print the
213 parameter name if available.
214
215 'step' - the step size to be used in calculating the numerical
216 derivatives. If set to zero, then the step size is
217 computed automatically. Ignored when AUTODERIVATIVE=0.
218
219 'mpside' - the sidedness of the finite difference when computing
220 numerical derivatives. This field can take four
221 values:
222
223 0 - one-sided derivative computed automatically
224 1 - one-sided derivative (f(x+h) - f(x) )/h
225 -1 - one-sided derivative (f(x) - f(x-h))/h
226 2 - two-sided derivative (f(x+h) - f(x-h))/(2*h)
227
228 Where H is the STEP parameter described above. The
229 "automatic" one-sided derivative method will chose a
230 direction for the finite difference which does not
231 violate any constraints. The other methods do not
232 perform this check. The two-sided method is in
233 principle more precise, but requires twice as many
234 function evaluations. Default: 0.
235
236 'mpmaxstep' - the maximum change to be made in the parameter
237 value. During the fitting process, the parameter
238 will never be changed by more than this value in
239 one iteration.
240
241 A value of 0 indicates no maximum. Default: 0.
242
243 'tied' - a string expression which "ties" the parameter to other
244 free or fixed parameters. Any expression involving
245 constants and the parameter array P are permitted.
246 Example: if parameter 2 is always to be twice parameter
247 1 then use the following: parinfo(2).tied = '2 * p(1)'.
248 Since they are totally constrained, tied parameters are
249 considered to be fixed; no errors are computed for them.
250 [ NOTE: the PARNAME can't be used in expressions. ]
251
252 'mpprint' - if set to 1, then the default iterfunct will print the
253 parameter value. If set to 0, the parameter value
254 will not be printed. This tag can be used to
255 selectively print only a few parameter values out of
256 many. Default: 1 (all parameters printed)
257
258
259 Future modifications to the PARINFO structure, if any, will involve
260 adding dictionary tags beginning with the two letters "MP".
261 Therefore programmers are urged to avoid using tags starting with
262 the same letters; otherwise they are free to include their own
263 fields within the PARINFO structure, and they will be ignored.
264
265 PARINFO Example:
266 parinfo = [{'value':0., 'fixed':0, 'limited':[0,0], 'limits':[0.,0.]}]*5
267 parinfo[0]['fixed'] = 1
268 parinfo[4]['limited'][0] = 1
269 parinfo[4]['limits'][0] = 50.
270 values = [5.7, 2.2, 500., 1.5, 2000.]
271 for i in range(5): parinfo[i]['value']=values[i]
272
273 A total of 5 parameters, with starting values of 5.7,
274 2.2, 500, 1.5, and 2000 are given. The first parameter
275 is fixed at a value of 5.7, and the last parameter is
276 constrained to be above 50.
277
278
279 EXAMPLE
280
281 import mpfit
282 import numpy.oldnumeric as Numeric
283 x = arange(100, float)
284 p0 = [5.7, 2.2, 500., 1.5, 2000.]
285 y = ( p[0] + p[1]*[x] + p[2]*[x**2] + p[3]*sqrt(x) +
286 p[4]*log(x))
287 fa = {'x':x, 'y':y, 'err':err}
288 m = mpfit('myfunct', p0, functkw=fa)
289 print 'status = ', m.status
290 if (m.status <= 0): print 'error message = ', m.errmsg
291 print 'parameters = ', m.params
292
293 Minimizes sum of squares of MYFUNCT. MYFUNCT is called with the X,
294 Y, and ERR keyword parameters that are given by FUNCTKW. The
295 results can be obtained from the returned object m.
296
297
298 THEORY OF OPERATION
299
300 There are many specific strategies for function minimization. One
301 very popular technique is to use function gradient information to
302 realize the local structure of the function. Near a local minimum
303 the function value can be taylor expanded about x0 as follows:
304
305 f(x) = f(x0) + f'(x0) . (x-x0) + (1/2) (x-x0) . f''(x0) . (x-x0)
306 ----- --------------- ------------------------------- (1)
307 Order 0th 1st 2nd
308
309 Here f'(x) is the gradient vector of f at x, and f''(x) is the
310 Hessian matrix of second derivatives of f at x. The vector x is
311 the set of function parameters, not the measured data vector. One
312 can find the minimum of f, f(xm) using Newton's method, and
313 arrives at the following linear equation:
314
315 f''(x0) . (xm-x0) = - f'(x0) (2)
316
317 If an inverse can be found for f''(x0) then one can solve for
318 (xm-x0), the step vector from the current position x0 to the new
319 projected minimum. Here the problem has been linearized (ie, the
320 gradient information is known to first order). f''(x0) is
321 symmetric n x n matrix, and should be positive definite.
322
323 The Levenberg - Marquardt technique is a variation on this theme.
324 It adds an additional diagonal term to the equation which may aid the
325 convergence properties:
326
327 (f''(x0) + nu I) . (xm-x0) = -f'(x0) (2a)
328
329 where I is the identity matrix. When nu is large, the overall
330 matrix is diagonally dominant, and the iterations follow steepest
331 descent. When nu is small, the iterations are quadratically
332 convergent.
333
334 In principle, if f''(x0) and f'(x0) are known then xm-x0 can be
335 determined. However the Hessian matrix is often difficult or
336 impossible to compute. The gradient f'(x0) may be easier to
337 compute, if even by finite difference techniques. So-called
338 quasi-Newton techniques attempt to successively estimate f''(x0)
339 by building up gradient information as the iterations proceed.
340
341 In the least squares problem there are further simplifications
342 which assist in solving eqn (2). The function to be minimized is
343 a sum of squares:
344
345 f = Sum(hi^2) (3)
346
347 where hi is the ith residual out of m residuals as described
348 above. This can be substituted back into eqn (2) after computing
349 the derivatives:
350
351 f' = 2 Sum(hi hi')
352 f'' = 2 Sum(hi' hj') + 2 Sum(hi hi'') (4)
353
354 If one assumes that the parameters are already close enough to a
355 minimum, then one typically finds that the second term in f'' is
356 negligible [or, in any case, is too difficult to compute]. Thus,
357 equation (2) can be solved, at least approximately, using only
358 gradient information.
359
360 In matrix notation, the combination of eqns (2) and (4) becomes:
361
362 hT' . h' . dx = - hT' . h (5)
363
364 Where h is the residual vector (length m), hT is its transpose, h'
365 is the Jacobian matrix (dimensions n x m), and dx is (xm-x0). The
366 user function supplies the residual vector h, and in some cases h'
367 when it is not found by finite differences (see MPFIT_FDJAC2,
368 which finds h and hT'). Even if dx is not the best absolute step
369 to take, it does provide a good estimate of the best *direction*,
370 so often a line minimization will occur along the dx vector
371 direction.
372
373 The method of solution employed by MINPACK is to form the Q . R
374 factorization of h', where Q is an orthogonal matrix such that QT .
375 Q = I, and R is upper right triangular. Using h' = Q . R and the
376 ortogonality of Q, eqn (5) becomes
377
378 (RT . QT) . (Q . R) . dx = - (RT . QT) . h
379 RT . R . dx = - RT . QT . h (6)
380 R . dx = - QT . h
381
382 where the last statement follows because R is upper triangular.
383 Here, R, QT and h are known so this is a matter of solving for dx.
384 The routine MPFIT_QRFAC provides the QR factorization of h, with
385 pivoting, and MPFIT_QRSOLV provides the solution for dx.
386
387
388 REFERENCES
389
390 MINPACK-1, Jorge More', available from netlib (www.netlib.org).
391 "Optimization Software Guide," Jorge More' and Stephen Wright,
392 SIAM, *Frontiers in Applied Mathematics*, Number 14.
393 More', Jorge J., "The Levenberg-Marquardt Algorithm:
394 Implementation and Theory," in *Numerical Analysis*, ed. Watson,
395 G. A., Lecture Notes in Mathematics 630, Springer-Verlag, 1977.
396
397
398 MODIFICATION HISTORY
399
400 Translated from MINPACK-1 in FORTRAN, Apr-Jul 1998, CM
401 Copyright (C) 1997-2002, Craig Markwardt
402 This software is provided as is without any warranty whatsoever.
403 Permission to use, copy, modify, and distribute modified or
404 unmodified copies is granted, provided this copyright and disclaimer
405 are included unchanged.
406
407 Translated from MPFIT (Craig Markwardt's IDL package) to Python,
408 August, 2002. Mark Rivers
409 Converted from Numeric to numpy (Sergey Koposov, July 2008)
410 """
411
412 from numpy import *
413 import numpy
414 import types
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
599 - def __init__(self, fcn, xall=None, functkw=None, parinfo=None,
600 ftol=1.e-10, xtol=1.e-10, gtol=1.e-10,
601 damp=0., maxiter=200, factor=100., nprint=1,
602 iterfunct='default', iterkw=None, nocovar=0,
603 fastnorm=0, rescale=0, autoderivative=1, quiet=0,
604 diag=None, epsfcn=None, debug=0):
605 """
606 Inputs:
607 fcn:
608 The function to be minimized. The function should return the weighted
609 deviations between the model and the data, as described above.
610
611 xall:
612 An array of starting values for each of the parameters of the model.
613 The number of parameters should be fewer than the number of measurements.
614
615 This parameter is optional if the parinfo keyword is used (but see
616 parinfo). The parinfo keyword provides a mechanism to fix or constrain
617 individual parameters.
618
619 Keywords:
620
621 autoderivative:
622 If this is set, derivatives of the function will be computed
623 automatically via a finite differencing procedure. If not set, then
624 fcn must provide the (analytical) derivatives.
625 Default: set (=1)
626 NOTE: to supply your own analytical derivatives,
627 explicitly pass autoderivative=0
628
629 fastnorm:
630 Set this keyword to select a faster algorithm to compute sum-of-square
631 values internally. For systems with large numbers of data points, the
632 standard algorithm can become prohibitively slow because it cannot be
633 vectorized well. By setting this keyword, MPFIT will run faster, but
634 it will be more prone to floating point overflows and underflows. Thus, setting
635 this keyword may sacrifice some stability in the fitting process.
636 Default: clear (=0)
637
638 ftol:
639 A nonnegative input variable. Termination occurs when both the actual
640 and predicted relative reductions in the sum of squares are at most
641 ftol (and status is accordingly set to 1 or 3). Therefore, ftol
642 measures the relative error desired in the sum of squares.
643 Default: 1E-10
644
645 functkw:
646 A dictionary which contains the parameters to be passed to the
647 user-supplied function specified by fcn via the standard Python
648 keyword dictionary mechanism. This is the way you can pass additional
649 data to your user-supplied function without using global variables.
650
651 Consider the following example:
652 if functkw = {'xval':[1.,2.,3.], 'yval':[1.,4.,9.],
653 'errval':[1.,1.,1.] }
654 then the user supplied function should be declared like this:
655 def myfunct(p, fjac=None, xval=None, yval=None, errval=None):
656
657 Default: {} No extra parameters are passed to the user-supplied
658 function.
659
660 gtol:
661 A nonnegative input variable. Termination occurs when the cosine of
662 the angle between fvec and any column of the jacobian is at most gtol
663 in absolute value (and status is accordingly set to 4). Therefore,
664 gtol measures the orthogonality desired between the function vector
665 and the columns of the jacobian.
666 Default: 1e-10
667
668 iterkw:
669 The keyword arguments to be passed to iterfunct via the dictionary
670 keyword mechanism. This should be a dictionary and is similar in
671 operation to FUNCTKW.
672 Default: {} No arguments are passed.
673
674 iterfunct:
675 The name of a function to be called upon each NPRINT iteration of the
676 MPFIT routine. It should be declared in the following way:
677 def iterfunct(myfunct, p, iter, fnorm, functkw=None,
678 parinfo=None, quiet=0, dof=None, [iterkw keywords here])
679 # perform custom iteration update
680
681 iterfunct must accept all three keyword parameters (FUNCTKW, PARINFO
682 and QUIET).
683
684 myfunct: The user-supplied function to be minimized,
685 p: The current set of model parameters
686 iter: The iteration number
687 functkw: The arguments to be passed to myfunct.
688 fnorm: The chi-squared value.
689 quiet: Set when no textual output should be printed.
690 dof: The number of degrees of freedom, normally the number of points
691 less the number of free parameters.
692 See below for documentation of parinfo.
693
694 In implementation, iterfunct can perform updates to the terminal or
695 graphical user interface, to provide feedback while the fit proceeds.
696 If the fit is to be stopped for any reason, then iterfunct should return a
697 a status value between -15 and -1. Otherwise it should return None
698 (e.g. no return statement) or 0.
699 In principle, iterfunct should probably not modify the parameter values,
700 because it may interfere with the algorithm's stability. In practice it
701 is allowed.
702
703 Default: an internal routine is used to print the parameter values.
704
705 Set iterfunct=None if there is no user-defined routine and you don't
706 want the internal default routine be called.
707
708 maxiter:
709 The maximum number of iterations to perform. If the number is exceeded,
710 then the status value is set to 5 and MPFIT returns.
711 Default: 200 iterations
712
713 nocovar:
714 Set this keyword to prevent the calculation of the covariance matrix
715 before returning (see COVAR)
716 Default: clear (=0) The covariance matrix is returned
717
718 nprint:
719 The frequency with which iterfunct is called. A value of 1 indicates
720 that iterfunct is called with every iteration, while 2 indicates every
721 other iteration, etc. Note that several Levenberg-Marquardt attempts
722 can be made in a single iteration.
723 Default value: 1
724
725 parinfo
726 Provides a mechanism for more sophisticated constraints to be placed on
727 parameter values. When parinfo is not passed, then it is assumed that
728 all parameters are free and unconstrained. Values in parinfo are never
729 modified during a call to MPFIT.
730
731 See description above for the structure of PARINFO.
732
733 Default value: None All parameters are free and unconstrained.
734
735 quiet:
736 Set this keyword when no textual output should be printed by MPFIT
737
738 damp:
739 A scalar number, indicating the cut-off value of residuals where
740 "damping" will occur. Residuals with magnitudes greater than this
741 number will be replaced by their hyperbolic tangent. This partially
742 mitigates the so-called large residual problem inherent in
743 least-squares solvers (as for the test problem CURVI,
744 http://www.maxthis.com/curviex.htm).
745 A value of 0 indicates no damping.
746 Default: 0
747
748 Note: DAMP doesn't work with autoderivative=0
749
750 xtol:
751 A nonnegative input variable. Termination occurs when the relative error
752 between two consecutive iterates is at most xtol (and status is
753 accordingly set to 2 or 3). Therefore, xtol measures the relative error
754 desired in the approximate solution.
755 Default: 1E-10
756
757 Outputs:
758
759 Returns an object of type mpfit. The results are attributes of this class,
760 e.g. mpfit.status, mpfit.errmsg, mpfit.params, npfit.niter, mpfit.covar.
761
762 .status
763 An integer status code is returned. All values greater than zero can
764 represent success (however .status == 5 may indicate failure to
765 converge). It can have one of the following values:
766
767 -16
768 A parameter or function value has become infinite or an undefined
769 number. This is usually a consequence of numerical overflow in the
770 user's model function, which must be avoided.
771
772 -15 to -1
773 These are error codes that either MYFUNCT or iterfunct may return to
774 terminate the fitting process. Values from -15 to -1 are reserved
775 for the user functions and will not clash with MPFIT.
776
777 0 Improper input parameters.
778
779 1 Both actual and predicted relative reductions in the sum of squares
780 are at most ftol.
781
782 2 Relative error between two consecutive iterates is at most xtol
783
784 3 Conditions for status = 1 and status = 2 both hold.
785
786 4 The cosine of the angle between fvec and any column of the jacobian
787 is at most gtol in absolute value.
788
789 5 The maximum number of iterations has been reached.
790
791 6 ftol is too small. No further reduction in the sum of squares is
792 possible.
793
794 7 xtol is too small. No further improvement in the approximate solution
795 x is possible.
796
797 8 gtol is too small. fvec is orthogonal to the columns of the jacobian
798 to machine precision.
799
800 .fnorm
801 The value of the summed squared residuals for the returned parameter
802 values.
803
804 .covar
805 The covariance matrix for the set of parameters returned by MPFIT.
806 The matrix is NxN where N is the number of parameters. The square root
807 of the diagonal elements gives the formal 1-sigma statistical errors on
808 the parameters if errors were treated "properly" in fcn.
809 Parameter errors are also returned in .perror.
810
811 To compute the correlation matrix, pcor, use this example:
812 cov = mpfit.covar
813 pcor = cov * 0.
814 for i in range(n):
815 for j in range(n):
816 pcor[i,j] = cov[i,j]/sqrt(cov[i,i]*cov[j,j])
817
818 If nocovar is set or MPFIT terminated abnormally, then .covar is set to
819 a scalar with value None.
820
821 .errmsg
822 A string error or warning message is returned.
823
824 .nfev
825 The number of calls to MYFUNCT performed.
826
827 .niter
828 The number of iterations completed.
829
830 .perror
831 The formal 1-sigma errors in each parameter, computed from the
832 covariance matrix. If a parameter is held fixed, or if it touches a
833 boundary, then the error is reported as zero.
834
835 If the fit is unweighted (i.e. no errors were given, or the weights
836 were uniformly set to unity), then .perror will probably not represent
837 the true parameter uncertainties.
838
839 *If* you can assume that the true reduced chi-squared value is unity --
840 meaning that the fit is implicitly assumed to be of good quality --
841 then the estimated parameter uncertainties can be computed by scaling
842 .perror by the measured chi-squared value.
843
844 dof = len(x) - len(mpfit.params) # deg of freedom
845 # scaled uncertainties
846 pcerror = mpfit.perror * sqrt(mpfit.fnorm / dof)
847
848 """
849
850
851 if functkw is None:
852 functkw = {}
853 if iterkw is None:
854 iterkw = {}
855
856 self.niter = 0
857 self.params = None
858 self.covar = None
859 self.perror = None
860 self.status = 0
861 self.debug = debug
862 self.errmsg = ''
863 self.fastnorm = fastnorm
864 self.nfev = 0
865 self.damp = damp
866 self.machar = machar(double=1)
867 machep = self.machar.machep
868 self.dof=0
869
870 if (fcn==None):
871 self.errmsg = "Usage: parms = mpfit('myfunt', ... )"
872 return
873
874 if (iterfunct == 'default'): iterfunct = self.defiter
875
876
877
878 if (self.damp != 0) and (autoderivative == 0):
879 self.errmsg = 'ERROR: keywords DAMP and AUTODERIVATIVE are mutually exclusive'
880 return
881
882
883 if (xall == None) and (parinfo == None):
884 self.errmsg = 'ERROR: must pass parameters in P or PARINFO'
885 return
886
887
888 if (parinfo != None):
889 if (type(parinfo) != types.ListType):
890 self.errmsg = 'ERROR: PARINFO must be a list of dictionaries.'
891 return
892 else:
893 if (type(parinfo[0]) != types.DictionaryType):
894 self.errmsg = 'ERROR: PARINFO must be a list of dictionaries.'
895 return
896 if ((xall != None) and (len(xall) != len(parinfo))):
897 self.errmsg = 'ERROR: number of elements in PARINFO and P must agree'
898 return
899
900
901
902 if (xall == None):
903 xall = self.parinfo(parinfo, 'value')
904 if (xall == None):
905 self.errmsg = 'ERROR: either P or PARINFO(*)["value"] must be supplied.'
906 return
907
908
909 xall = asarray(xall, float)
910
911 npar = len(xall)
912 self.fnorm = -1.
913 fnorm1 = -1.
914
915
916 ptied = self.parinfo(parinfo, 'tied', default='', n=npar)
917 self.qanytied = 0
918 for i in range(npar):
919 ptied[i] = ptied[i].strip()
920 if (ptied[i] != ''): self.qanytied = 1
921 self.ptied = ptied
922
923
924 pfixed = self.parinfo(parinfo, 'fixed', default=0, n=npar)
925 pfixed = (pfixed == 1)
926 for i in range(npar):
927 pfixed[i] = pfixed[i] or (ptied[i] != '')
928
929
930 step = self.parinfo(parinfo, 'step', default=0., n=npar)
931 dstep = self.parinfo(parinfo, 'relstep', default=0., n=npar)
932 dside = self.parinfo(parinfo, 'mpside', default=0, n=npar)
933
934
935 maxstep = self.parinfo(parinfo, 'mpmaxstep', default=0., n=npar)
936 minstep = self.parinfo(parinfo, 'mpminstep', default=0., n=npar)
937 qmin = minstep * 0
938 qmax = maxstep != 0
939 wh = (nonzero(((qmin!=0.) & (qmax!=0.)) & (maxstep < minstep)))[0]
940 if (len(wh) > 0):
941 self.errmsg = 'ERROR: MPMINSTEP is greater than MPMAXSTEP'
942 return
943 wh = (nonzero((qmin!=0.) | (qmax!=0.)))[0]
944 qminmax = len(wh > 0)
945
946
947 ifree = (nonzero(pfixed != 1))[0]
948 nfree = len(ifree)
949 if nfree == 0:
950 self.errmsg = 'ERROR: no free parameters'
951 return
952
953
954 self.params = xall
955 x = take(self.params, ifree)
956
957
958 limited = self.parinfo(parinfo, 'limited', default=[0,0], n=npar)
959 limits = self.parinfo(parinfo, 'limits', default=[0.,0.], n=npar)
960 if (limited != None) and (limits != None):
961
962 wh = (nonzero((limited[:,0] & (xall < limits[:,0])) |
963 (limited[:,1] & (xall > limits[:,1]))))[0]
964 if (len(wh) > 0):
965 self.errmsg = 'ERROR: parameters are not within PARINFO limits'
966 return
967 wh = (nonzero((limited[:,0] & limited[:,1]) &
968 (limits[:,0] >= limits[:,1]) &
969 (pfixed == 0)))[0]
970 if (len(wh) > 0):
971 self.errmsg = 'ERROR: PARINFO parameter limits are not consistent'
972 return
973
974
975 qulim = take(limited[:,1], ifree)
976 ulim = take(limits [:,1], ifree)
977 qllim = take(limited[:,0], ifree)
978 llim = take(limits [:,0], ifree)
979
980 wh = (nonzero((qulim!=0.) | (qllim!=0.)))[0]
981 if (len(wh) > 0): qanylim = 1
982 else: qanylim = 0
983 else:
984
985 qulim = zeros(nfree)
986 ulim = x * 0.
987 qllim = qulim
988 llim = x * 0.
989 qanylim = 0
990
991 n = len(x)
992
993 if ((n < 0) or (ftol <= 0) or (xtol <= 0) or (gtol <= 0)
994 or (maxiter < 0) or (factor <= 0)):
995 self.errmsg = 'ERROR: input keywords are inconsistent'
996 return
997
998 if (rescale != 0):
999 self.errmsg = 'ERROR: DIAG parameter scales are inconsistent'
1000 if (len(diag) < n): return
1001 wh = (nonzero(diag <= 0))[0]
1002 if (len(wh) > 0): return
1003 self.errmsg = ''
1004
1005
1006 x = asarray(x, float)
1007
1008 [self.status, fvec] = self.call(fcn, self.params, functkw)
1009 if (self.status < 0):
1010 self.errmsg = 'ERROR: first call to "'+str(fcn)+'" failed'
1011 return
1012
1013 m = len(fvec)
1014 if (m < n):
1015 self.errmsg = 'ERROR: number of parameters must not exceed data'
1016 return
1017 self.dof = m-nfree
1018 self.fnorm = self.enorm(fvec)
1019
1020
1021
1022 par = 0.
1023 self.niter = 1
1024 qtf = x * 0.
1025 self.status = 0
1026
1027
1028
1029 while(1):
1030
1031
1032 put(self.params, ifree, x)
1033 if (self.qanytied): self.params = self.tie(self.params, ptied)
1034
1035 if (nprint > 0) and (iterfunct != None):
1036 if (((self.niter-1) % nprint) == 0):
1037 mperr = 0
1038 xnew0 = self.params.copy()
1039
1040 dof = numpy.max([len(fvec) - len(x), 0])
1041 status = iterfunct(fcn, self.params, self.niter, self.fnorm**2,
1042 functkw=functkw, parinfo=parinfo, quiet=quiet,
1043 dof=dof, **iterkw)
1044 if (status != None): self.status = status
1045
1046
1047 if (self.status < 0):
1048 self.errmsg = 'WARNING: premature termination by ' + str(iterfunct)
1049 return
1050
1051
1052 if (numpy.max(abs(xnew0-self.params)) > 0):
1053 if (self.qanytied): self.params = self.tie(self.params, ptied)
1054 x = take(self.params, ifree)
1055
1056
1057
1058 self.status = 2
1059 catch_msg = 'calling MPFIT_FDJAC2'
1060 fjac = self.fdjac2(fcn, x, fvec, step, qulim, ulim, dside,
1061 epsfcn=epsfcn,
1062 autoderivative=autoderivative, dstep=dstep,
1063 functkw=functkw, ifree=ifree, xall=self.params)
1064 if (fjac == None):
1065 self.errmsg = 'WARNING: premature termination by FDJAC2'
1066 return
1067
1068
1069 if (qanylim):
1070 catch_msg = 'zeroing derivatives of pegged parameters'
1071 whlpeg = (nonzero(qllim & (x == llim)))[0]
1072 nlpeg = len(whlpeg)
1073 whupeg = (nonzero(qulim & (x == ulim)))[0]
1074 nupeg = len(whupeg)
1075
1076 if (nlpeg > 0):
1077
1078 for i in range(nlpeg):
1079 sum0 = sum(fvec * fjac[:,whlpeg[i]])
1080 if (sum0 > 0): fjac[:,whlpeg[i]] = 0
1081 if (nupeg > 0):
1082
1083 for i in range(nupeg):
1084 sum0 = sum(fvec * fjac[:,whupeg[i]])
1085 if (sum0 < 0): fjac[:,whupeg[i]] = 0
1086
1087
1088 [fjac, ipvt, wa1, wa2] = self.qrfac(fjac, pivot=1)
1089
1090
1091
1092 catch_msg = 'rescaling diagonal elements'
1093 if (self.niter == 1):
1094 if ((rescale==0) or (len(diag) < n)):
1095 diag = wa2.copy()
1096 wh = (nonzero(diag == 0))[0]
1097 put(diag, wh, 1.)
1098
1099
1100
1101 wa3 = diag * x
1102 xnorm = self.enorm(wa3)
1103 delta = factor*xnorm
1104 if (delta == 0.): delta = factor
1105
1106
1107 catch_msg = 'forming (q transpose)*fvec'
1108 wa4 = fvec.copy()
1109 for j in range(n):
1110 lj = ipvt[j]
1111 temp3 = fjac[j,lj]
1112 if (temp3 != 0):
1113 fj = fjac[j:,lj]
1114 wj = wa4[j:]
1115
1116 wa4[j:] = wj - fj * sum(fj*wj) / temp3
1117 fjac[j,lj] = wa1[j]
1118 qtf[j] = wa4[j]
1119
1120
1121 fjac = fjac[0:n, 0:n]
1122 fjac.shape = [n, n]
1123 temp = fjac.copy()
1124 for i in range(n):
1125 temp[:,i] = fjac[:, ipvt[i]]
1126 fjac = temp.copy()
1127
1128
1129
1130
1131
1132
1133
1134
1135 catch_msg = 'computing the scaled gradient'
1136 gnorm = 0.
1137 if (self.fnorm != 0):
1138 for j in range(n):
1139 l = ipvt[j]
1140 if (wa2[l] != 0):
1141 sum0 = sum(fjac[0:j+1,j]*qtf[0:j+1])/self.fnorm
1142 gnorm = numpy.max([gnorm,abs(sum0/wa2[l])])
1143
1144
1145 if (gnorm <= gtol):
1146 self.status = 4
1147 break
1148 if maxiter == 0:
1149 self.status = 5
1150 break
1151
1152
1153 if (rescale == 0):
1154 diag = choose(diag>wa2, (wa2, diag))
1155
1156
1157 while(1):
1158
1159
1160 catch_msg = 'calculating LM parameter (MPFIT_)'
1161 [fjac, par, wa1, wa2] = self.lmpar(fjac, ipvt, diag, qtf,
1162 delta, wa1, wa2, par=par)
1163
1164 wa1 = -wa1
1165
1166 if (qanylim == 0) and (qminmax == 0):
1167
1168 alpha = 1.
1169 wa2 = x + wa1
1170
1171 else:
1172
1173
1174
1175
1176 alpha = 1.
1177
1178 if (qanylim):
1179
1180 catch_msg = 'checking for a step out of bounds'
1181 if (nlpeg > 0):
1182 put(wa1, whlpeg, clip(
1183 take(wa1, whlpeg), 0., numpy.max(wa1)))
1184 if (nupeg > 0):
1185 put(wa1, whupeg, clip(
1186 take(wa1, whupeg), numpy.min(wa1), 0.))
1187
1188 dwa1 = abs(wa1) > machep
1189 whl = (nonzero(((dwa1!=0.) & qllim) & ((x + wa1) < llim)))[0]
1190 if (len(whl) > 0):
1191 t = ((take(llim, whl) - take(x, whl)) /
1192 take(wa1, whl))
1193 alpha = numpy.min([alpha, numpy.min(t)])
1194 whu = (nonzero(((dwa1!=0.) & qulim) & ((x + wa1) > ulim)))[0]
1195 if (len(whu) > 0):
1196 t = ((take(ulim, whu) - take(x, whu)) /
1197 take(wa1, whu))
1198 alpha = numpy.min([alpha, numpy.min(t)])
1199
1200
1201 if (qminmax):
1202 nwa1 = wa1 * alpha
1203 whmax = (nonzero((qmax != 0.) & (maxstep > 0)))[0]
1204 if (len(whmax) > 0):
1205 mrat = numpy.max(numpy.abs(nwa1[whmax]) /
1206 numpy.abs(maxstep[ifree[whmax]]))
1207 if (mrat > 1): alpha = alpha / mrat
1208
1209
1210 wa1 = wa1 * alpha
1211 wa2 = x + wa1
1212
1213
1214
1215 sgnu = (ulim >= 0) * 2. - 1.
1216 sgnl = (llim >= 0) * 2. - 1.
1217
1218
1219 ulim1 = ulim * (1 - sgnu * machep) - (ulim == 0) * machep
1220 llim1 = llim * (1 + sgnl * machep) + (llim == 0) * machep
1221 wh = (nonzero((qulim!=0) & (wa2 >= ulim1)))[0]
1222 if (len(wh) > 0): put(wa2, wh, take(ulim, wh))
1223 wh = (nonzero((qllim!=0.) & (wa2 <= llim1)))[0]
1224 if (len(wh) > 0): put(wa2, wh, take(llim, wh))
1225
1226 wa3 = diag * wa1
1227 pnorm = self.enorm(wa3)
1228
1229
1230 if (self.niter == 1): delta = numpy.min([delta,pnorm])
1231
1232 put(self.params, ifree, wa2)
1233
1234
1235 mperr = 0
1236 catch_msg = 'calling '+str(fcn)
1237 [self.status, wa4] = self.call(fcn, self.params, functkw)
1238 if (self.status < 0):
1239 self.errmsg = 'WARNING: premature termination by "'+fcn+'"'
1240 return
1241 fnorm1 = self.enorm(wa4)
1242
1243
1244 catch_msg = 'computing convergence criteria'
1245 actred = -1.
1246 if ((0.1 * fnorm1) < self.fnorm): actred = - (fnorm1/self.fnorm)**2 + 1.
1247
1248
1249
1250 for j in range(n):
1251 wa3[j] = 0
1252 wa3[0:j+1] = wa3[0:j+1] + fjac[0:j+1,j]*wa1[ipvt[j]]
1253
1254
1255
1256 temp1 = self.enorm(alpha*wa3)/self.fnorm
1257 temp2 = (sqrt(alpha*par)*pnorm)/self.fnorm
1258 prered = temp1*temp1 + (temp2*temp2)/0.5
1259 dirder = -(temp1*temp1 + temp2*temp2)
1260
1261
1262 ratio = 0.
1263 if (prered != 0): ratio = actred/prered
1264
1265
1266 if (ratio <= 0.25):
1267 if (actred >= 0): temp = .5
1268 else: temp = .5*dirder/(dirder + .5*actred)
1269 if ((0.1*fnorm1) >= self.fnorm) or (temp < 0.1): temp = 0.1
1270 delta = temp*numpy.min([delta,pnorm/0.1])
1271 par = par/temp
1272 else:
1273 if (par == 0) or (ratio >= 0.75):
1274 delta = pnorm/.5
1275 par = .5*par
1276
1277
1278 if (ratio >= 0.0001):
1279
1280 x = wa2
1281 wa2 = diag * x
1282 fvec = wa4
1283 xnorm = self.enorm(wa2)
1284 self.fnorm = fnorm1
1285 self.niter = self.niter + 1
1286
1287
1288 if ((abs(actred) <= ftol) and (prered <= ftol)
1289 and (0.5 * ratio <= 1)): self.status = 1
1290 if delta <= xtol*xnorm: self.status = 2
1291 if ((abs(actred) <= ftol) and (prered <= ftol)
1292 and (0.5 * ratio <= 1) and (self.status == 2)): self.status = 3
1293 if (self.status != 0): break
1294
1295
1296 if (self.niter >= maxiter): self.status = 5
1297 if ((abs(actred) <= machep) and (prered <= machep)
1298 and (0.5*ratio <= 1)): self.status = 6
1299 if delta <= machep*xnorm: self.status = 7
1300 if gnorm <= machep: self.status = 8
1301 if (self.status != 0): break
1302
1303
1304 if (ratio >= 0.0001): break
1305
1306
1307 if ~numpy.all(numpy.isfinite(wa1) & numpy.isfinite(wa2) &
1308 numpy.isfinite(x)) or ~numpy.isfinite(ratio):
1309 errmsg = ('''ERROR: parameter or function value(s) have become
1310 'infinite; check model function for over- 'and underflow''')
1311 self.status = -16
1312 break
1313
1314
1315
1316 if (self.status != 0): break;
1317
1318
1319 catch_msg = 'in the termination phase'
1320
1321 if (len(self.params) == 0):
1322 return
1323
1324
1325 if (nfree == 0): self.params = xall.copy()
1326 else: put(self.params, ifree, x)
1327 if (nprint > 0) and (self.status > 0):
1328 catch_msg = 'calling ' + str(fcn)
1329 [status, fvec] = self.call(fcn, self.params, functkw)
1330 catch_msg = 'in the termination phase'
1331 self.fnorm = self.enorm(fvec)
1332
1333 if ((self.fnorm != None) and (fnorm1 != None)):
1334 self.fnorm = numpy.max([self.fnorm, fnorm1])
1335 self.fnorm = self.fnorm**2.
1336
1337 self.covar = None
1338 self.perror = None
1339
1340 if ((self.status > 0) and (nocovar==0) and (n != None)
1341 and (fjac != None) and (ipvt != None)):
1342 sz = shape(fjac)
1343 if ((n > 0) and (sz[0] >= n) and (sz[1] >= n)
1344 and (len(ipvt) >= n)):
1345
1346 catch_msg = 'computing the covariance matrix'
1347 cv = self.calc_covar(fjac[0:n,0:n], ipvt[0:n])
1348 cv.shape = [n, n]
1349 nn = len(xall)
1350
1351
1352
1353 self.covar = zeros([nn, nn], dtype=float)
1354 for i in range(n):
1355 self.covar[ifree,ifree[i]] = cv[:,i]
1356
1357
1358 catch_msg = 'computing parameter errors'
1359 self.perror = zeros(nn, dtype=float)
1360 d = diagonal(self.covar)
1361 wh = (nonzero(d >= 0))[0]
1362 if len(wh) > 0:
1363 put(self.perror, wh, sqrt(take(d, wh)))
1364 return
1365
1366
1368 return {'params': self.params,
1369 'niter': self.niter,
1370 'params': self.params,
1371 'covar': self.covar,
1372 'perror': self.perror,
1373 'status': self.status,
1374 'debug': self.debug,
1375 'errmsg': self.errmsg,
1376 'fastnorm': self.fastnorm,
1377 'nfev': self.nfev,
1378 'damp': self.damp
1379
1380 }.__str__()
1381
1382
1383
1384 - def defiter(self, fcn, x, iter, fnorm=None, functkw=None,
1385 quiet=0, iterstop=None, parinfo=None,
1386 format=None, pformat='%.10g', dof=1):
1387
1388 if (self.debug): print 'Entering defiter...'
1389 if (quiet): return
1390 if (fnorm == None):
1391 [status, fvec] = self.call(fcn, x, functkw)
1392 fnorm = self.enorm(fvec)**2
1393
1394
1395 nprint = len(x)
1396 print "Iter ", ('%6i' % iter)," CHI-SQUARE = ",('%.10g' % fnorm)," DOF = ", ('%i' % dof)
1397 for i in range(nprint):
1398 if (parinfo != None) and (parinfo[i].has_key('parname')):
1399 p = ' ' + parinfo[i]['parname'] + ' = '
1400 else:
1401 p = ' P' + str(i) + ' = '
1402 if (parinfo != None) and (parinfo[i].has_key('mpprint')):
1403 iprint = parinfo[i]['mpprint']
1404 else:
1405 iprint = 1
1406 if (iprint):
1407 print p + (pformat % x[i]) + ' '
1408 return(0)
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428 - def parinfo(self, parinfo=None, key='a', default=None, n=0):
1449
1450
1451
1452 - def call(self, fcn, x, functkw, fjac=None):
1453 if (self.debug): print 'Entering call...'
1454 if (self.qanytied): x = self.tie(x, self.ptied)
1455 self.nfev = self.nfev + 1
1456 if (fjac == None):
1457 [status, f] = fcn(x, fjac=fjac, **functkw)
1458 if (self.damp > 0):
1459
1460
1461
1462 f = tanh(f/self.damp)
1463 return([status, f])
1464 else:
1465 return(fcn(x, fjac=fjac, **functkw))
1466
1467
1469
1470 if (self.debug): print 'Entering enorm...'
1471
1472
1473
1474
1475
1476
1477
1478
1479 if (self.fastnorm):
1480 ans = sqrt(sum(vec*vec))
1481 else:
1482 agiant = self.machar.rgiant / len(vec)
1483 adwarf = self.machar.rdwarf * len(vec)
1484
1485
1486
1487 mx = numpy.max(vec)
1488 mn = numpy.min(vec)
1489 mx = max([abs(mx), abs(mn)])
1490 if mx == 0: return(vec[0]*0.)
1491 if mx > agiant or mx < adwarf:
1492 ans = mx * sqrt(sum((vec/mx)*(vec/mx)))
1493 else:
1494 ans = sqrt(sum(vec*vec))
1495
1496 return(ans)
1497
1498
1499 - def fdjac2(self, fcn, x, fvec, step=None, ulimited=None, ulimit=None, dside=None,
1500 epsfcn=None, autoderivative=1,
1501 functkw=None, xall=None, ifree=None, dstep=None):
1502
1503 if (self.debug): print 'Entering fdjac2...'
1504 machep = self.machar.machep
1505 if epsfcn == None: epsfcn = machep
1506 if xall == None: xall = x
1507 if ifree == None: ifree = arange(len(xall))
1508 if step == None: step = x * 0.
1509 nall = len(xall)
1510
1511 eps = sqrt(numpy.max([epsfcn, machep]))
1512 m = len(fvec)
1513 n = len(x)
1514
1515
1516 if (autoderivative == 0):
1517 mperr = 0
1518 fjac = zeros(nall, dtype=float)
1519 Put(fjac, ifree, 1.0)
1520 [status, fp] = self.call(fcn, xall, functkw, fjac=fjac)
1521
1522 if len(fjac) != m*nall:
1523 print 'ERROR: Derivative matrix was not computed properly.'
1524 return(None)
1525
1526
1527
1528 fjac.shape = [m,nall]
1529 fjac = -fjac
1530
1531
1532 if len(ifree) < nall:
1533 fjac = fjac[:,ifree]
1534 fjac.shape = [m, n]
1535 return(fjac)
1536
1537 fjac = zeros([m, n], dtype=float)
1538
1539 h = eps * abs(x)
1540
1541
1542
1543 if step != None:
1544 stepi = take(step, ifree)
1545 wh = (nonzero(stepi > 0))[0]
1546 if (len(wh) > 0): put(h, wh, take(stepi, wh))
1547
1548
1549
1550 if (len(dstep) > 0):
1551 dstepi = take(dstep, ifree)
1552 wh = (nonzero(dstepi > 0))[0]
1553 if len(wh) > 0: put(h, wh, abs(take(dstepi,wh)*take(x,wh)))
1554
1555
1556 wh = (nonzero(h == 0))[0]
1557 if len(wh) > 0: put(h, wh, eps)
1558
1559
1560
1561
1562
1563 mask = dside[ifree] == -1
1564 if len(ulimited) > 0 and len(ulimit) > 0:
1565 mask = logical_or((mask!=0), logical_and((ulimited!=0),(x > ulimit-h)))
1566 wh = (nonzero(mask))[0]
1567 if len(wh) > 0: put(h, wh, -take(h, wh))
1568
1569 for j in range(n):
1570 xp = xall.copy()
1571 xp[ifree[j]] = xp[ifree[j]] + h[j]
1572 [status, fp] = self.call(fcn, xp, functkw)
1573 if (status < 0): return(None)
1574
1575 if abs(dside[ifree[j]]) <= 1:
1576
1577
1578 fjac[0:,j] = (fp-fvec)/h[j]
1579
1580 else:
1581
1582 xp[ifree[j]] = xall[ifree[j]] - h[j]
1583
1584 mperr = 0
1585 [status, fm] = self.call(fcn, xp, functkw)
1586 if (status < 0): return(None)
1587
1588
1589 fjac[0:,j] = (fp-fm)/(2*h[j])
1590 return(fjac)
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726 - def qrfac(self, a, pivot=0):
1727
1728 if (self.debug): print 'Entering qrfac...'
1729 machep = self.machar.machep
1730 sz = shape(a)
1731 m = sz[0]
1732 n = sz[1]
1733
1734
1735 acnorm = zeros(n, dtype=float)
1736 for j in range(n):
1737 acnorm[j] = self.enorm(a[:,j])
1738 rdiag = acnorm.copy()
1739 wa = rdiag.copy()
1740 ipvt = arange(n)
1741
1742
1743 minmn = numpy.min([m,n])
1744 for j in range(minmn):
1745 if (pivot != 0):
1746
1747 rmax = numpy.max(rdiag[j:])
1748 kmax = (nonzero(rdiag[j:] == rmax))[0]
1749 ct = len(kmax)
1750 kmax = kmax + j
1751 if ct > 0:
1752 kmax = kmax[0]
1753
1754
1755
1756
1757
1758 if kmax != j:
1759 temp = ipvt[j] ; ipvt[j] = ipvt[kmax] ; ipvt[kmax] = temp
1760 rdiag[kmax] = rdiag[j]
1761 wa[kmax] = wa[j]
1762
1763
1764
1765 lj = ipvt[j]
1766 ajj = a[j:,lj]
1767 ajnorm = self.enorm(ajj)
1768 if ajnorm == 0: break
1769 if a[j,lj] < 0: ajnorm = -ajnorm
1770
1771 ajj = ajj / ajnorm
1772 ajj[0] = ajj[0] + 1
1773
1774 a[j:,lj] = ajj
1775
1776
1777
1778
1779
1780
1781
1782 if (j+1 < n):
1783 for k in range(j+1, n):
1784 lk = ipvt[k]
1785 ajk = a[j:,lk]
1786
1787
1788 if a[j,lj] != 0:
1789 a[j:,lk] = ajk - ajj * sum(ajk*ajj)/a[j,lj]
1790 if ((pivot != 0) and (rdiag[k] != 0)):
1791 temp = a[j,lk]/rdiag[k]
1792 rdiag[k] = rdiag[k] * sqrt(numpy.max([(1.-temp**2), 0.]))
1793 temp = rdiag[k]/wa[k]
1794 if ((0.05*temp*temp) <= machep):
1795 rdiag[k] = self.enorm(a[j+1:,lk])
1796 wa[k] = rdiag[k]
1797 rdiag[j] = -ajnorm
1798 return([a, ipvt, rdiag, acnorm])
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879 - def qrsolv(self, r, ipvt, diag, qtb, sdiag):
1880 if (self.debug): print 'Entering qrsolv...'
1881 sz = shape(r)
1882 m = sz[0]
1883 n = sz[1]
1884
1885
1886
1887
1888 for j in range(n):
1889 r[j:n,j] = r[j,j:n]
1890 x = diagonal(r)
1891 wa = qtb.copy()
1892
1893
1894 for j in range(n):
1895 l = ipvt[j]
1896 if (diag[l] == 0): break
1897 sdiag[j:] = 0
1898 sdiag[j] = diag[l]
1899
1900
1901
1902
1903
1904 qtbpj = 0.
1905 for k in range(j,n):
1906 if (sdiag[k] == 0): break
1907 if (abs(r[k,k]) < abs(sdiag[k])):
1908 cotan = r[k,k]/sdiag[k]
1909 sine = 0.5/sqrt(.25 + .25*cotan*cotan)
1910 cosine = sine*cotan
1911 else:
1912 tang = sdiag[k]/r[k,k]
1913 cosine = 0.5/sqrt(.25 + .25*tang*tang)
1914 sine = cosine*tang
1915
1916
1917
1918 r[k,k] = cosine*r[k,k] + sine*sdiag[k]
1919 temp = cosine*wa[k] + sine*qtbpj
1920 qtbpj = -sine*wa[k] + cosine*qtbpj
1921 wa[k] = temp
1922
1923
1924 if (n > k+1):
1925 temp = cosine*r[k+1:n,k] + sine*sdiag[k+1:n]
1926 sdiag[k+1:n] = -sine*r[k+1:n,k] + cosine*sdiag[k+1:n]
1927 r[k+1:n,k] = temp
1928 sdiag[j] = r[j,j]
1929 r[j,j] = x[j]
1930
1931
1932
1933 nsing = n
1934 wh = (nonzero(sdiag == 0))[0]
1935 if (len(wh) > 0):
1936 nsing = wh[0]
1937 wa[nsing:] = 0
1938
1939 if (nsing >= 1):
1940 wa[nsing-1] = wa[nsing-1]/sdiag[nsing-1]
1941
1942 for j in range(nsing-2,-1,-1):
1943 sum0 = sum(r[j+1:nsing,j]*wa[j+1:nsing])
1944 wa[j] = (wa[j]-sum0)/sdiag[j]
1945
1946
1947 put(x, ipvt, wa)
1948 return(r, x, sdiag)
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047 - def lmpar(self, r, ipvt, diag, qtb, delta, x, sdiag, par=None):
2048
2049 if (self.debug): print 'Entering lmpar...'
2050 dwarf = self.machar.minnum
2051 machep = self.machar.machep
2052 sz = shape(r)
2053 m = sz[0]
2054 n = sz[1]
2055
2056
2057
2058 nsing = n
2059 wa1 = qtb.copy()
2060 rthresh = numpy.max(numpy.abs(diagonal(r))) * machep
2061 wh = (nonzero(numpy.abs(diagonal(r)) < rthresh))[0]
2062 if len(wh) > 0:
2063 nsing = wh[0]
2064 wa1[wh[0]:] = 0
2065 if nsing >= 1:
2066
2067 for j in range(nsing-1,-1,-1):
2068 wa1[j] = wa1[j]/r[j,j]
2069 if (j-1 >= 0):
2070 wa1[0:j] = wa1[0:j] - r[0:j,j]*wa1[j]
2071
2072
2073 put(x, ipvt, wa1)
2074
2075
2076
2077 iter = 0
2078 wa2 = diag * x
2079 dxnorm = self.enorm(wa2)
2080 fp = dxnorm - delta
2081 if (fp <= 0.1*delta):
2082 return[r, 0., x, sdiag]
2083
2084
2085
2086
2087
2088 parl = 0.
2089 if nsing >= n:
2090 wa1 = take(diag, ipvt)*take(wa2, ipvt)/dxnorm
2091 wa1[0] = wa1[0] / r[0,0]
2092 for j in range(1,n):
2093 sum0 = sum(r[0:j,j]*wa1[0:j])
2094 wa1[j] = (wa1[j] - sum0)/r[j,j]
2095
2096 temp = self.enorm(wa1)
2097 parl = ((fp/delta)/temp)/temp
2098
2099
2100 for j in range(n):
2101 sum0 = sum(r[0:j+1,j]*qtb[0:j+1])
2102 wa1[j] = sum0/diag[ipvt[j]]
2103 gnorm = self.enorm(wa1)
2104 paru = gnorm/delta
2105 if paru == 0: paru = dwarf/numpy.min([delta,0.1])
2106
2107
2108
2109
2110 par = numpy.max([par,parl])
2111 par = numpy.min([par,paru])
2112 if par == 0: par = gnorm/dxnorm
2113
2114
2115 while(1):
2116 iter = iter + 1
2117
2118
2119 if par == 0: par = numpy.max([dwarf, paru*0.001])
2120 temp = sqrt(par)
2121 wa1 = temp * diag
2122 [r, x, sdiag] = self.qrsolv(r, ipvt, wa1, qtb, sdiag)
2123 wa2 = diag*x
2124 dxnorm = self.enorm(wa2)
2125 temp = fp
2126 fp = dxnorm - delta
2127
2128 if ((abs(fp) <= 0.1*delta) or
2129 ((parl == 0) and (fp <= temp) and (temp < 0)) or
2130 (iter == 10)): break;
2131
2132
2133 wa1 = take(diag, ipvt)*take(wa2, ipvt)/dxnorm
2134
2135 for j in range(n-1):
2136 wa1[j] = wa1[j]/sdiag[j]
2137 wa1[j+1:n] = wa1[j+1:n] - r[j+1:n,j]*wa1[j]
2138 wa1[n-1] = wa1[n-1]/sdiag[n-1]
2139
2140 temp = self.enorm(wa1)
2141 parc = ((fp/delta)/temp)/temp
2142
2143
2144 if fp > 0: parl = numpy.max([parl,par])
2145 if fp < 0: paru = numpy.min([paru,par])
2146
2147
2148 par = numpy.max([parl, par+parc])
2149
2150
2151
2152
2153 return[r, par, x, sdiag]
2154
2155
2156
2157 - def tie(self, p, ptied=None):
2158 if (self.debug): print 'Entering tie...'
2159 if (ptied == None): return
2160 for i in range(len(ptied)):
2161 if ptied[i] == '': continue
2162 cmd = 'p[' + str(i) + '] = ' + ptied[i]
2163 exec(cmd)
2164 return(p)
2165
2166
2167
2168
2169
2170
2171
2172
2173
2174
2175
2176
2177
2178
2179
2180
2181
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191
2192
2193
2194
2195
2196
2197
2198
2199
2200
2201
2202
2203
2204
2205
2206
2207
2208
2209
2210
2211
2212
2213
2214
2215
2216
2217
2218
2219
2220
2221
2222
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234 - def calc_covar(self, rr, ipvt=None, tol=1.e-14):
2235
2236 if (self.debug): print 'Entering calc_covar...'
2237 if rank(rr) != 2:
2238 print 'ERROR: r must be a two-dimensional matrix'
2239 return(-1)
2240 s = shape(rr)
2241 n = s[0]
2242 if s[0] != s[1]:
2243 print 'ERROR: r must be a square matrix'
2244 return(-1)
2245
2246 if (ipvt == None): ipvt = arange(n)
2247 r = rr.copy()
2248 r.shape = [n,n]
2249
2250
2251 l = -1
2252 tolr = tol * abs(r[0,0])
2253 for k in range(n):
2254 if (abs(r[k,k]) <= tolr): break
2255 r[k,k] = 1./r[k,k]
2256 for j in range(k):
2257 temp = r[k,k] * r[j,k]
2258 r[j,k] = 0.
2259 r[0:j+1,k] = r[0:j+1,k] - temp*r[0:j+1,j]
2260 l = k
2261
2262
2263
2264 if l >= 0:
2265 for k in range(l+1):
2266 for j in range(k):
2267 temp = r[j,k]
2268 r[0:j+1,j] = r[0:j+1,j] + temp*r[0:j+1,k]
2269 temp = r[k,k]
2270 r[0:k+1,k] = temp * r[0:k+1,k]
2271
2272
2273
2274 wa = repeat([r[0,0]], n)
2275 for j in range(n):
2276 jj = ipvt[j]
2277 sing = j > l
2278 for i in range(j+1):
2279 if sing: r[i,j] = 0.
2280 ii = ipvt[i]
2281 if ii > jj: r[ii,jj] = r[i,j]
2282 if ii < jj: r[jj,ii] = r[i,j]
2283 wa[jj] = r[j,j]
2284
2285
2286 for j in range(n):
2287 r[0:j+1,j] = r[j,0:j+1]
2288 r[j,j] = wa[j]
2289
2290 return(r)
2291
2294 if (double == 0):
2295 self.machep = 1.19209e-007
2296 self.maxnum = 3.40282e+038
2297 self.minnum = 1.17549e-038
2298 self.maxgam = 171.624376956302725
2299 else:
2300 self.machep = 2.2204460e-016
2301 self.maxnum = 1.7976931e+308
2302 self.minnum = 2.2250739e-308
2303 self.maxgam = 171.624376956302725
2304
2305 self.maxlog = log(self.maxnum)
2306 self.minlog = log(self.minnum)
2307 self.rdwarf = sqrt(self.minnum*1.5) * 10
2308 self.rgiant = sqrt(self.maxnum) * 0.1
2309