Package omero :: Package util :: Module mpfit
[hide private]
[frames] | no frames]

Source Code for Module omero.util.mpfit

   1  #!/usr/bin/env python 
   2  # -*- coding: utf-8 -*- 
   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  #        Original FORTRAN documentation 
 417  #        ********** 
 418  # 
 419  #        subroutine lmdif 
 420  # 
 421  #        the purpose of lmdif is to minimize the sum of the squares of 
 422  #        m nonlinear functions in n variables by a modification of 
 423  #        the levenberg-marquardt algorithm. the user must provide a 
 424  #        subroutine which calculates the functions. the jacobian is 
 425  #        then calculated by a forward-difference approximation. 
 426  # 
 427  #        the subroutine statement is 
 428  # 
 429  #          subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn, 
 430  #                                               diag,mode,factor,nprint,info,nfev,fjac, 
 431  #                                               ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4) 
 432  # 
 433  #        where 
 434  # 
 435  #          fcn is the name of the user-supplied subroutine which 
 436  #                calculates the functions. fcn must be declared 
 437  #                in an external statement in the user calling 
 438  #                program, and should be written as follows. 
 439  # 
 440  #                subroutine fcn(m,n,x,fvec,iflag) 
 441  #                integer m,n,iflag 
 442  #                double precision x(n),fvec(m) 
 443  #                ---------- 
 444  #                calculate the functions at x and 
 445  #                return this vector in fvec. 
 446  #                ---------- 
 447  #                return 
 448  #                end 
 449  # 
 450  #                the value of iflag should not be changed by fcn unless 
 451  #                the user wants to terminate execution of lmdif. 
 452  #                in this case set iflag to a negative integer. 
 453  # 
 454  #          m is a positive integer input variable set to the number 
 455  #                of functions. 
 456  # 
 457  #          n is a positive integer input variable set to the number 
 458  #                of variables. n must not exceed m. 
 459  # 
 460  #          x is an array of length n. on input x must contain 
 461  #                an initial estimate of the solution vector. on output x 
 462  #                contains the final estimate of the solution vector. 
 463  # 
 464  #          fvec is an output array of length m which contains 
 465  #                the functions evaluated at the output x. 
 466  # 
 467  #          ftol is a nonnegative input variable. termination 
 468  #                occurs when both the actual and predicted relative 
 469  #                reductions in the sum of squares are at most ftol. 
 470  #                therefore, ftol measures the relative error desired 
 471  #                in the sum of squares. 
 472  # 
 473  #          xtol is a nonnegative input variable. termination 
 474  #                occurs when the relative error between two consecutive 
 475  #                iterates is at most xtol. therefore, xtol measures the 
 476  #                relative error desired in the approximate solution. 
 477  # 
 478  #          gtol is a nonnegative input variable. termination 
 479  #                occurs when the cosine of the angle between fvec and 
 480  #                any column of the jacobian is at most gtol in absolute 
 481  #                value. therefore, gtol measures the orthogonality 
 482  #                desired between the function vector and the columns 
 483  #                of the jacobian. 
 484  # 
 485  #          maxfev is a positive integer input variable. termination 
 486  #                occurs when the number of calls to fcn is at least 
 487  #                maxfev by the end of an iteration. 
 488  # 
 489  #          epsfcn is an input variable used in determining a suitable 
 490  #                step length for the forward-difference approximation. this 
 491  #                approximation assumes that the relative errors in the 
 492  #                functions are of the order of epsfcn. if epsfcn is less 
 493  #                than the machine precision, it is assumed that the relative 
 494  #                errors in the functions are of the order of the machine 
 495  #                precision. 
 496  # 
 497  #          diag is an array of length n. if mode = 1 (see 
 498  #                below), diag is internally set. if mode = 2, diag 
 499  #                must contain positive entries that serve as 
 500  #                multiplicative scale factors for the variables. 
 501  # 
 502  #          mode is an integer input variable. if mode = 1, the 
 503  #                variables will be scaled internally. if mode = 2, 
 504  #                the scaling is specified by the input diag. other 
 505  #                values of mode are equivalent to mode = 1. 
 506  # 
 507  #          factor is a positive input variable used in determining the 
 508  #                initial step bound. this bound is set to the product of 
 509  #                factor and the euclidean norm of diag*x if nonzero, or else 
 510  #                to factor itself. in most cases factor should lie in the 
 511  #                interval (.1,100.). 100. is a generally recommended value. 
 512  # 
 513  #          nprint is an integer input variable that enables controlled 
 514  #                printing of iterates if it is positive. in this case, 
 515  #                fcn is called with iflag = 0 at the beginning of the first 
 516  #                iteration and every nprint iterations thereafter and 
 517  #                immediately prior to return, with x and fvec available 
 518  #                for printing. if nprint is not positive, no special calls 
 519  #                of fcn with iflag = 0 are made. 
 520  # 
 521  #          info is an integer output variable. if the user has 
 522  #                terminated execution, info is set to the (negative) 
 523  #                value of iflag. see description of fcn. otherwise, 
 524  #                info is set as follows. 
 525  # 
 526  #                info = 0  improper input parameters. 
 527  # 
 528  #                info = 1  both actual and predicted relative reductions 
 529  #                                  in the sum of squares are at most ftol. 
 530  # 
 531  #                info = 2  relative error between two consecutive iterates 
 532  #                                  is at most xtol. 
 533  # 
 534  #                info = 3  conditions for info = 1 and info = 2 both hold. 
 535  # 
 536  #                info = 4  the cosine of the angle between fvec and any 
 537  #                                  column of the jacobian is at most gtol in 
 538  #                                  absolute value. 
 539  # 
 540  #                info = 5  number of calls to fcn has reached or 
 541  #                                  exceeded maxfev. 
 542  # 
 543  #                info = 6  ftol is too small. no further reduction in 
 544  #                                  the sum of squares is possible. 
 545  # 
 546  #                info = 7  xtol is too small. no further improvement in 
 547  #                                  the approximate solution x is possible. 
 548  # 
 549  #                info = 8  gtol is too small. fvec is orthogonal to the 
 550  #                                  columns of the jacobian to machine precision. 
 551  # 
 552  #          nfev is an integer output variable set to the number of 
 553  #                calls to fcn. 
 554  # 
 555  #          fjac is an output m by n array. the upper n by n submatrix 
 556  #                of fjac contains an upper triangular matrix r with 
 557  #                diagonal elements of nonincreasing magnitude such that 
 558  # 
 559  #                               t        t                 t 
 560  #                          p *(jac *jac)*p = r *r, 
 561  # 
 562  #                where p is a permutation matrix and jac is the final 
 563  #                calculated jacobian. column j of p is column ipvt(j) 
 564  #                (see below) of the identity matrix. the lower trapezoidal 
 565  #                part of fjac contains information generated during 
 566  #                the computation of r. 
 567  # 
 568  #          ldfjac is a positive integer input variable not less than m 
 569  #                which specifies the leading dimension of the array fjac. 
 570  # 
 571  #          ipvt is an integer output array of length n. ipvt 
 572  #                defines a permutation matrix p such that jac*p = q*r, 
 573  #                where jac is the final calculated jacobian, q is 
 574  #                orthogonal (not stored), and r is upper triangular 
 575  #                with diagonal elements of nonincreasing magnitude. 
 576  #                column j of p is column ipvt(j) of the identity matrix. 
 577  # 
 578  #          qtf is an output array of length n which contains 
 579  #                the first n elements of the vector (q transpose)*fvec. 
 580  # 
 581  #          wa1, wa2, and wa3 are work arrays of length n. 
 582  # 
 583  #          wa4 is a work array of length m. 
 584  # 
 585  #        subprograms called 
 586  # 
 587  #          user-supplied ...... fcn 
 588  # 
 589  #          minpack-supplied ... dpmpar,enorm,fdjac2,,qrfac 
 590  # 
 591  #          fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod 
 592  # 
 593  #        argonne national laboratory. minpack project. march 1980. 
 594  #        burton s. garbow, kenneth e. hillstrom, jorge j. more 
 595  # 
 596  #        ********** 
 597   
598 -class mpfit:
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 # See #4749 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 # Invalid input flag set while we check inputs 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 ## Parameter damping doesn't work when user is providing their own 877 ## gradients. 878 if (self.damp != 0) and (autoderivative == 0): 879 self.errmsg = 'ERROR: keywords DAMP and AUTODERIVATIVE are mutually exclusive' 880 return 881 882 ## Parameters can either be stored in parinfo, or x. x takes precedence if it exists 883 if (xall == None) and (parinfo == None): 884 self.errmsg = 'ERROR: must pass parameters in P or PARINFO' 885 return 886 887 ## Be sure that PARINFO is of the right type 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 ## If the parameters were not specified at the command line, then 901 ## extract them from PARINFO 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 ## Make sure parameters are Numeric arrays of type Float 909 xall = asarray(xall, float) 910 911 npar = len(xall) 912 self.fnorm = -1. 913 fnorm1 = -1. 914 915 ## TIED parameters? 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 ## FIXED parameters ? 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] != '') ## Tied parameters are also effectively fixed 928 929 ## Finite differencing step, absolute and relative, and sidedness of deriv. 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 ## Maximum and minimum steps allowed to be taken in one iteration 935 maxstep = self.parinfo(parinfo, 'mpmaxstep', default=0., n=npar) 936 minstep = self.parinfo(parinfo, 'mpminstep', default=0., n=npar) 937 qmin = minstep * 0 ## Remove minstep for now!! 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 ## Finish up the free parameters 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 ## Compose only VARYING parameters 954 self.params = xall ## self.params is the set of parameters to be returned 955 x = take(self.params, ifree) ## x is the set of free parameters 956 957 ## LIMITED parameters ? 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 ## Error checking on limits in parinfo 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 ## Transfer structure values to local variables 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 ## Fill in local variables with dummy values 985 qulim = zeros(nfree) 986 ulim = x * 0. 987 qllim = qulim 988 llim = x * 0. 989 qanylim = 0 990 991 n = len(x) 992 ## Check input parameters for errors 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 # Make sure x is a Numeric array of type Float 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 ## Initialize Levelberg-Marquardt parameter and iteration counter 1021 1022 par = 0. 1023 self.niter = 1 1024 qtf = x * 0. 1025 self.status = 0 1026 1027 ## Beginning of the outer loop 1028 1029 while(1): 1030 1031 ## If requested, call fcn to enable printing of iterates 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 ## Check for user termination 1047 if (self.status < 0): 1048 self.errmsg = 'WARNING: premature termination by ' + str(iterfunct) 1049 return 1050 1051 ## If parameters were changed (grrr..) then re-tie 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 ## Calculate the jacobian matrix 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 ## Determine if any of the parameters are pegged at the limits 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 ## See if any "pegged" values should keep their derivatives 1076 if (nlpeg > 0): 1077 ## Total derivative of sum wrt lower pegged parameters 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 ## Total derivative of sum wrt upper pegged parameters 1083 for i in range(nupeg): 1084 sum0 = sum(fvec * fjac[:,whupeg[i]]) 1085 if (sum0 < 0): fjac[:,whupeg[i]] = 0 1086 1087 ## Compute the QR factorization of the jacobian 1088 [fjac, ipvt, wa1, wa2] = self.qrfac(fjac, pivot=1) 1089 1090 ## On the first iteration if "diag" is unspecified, scale 1091 ## according to the norms of the columns of the initial jacobian 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 ## On the first iteration, calculate the norm of the scaled x 1100 ## and initialize the step bound delta 1101 wa3 = diag * x 1102 xnorm = self.enorm(wa3) 1103 delta = factor*xnorm 1104 if (delta == 0.): delta = factor 1105 1106 ## Form (q transpose)*fvec and store the first n components in qtf 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 ## *** optimization wa4(j:*) 1116 wa4[j:] = wj - fj * sum(fj*wj) / temp3 1117 fjac[j,lj] = wa1[j] 1118 qtf[j] = wa4[j] 1119 ## From this point on, only the square matrix, consisting of the 1120 ## triangle of R, is needed. 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 ## Check for overflow. This should be a cheap test here since FJAC 1129 ## has been reduced to a (small) square matrix, and the test is 1130 ## O(N^2). 1131 #wh = where(finite(fjac) EQ 0, ct) 1132 #if ct GT 0 then goto, FAIL_OVERFLOW 1133 1134 ## Compute the norm of the scaled gradient 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 ## Test for convergence of the gradient norm 1145 if (gnorm <= gtol): 1146 self.status = 4 1147 break 1148 if maxiter == 0: 1149 self.status = 5 1150 break 1151 1152 ## Rescale if necessary 1153 if (rescale == 0): 1154 diag = choose(diag>wa2, (wa2, diag)) 1155 1156 ## Beginning of the inner loop 1157 while(1): 1158 1159 ## Determine the levenberg-marquardt parameter 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 ## Store the direction p and x+p. Calculate the norm of p 1164 wa1 = -wa1 1165 1166 if (qanylim == 0) and (qminmax == 0): 1167 ## No parameter limits, so just move to new position WA2 1168 alpha = 1. 1169 wa2 = x + wa1 1170 1171 else: 1172 1173 ## Respect the limits. If a step were to go out of bounds, then 1174 ## we should take a step in the same direction but shorter distance. 1175 ## The step should take us right to the limit in that case. 1176 alpha = 1. 1177 1178 if (qanylim): 1179 ## Do not allow any steps out of bounds 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 ## Obey any max step values. 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 ## Scale the resulting vector 1210 wa1 = wa1 * alpha 1211 wa2 = x + wa1 1212 1213 ## Adjust the final output values. If the step put us exactly 1214 ## on a boundary, make sure it is exact. 1215 sgnu = (ulim >= 0) * 2. - 1. 1216 sgnl = (llim >= 0) * 2. - 1. 1217 ## Handles case of 1218 ## ... nonzero *LIM ... ...zero * LIM 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 # endelse 1226 wa3 = diag * wa1 1227 pnorm = self.enorm(wa3) 1228 1229 ## On the first iteration, adjust the initial step bound 1230 if (self.niter == 1): delta = numpy.min([delta,pnorm]) 1231 1232 put(self.params, ifree, wa2) 1233 1234 ## Evaluate the function at x+p and calculate its norm 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 ## Compute the scaled actual reduction 1244 catch_msg = 'computing convergence criteria' 1245 actred = -1. 1246 if ((0.1 * fnorm1) < self.fnorm): actred = - (fnorm1/self.fnorm)**2 + 1. 1247 1248 ## Compute the scaled predicted reduction and the scaled directional 1249 ## derivative 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 ## Remember, alpha is the fraction of the full LM step actually 1255 ## taken 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 ## Compute the ratio of the actual to the predicted reduction. 1262 ratio = 0. 1263 if (prered != 0): ratio = actred/prered 1264 1265 ## Update the step bound 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 ## Test for successful iteration 1278 if (ratio >= 0.0001): 1279 ## Successful iteration. Update x, fvec, and their norms 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 ## Tests for convergence 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 ## Tests for termination and stringent tolerances 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 ## End of inner loop. Repeat if iteration unsuccessful 1304 if (ratio >= 0.0001): break 1305 1306 ## Check for over/underflow 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 ##wh = where(finite(wa1) EQ 0 OR finite(wa2) EQ 0 OR finite(x) EQ 0, ct) 1314 ##if ct GT 0 OR finite(ratio) EQ 0 then begin 1315 1316 if (self.status != 0): break; 1317 ## End of outer loop. 1318 1319 catch_msg = 'in the termination phase' 1320 ## Termination, either normal or user imposed. 1321 if (len(self.params) == 0): 1322 return 1323 #if xnew==None: # failed patch WR 1324 # return 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 ## (very carefully) set the covariance matrix COVAR 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 ## Fill in actual covariance matrix, accounting for fixed 1352 ## parameters. 1353 self.covar = zeros([nn, nn], dtype=float) 1354 for i in range(n): 1355 self.covar[ifree,ifree[i]] = cv[:,i] 1356 1357 ## Compute errors in parameters 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
1367 - def __str__(self):
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 #,'machar':self.machar 1380 }.__str__()
1381 1382 ## Default procedure to be called every iteration. It simply prints 1383 ## the parameter values.
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 ## Determine which parameters to print 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 ## DO_ITERSTOP: 1411 ## if keyword_set(iterstop) then begin 1412 ## k = get_kbrd(0) 1413 ## if k EQ string(byte(7)) then begin 1414 ## message, 'WARNING: minimization not complete', /info 1415 ## print, 'Do you want to terminate this procedure? (y/n)', $ 1416 ## format='(A,$)' 1417 ## k = '' 1418 ## read, k 1419 ## if strupcase(strmid(k,0,1)) EQ 'Y' then begin 1420 ## message, 'WARNING: Procedure is terminating.', /info 1421 ## mperr = -1 1422 ## endif 1423 ## endif 1424 ## endif 1425 1426 1427 ## Procedure to parse the parameter values in PARINFO, which is a list of dictionaries
1428 - def parinfo(self, parinfo=None, key='a', default=None, n=0):
1429 if (self.debug): print 'Entering parinfo...' 1430 if (n == 0) and (parinfo != None): n = len(parinfo) 1431 if (n == 0): 1432 values = default 1433 return(values) 1434 values = [] 1435 for i in range(n): 1436 if ((parinfo != None) and (parinfo[i].has_key(key))): 1437 values.append(parinfo[i][key]) 1438 else: 1439 values.append(default) 1440 1441 # Convert to numeric arrays if possible 1442 test = default 1443 if (type(default) == types.ListType): test=default[0] 1444 if (isinstance(test, types.IntType)): 1445 values = asarray(values, int) 1446 elif (isinstance(test, types.FloatType)): 1447 values = asarray(values, float) 1448 return(values)
1449 1450 ## Call user function or procedure, with _EXTRA or not, with 1451 ## derivatives or not.
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 ## Apply the damping if requested. This replaces the residuals 1460 ## with their hyperbolic tangent. Thus residuals larger than 1461 ## DAMP are essentially clipped. 1462 f = tanh(f/self.damp) 1463 return([status, f]) 1464 else: 1465 return(fcn(x, fjac=fjac, **functkw))
1466 1467
1468 - def enorm(self, vec):
1469 1470 if (self.debug): print 'Entering enorm...' 1471 ## NOTE: it turns out that, for systems that have a lot of data 1472 ## points, this routine is a big computing bottleneck. The extended 1473 ## computations that need to be done cannot be effectively 1474 ## vectorized. The introduction of the FASTNORM configuration 1475 ## parameter allows the user to select a faster routine, which is 1476 ## based on TOTAL() alone. 1477 1478 # Very simple-minded sum-of-squares 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 ## This is hopefully a compromise between speed and robustness. 1486 ## Need to do this because of the possibility of over- or underflow. 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 ## Compute analytical derivative if requested 1516 if (autoderivative == 0): 1517 mperr = 0 1518 fjac = zeros(nall, dtype=float) 1519 Put(fjac, ifree, 1.0) ## Specify which parameters need derivatives 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 ## This definition is consistent with CURVEFIT 1527 ## Sign error found (thanks Jesus Fernandez <fernande@irm.chu-caen.fr>) 1528 fjac.shape = [m,nall] 1529 fjac = -fjac 1530 1531 ## Select only the free parameters 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 ## if STEP is given, use that 1542 ## STEP includes the fixed parameters 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 ## if relative step is given, use that 1549 ## DSTEP includes the fixed parameters 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 ## In case any of the step values are zero 1556 wh = (nonzero(h == 0))[0] 1557 if len(wh) > 0: put(h, wh, eps) 1558 1559 ## Reverse the sign of the step if we are up against the parameter 1560 ## limit, or if the user requested it. 1561 ## DSIDE includes the fixed parameters (ULIMITED/ULIMIT have only 1562 ## varying ones) 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 ## Loop through parameters, computing the derivative for each 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 ## COMPUTE THE ONE-SIDED DERIVATIVE 1577 ## Note optimization fjac(0:*,j) 1578 fjac[0:,j] = (fp-fvec)/h[j] 1579 1580 else: 1581 ## COMPUTE THE TWO-SIDED DERIVATIVE 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 ## Note optimization fjac(0:*,j) 1589 fjac[0:,j] = (fp-fm)/(2*h[j]) 1590 return(fjac)
1591 1592 1593 1594 # Original FORTRAN documentation 1595 # ********** 1596 # 1597 # subroutine qrfac 1598 # 1599 # this subroutine uses householder transformations with column 1600 # pivoting (optional) to compute a qr factorization of the 1601 # m by n matrix a. that is, qrfac determines an orthogonal 1602 # matrix q, a permutation matrix p, and an upper trapezoidal 1603 # matrix r with diagonal elements of nonincreasing magnitude, 1604 # such that a*p = q*r. the householder transformation for 1605 # column k, k = 1,2,...,min(m,n), is of the form 1606 # 1607 # t 1608 # i - (1/u(k))*u*u 1609 # 1610 # where u has zeros in the first k-1 positions. the form of 1611 # this transformation and the method of pivoting first 1612 # appeared in the corresponding linpack subroutine. 1613 # 1614 # the subroutine statement is 1615 # 1616 # subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa) 1617 # 1618 # where 1619 # 1620 # m is a positive integer input variable set to the number 1621 # of rows of a. 1622 # 1623 # n is a positive integer input variable set to the number 1624 # of columns of a. 1625 # 1626 # a is an m by n array. on input a contains the matrix for 1627 # which the qr factorization is to be computed. on output 1628 # the strict upper trapezoidal part of a contains the strict 1629 # upper trapezoidal part of r, and the lower trapezoidal 1630 # part of a contains a factored form of q (the non-trivial 1631 # elements of the u vectors described above). 1632 # 1633 # lda is a positive integer input variable not less than m 1634 # which specifies the leading dimension of the array a. 1635 # 1636 # pivot is a logical input variable. if pivot is set true, 1637 # then column pivoting is enforced. if pivot is set false, 1638 # then no column pivoting is done. 1639 # 1640 # ipvt is an integer output array of length lipvt. ipvt 1641 # defines the permutation matrix p such that a*p = q*r. 1642 # column j of p is column ipvt(j) of the identity matrix. 1643 # if pivot is false, ipvt is not referenced. 1644 # 1645 # lipvt is a positive integer input variable. if pivot is false, 1646 # then lipvt may be as small as 1. if pivot is true, then 1647 # lipvt must be at least n. 1648 # 1649 # rdiag is an output array of length n which contains the 1650 # diagonal elements of r. 1651 # 1652 # acnorm is an output array of length n which contains the 1653 # norms of the corresponding columns of the input matrix a. 1654 # if this information is not needed, then acnorm can coincide 1655 # with rdiag. 1656 # 1657 # wa is a work array of length n. if pivot is false, then wa 1658 # can coincide with rdiag. 1659 # 1660 # subprograms called 1661 # 1662 # minpack-supplied ... dpmpar,enorm 1663 # 1664 # fortran-supplied ... dmax1,dsqrt,min0 1665 # 1666 # argonne national laboratory. minpack project. march 1980. 1667 # burton s. garbow, kenneth e. hillstrom, jorge j. more 1668 # 1669 # ********** 1670 # 1671 # PIVOTING / PERMUTING: 1672 # 1673 # Upon return, A(*,*) is in standard parameter order, A(*,IPVT) is in 1674 # permuted order. 1675 # 1676 # RDIAG is in permuted order. 1677 # ACNORM is in standard parameter order. 1678 # 1679 # 1680 # NOTE: in IDL the factors appear slightly differently than described 1681 # above. The matrix A is still m x n where m >= n. 1682 # 1683 # The "upper" triangular matrix R is actually stored in the strict 1684 # lower left triangle of A under the standard notation of IDL. 1685 # 1686 # The reflectors that generate Q are in the upper trapezoid of A upon 1687 # output. 1688 # 1689 # EXAMPLE: decompose the matrix [[9.,2.,6.],[4.,8.,7.]] 1690 # aa = [[9.,2.,6.],[4.,8.,7.]] 1691 # mpfit_qrfac, aa, aapvt, rdiag, aanorm 1692 # IDL> print, aa 1693 # 1.81818* 0.181818* 0.545455* 1694 # -8.54545+ 1.90160* 0.432573* 1695 # IDL> print, rdiag 1696 # -11.0000+ -7.48166+ 1697 # 1698 # The components marked with a * are the components of the 1699 # reflectors, and those marked with a + are components of R. 1700 # 1701 # To reconstruct Q and R we proceed as follows. First R. 1702 # r = fltarr(m, n) 1703 # for i = 0, n-1 do r(0:i,i) = aa(0:i,i) # fill in lower diag 1704 # r(lindgen(n)*(m+1)) = rdiag 1705 # 1706 # Next, Q, which are composed from the reflectors. Each reflector v 1707 # is taken from the upper trapezoid of aa, and converted to a matrix 1708 # via (I - 2 vT . v / (v . vT)). 1709 # 1710 # hh = ident ## identity matrix 1711 # for i = 0, n-1 do begin 1712 # v = aa(*,i) & if i GT 0 then v(0:i-1) = 0 ## extract reflector 1713 # hh = hh ## (ident - 2*(v # v)/total(v * v)) ## generate matrix 1714 # endfor 1715 # 1716 # Test the result: 1717 # IDL> print, hh ## transpose(r) 1718 # 9.00000 4.00000 1719 # 2.00000 8.00000 1720 # 6.00000 7.00000 1721 # 1722 # Note that it is usually never necessary to form the Q matrix 1723 # explicitly, and MPFIT does not. 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 ## Compute the initial column norms and initialize arrays 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 ## Reduce a to r with householder transformations 1743 minmn = numpy.min([m,n]) 1744 for j in range(minmn): 1745 if (pivot != 0): 1746 ## Bring the column of largest norm into the pivot position 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 ## Exchange rows via the pivot only. Avoid actually exchanging 1755 ## the rows, in case there is lots of memory transfer. The 1756 ## exchange occurs later, within the body of MPFIT, after the 1757 ## extraneous columns of the matrix have been shed. 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 ## Compute the householder transformation to reduce the jth 1764 ## column of A to a multiple of the jth unit vector 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 ## *** Note optimization a(j:*,j) 1774 a[j:,lj] = ajj 1775 1776 ## Apply the transformation to the remaining columns 1777 ## and update the norms 1778 1779 ## NOTE to SELF: tried to optimize this by removing the loop, 1780 ## but it actually got slower. Reverted to "for" loop to keep 1781 ## it simple. 1782 if (j+1 < n): 1783 for k in range(j+1, n): 1784 lk = ipvt[k] 1785 ajk = a[j:,lk] 1786 ## *** Note optimization a(j:*,lk) 1787 ## (corrected 20 Jul 2000) 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 # Original FORTRAN documentation 1802 # ********** 1803 # 1804 # subroutine qrsolv 1805 # 1806 # given an m by n matrix a, an n by n diagonal matrix d, 1807 # and an m-vector b, the problem is to determine an x which 1808 # solves the system 1809 # 1810 # a*x = b , d*x = 0 , 1811 # 1812 # in the least squares sense. 1813 # 1814 # this subroutine completes the solution of the problem 1815 # if it is provided with the necessary information from the 1816 # factorization, with column pivoting, of a. that is, if 1817 # a*p = q*r, where p is a permutation matrix, q has orthogonal 1818 # columns, and r is an upper triangular matrix with diagonal 1819 # elements of nonincreasing magnitude, then qrsolv expects 1820 # the full upper triangle of r, the permutation matrix p, 1821 # and the first n components of (q transpose)*b. the system 1822 # a*x = b, d*x = 0, is then equivalent to 1823 # 1824 # t t 1825 # r*z = q *b , p *d*p*z = 0 , 1826 # 1827 # where x = p*z. if this system does not have full rank, 1828 # then a least squares solution is obtained. on output qrsolv 1829 # also provides an upper triangular matrix s such that 1830 # 1831 # t t t 1832 # p *(a *a + d*d)*p = s *s . 1833 # 1834 # s is computed within qrsolv and may be of separate interest. 1835 # 1836 # the subroutine statement is 1837 # 1838 # subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa) 1839 # 1840 # where 1841 # 1842 # n is a positive integer input variable set to the order of r. 1843 # 1844 # r is an n by n array. on input the full upper triangle 1845 # must contain the full upper triangle of the matrix r. 1846 # on output the full upper triangle is unaltered, and the 1847 # strict lower triangle contains the strict upper triangle 1848 # (transposed) of the upper triangular matrix s. 1849 # 1850 # ldr is a positive integer input variable not less than n 1851 # which specifies the leading dimension of the array r. 1852 # 1853 # ipvt is an integer input array of length n which defines the 1854 # permutation matrix p such that a*p = q*r. column j of p 1855 # is column ipvt(j) of the identity matrix. 1856 # 1857 # diag is an input array of length n which must contain the 1858 # diagonal elements of the matrix d. 1859 # 1860 # qtb is an input array of length n which must contain the first 1861 # n elements of the vector (q transpose)*b. 1862 # 1863 # x is an output array of length n which contains the least 1864 # squares solution of the system a*x = b, d*x = 0. 1865 # 1866 # sdiag is an output array of length n which contains the 1867 # diagonal elements of the upper triangular matrix s. 1868 # 1869 # wa is a work array of length n. 1870 # 1871 # subprograms called 1872 # 1873 # fortran-supplied ... dabs,dsqrt 1874 # 1875 # argonne national laboratory. minpack project. march 1980. 1876 # burton s. garbow, kenneth e. hillstrom, jorge j. more 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 ## copy r and (q transpose)*b to preserve input and initialize s. 1886 ## in particular, save the diagonal elements of r in x. 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 ## Eliminate the diagonal matrix d using a givens rotation 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 ## The transformations to eliminate the row of d modify only a 1901 ## single element of (q transpose)*b beyond the first n, which 1902 ## is initially zero. 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 ## Compute the modified diagonal element of r and the 1917 ## modified element of ((q transpose)*b,0). 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 ## Accumulate the transformation in the row of s 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 ## Solve the triangular system for z. If the system is singular 1932 ## then obtain a least squares solution 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] ## Degenerate case 1941 ## *** Reverse loop *** 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 ## Permute the components of z back to components of x 1947 put(x, ipvt, wa) 1948 return(r, x, sdiag)
1949 1950 1951 1952 1953 # Original FORTRAN documentation 1954 # 1955 # subroutine lmpar 1956 # 1957 # given an m by n matrix a, an n by n nonsingular diagonal 1958 # matrix d, an m-vector b, and a positive number delta, 1959 # the problem is to determine a value for the parameter 1960 # par such that if x solves the system 1961 # 1962 # a*x = b , sqrt(par)*d*x = 0 , 1963 # 1964 # in the least squares sense, and dxnorm is the euclidean 1965 # norm of d*x, then either par is zero and 1966 # 1967 # (dxnorm-delta) .le. 0.1*delta , 1968 # 1969 # or par is positive and 1970 # 1971 # abs(dxnorm-delta) .le. 0.1*delta . 1972 # 1973 # this subroutine completes the solution of the problem 1974 # if it is provided with the necessary information from the 1975 # qr factorization, with column pivoting, of a. that is, if 1976 # a*p = q*r, where p is a permutation matrix, q has orthogonal 1977 # columns, and r is an upper triangular matrix with diagonal 1978 # elements of nonincreasing magnitude, then lmpar expects 1979 # the full upper triangle of r, the permutation matrix p, 1980 # and the first n components of (q transpose)*b. on output 1981 # lmpar also provides an upper triangular matrix s such that 1982 # 1983 # t t t 1984 # p *(a *a + par*d*d)*p = s *s . 1985 # 1986 # s is employed within lmpar and may be of separate interest. 1987 # 1988 # only a few iterations are generally needed for convergence 1989 # of the algorithm. if, however, the limit of 10 iterations 1990 # is reached, then the output par will contain the best 1991 # value obtained so far. 1992 # 1993 # the subroutine statement is 1994 # 1995 # subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, 1996 # wa1,wa2) 1997 # 1998 # where 1999 # 2000 # n is a positive integer input variable set to the order of r. 2001 # 2002 # r is an n by n array. on input the full upper triangle 2003 # must contain the full upper triangle of the matrix r. 2004 # on output the full upper triangle is unaltered, and the 2005 # strict lower triangle contains the strict upper triangle 2006 # (transposed) of the upper triangular matrix s. 2007 # 2008 # ldr is a positive integer input variable not less than n 2009 # which specifies the leading dimension of the array r. 2010 # 2011 # ipvt is an integer input array of length n which defines the 2012 # permutation matrix p such that a*p = q*r. column j of p 2013 # is column ipvt(j) of the identity matrix. 2014 # 2015 # diag is an input array of length n which must contain the 2016 # diagonal elements of the matrix d. 2017 # 2018 # qtb is an input array of length n which must contain the first 2019 # n elements of the vector (q transpose)*b. 2020 # 2021 # delta is a positive input variable which specifies an upper 2022 # bound on the euclidean norm of d*x. 2023 # 2024 # par is a nonnegative variable. on input par contains an 2025 # initial estimate of the levenberg-marquardt parameter. 2026 # on output par contains the final estimate. 2027 # 2028 # x is an output array of length n which contains the least 2029 # squares solution of the system a*x = b, sqrt(par)*d*x = 0, 2030 # for the output par. 2031 # 2032 # sdiag is an output array of length n which contains the 2033 # diagonal elements of the upper triangular matrix s. 2034 # 2035 # wa1 and wa2 are work arrays of length n. 2036 # 2037 # subprograms called 2038 # 2039 # minpack-supplied ... dpmpar,enorm,qrsolv 2040 # 2041 # fortran-supplied ... dabs,dmax1,dmin1,dsqrt 2042 # 2043 # argonne national laboratory. minpack project. march 1980. 2044 # burton s. garbow, kenneth e. hillstrom, jorge j. more 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 ## Compute and store in x the gauss-newton direction. If the 2057 ## jacobian is rank-deficient, obtain a least-squares solution 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] #patched from lthresh which doesn't exist WR 2062 if len(wh) > 0: 2063 nsing = wh[0] 2064 wa1[wh[0]:] = 0 2065 if nsing >= 1: 2066 ## *** Reverse loop *** 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 ## Note: ipvt here is a permutation array 2073 put(x, ipvt, wa1) 2074 2075 ## Initialize the iteration counter. Evaluate the function at the 2076 ## origin, and test for acceptance of the gauss-newton direction 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 ## If the jacobian is not rank deficient, the newton step provides a 2085 ## lower bound, parl, for the zero of the function. Otherwise set 2086 ## this bound to zero. 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] ## Degenerate case 2092 for j in range(1,n): ## Note "1" here, not zero 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 ## Calculate an upper bound, paru, for the zero of the function 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 ## If the input par lies outside of the interval (parl,paru), set 2108 ## par to the closer endpoint 2109 2110 par = numpy.max([par,parl]) 2111 par = numpy.min([par,paru]) 2112 if par == 0: par = gnorm/dxnorm 2113 2114 ## Beginning of an interation 2115 while(1): 2116 iter = iter + 1 2117 2118 ## Evaluate the function at the current value of par 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 ## Compute the newton correction 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] ## Degenerate case 2139 2140 temp = self.enorm(wa1) 2141 parc = ((fp/delta)/temp)/temp 2142 2143 ## Depending on the sign of the function, update parl or paru 2144 if fp > 0: parl = numpy.max([parl,par]) 2145 if fp < 0: paru = numpy.min([paru,par]) 2146 2147 ## Compute an improved estimate for par 2148 par = numpy.max([parl, par+parc]) 2149 2150 ## End of an iteration 2151 2152 ## Termination 2153 return[r, par, x, sdiag]
2154 2155 2156 ## Procedure to tie one parameter to another.
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 # Original FORTRAN documentation 2168 # ********** 2169 # 2170 # subroutine covar 2171 # 2172 # given an m by n matrix a, the problem is to determine 2173 # the covariance matrix corresponding to a, defined as 2174 # 2175 # t 2176 # inverse(a *a) . 2177 # 2178 # this subroutine completes the solution of the problem 2179 # if it is provided with the necessary information from the 2180 # qr factorization, with column pivoting, of a. that is, if 2181 # a*p = q*r, where p is a permutation matrix, q has orthogonal 2182 # columns, and r is an upper triangular matrix with diagonal 2183 # elements of nonincreasing magnitude, then covar expects 2184 # the full upper triangle of r and the permutation matrix p. 2185 # the covariance matrix is then computed as 2186 # 2187 # t t 2188 # p*inverse(r *r)*p . 2189 # 2190 # if a is nearly rank deficient, it may be desirable to compute 2191 # the covariance matrix corresponding to the linearly independent 2192 # columns of a. to define the numerical rank of a, covar uses 2193 # the tolerance tol. if l is the largest integer such that 2194 # 2195 # abs(r(l,l)) .gt. tol*abs(r(1,1)) , 2196 # 2197 # then covar computes the covariance matrix corresponding to 2198 # the first l columns of r. for k greater than l, column 2199 # and row ipvt(k) of the covariance matrix are set to zero. 2200 # 2201 # the subroutine statement is 2202 # 2203 # subroutine covar(n,r,ldr,ipvt,tol,wa) 2204 # 2205 # where 2206 # 2207 # n is a positive integer input variable set to the order of r. 2208 # 2209 # r is an n by n array. on input the full upper triangle must 2210 # contain the full upper triangle of the matrix r. on output 2211 # r contains the square symmetric covariance matrix. 2212 # 2213 # ldr is a positive integer input variable not less than n 2214 # which specifies the leading dimension of the array r. 2215 # 2216 # ipvt is an integer input array of length n which defines the 2217 # permutation matrix p such that a*p = q*r. column j of p 2218 # is column ipvt(j) of the identity matrix. 2219 # 2220 # tol is a nonnegative input variable used to define the 2221 # numerical rank of a in the manner described above. 2222 # 2223 # wa is a work array of length n. 2224 # 2225 # subprograms called 2226 # 2227 # fortran-supplied ... dabs 2228 # 2229 # argonne national laboratory. minpack project. august 1980. 2230 # burton s. garbow, kenneth e. hillstrom, jorge j. more 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 ## For the inverse of r in the full upper triangle of r 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 ## Form the full upper triangle of the inverse of (r transpose)*r 2263 ## in the full upper triangle of r 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 ## For the full lower triangle of the covariance matrix 2273 ## in the strict lower triangle or and in wa 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 ## Symmetrize the covariance matrix in r 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
2292 -class machar:
2293 - def __init__(self, double=1):
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