From ad1d75adfdae863cfe0c58148665a550cf7f7690 Mon Sep 17 00:00:00 2001
From: "W. Trevor King" <wking@drexel.edu>
Date: Wed, 11 Aug 2010 09:11:36 -0400
Subject: [PATCH] Oops, correct my faulty interpretation of leastsq's 'diag'
 argument in fit.py

---
 hooke/util/fit.py | 60 +++++++++++++++++++++++++++++++++--------------
 1 file changed, 43 insertions(+), 17 deletions(-)

diff --git a/hooke/util/fit.py b/hooke/util/fit.py
index efea018..1360677 100644
--- a/hooke/util/fit.py
+++ b/hooke/util/fit.py
@@ -42,7 +42,7 @@ class ModelFitter (object):
         Store any extra information useful inside your overridden
         methods.
     rescale : boolean
-        Rescale parameters so the scale of each is 1.0.  Also rescale
+        Rescale parameters so the guess for each is 1.0.  Also rescale
         the data so data.std() == 1.0.
 
     Examples
@@ -76,14 +76,13 @@ class ModelFitter (object):
     ...         return [float(self._data[-1] - self._data[0])/len(self._data),
     ...                 self._data[0]]
     ...     def guess_scale(self, params, outqueue=None):
-    ...         slope_scale = params[0]/10.
-    ...         if slope_scale == 0:  # data is expected to be flat
-    ...             slope_scale = float(self._data.max()-self._data.min())/len(self._data)
-    ...             if slope_scale == 0:  # data is completely flat
-    ...                 slope_scale = 1.
-    ...         offset_scale = self._data.std()/10.0
-    ...         if offset_scale == 0:  # data is completely flat
-    ...             offset_scale = 1.
+    ...         slope_scale = 0.1
+    ...         if params[1] == 0:
+    ...             offset_scale = 1
+    ...         else:
+    ...             offset_scale = 0.1*self._data.std()/abs(params[1])
+    ...             if offset_scale == 0:  # data is completely flat
+    ...                 offset_scale = 1.
     ...         return [slope_scale, offset_scale]
     >>> data = 20*numpy.sin(arange(1000)) + 7.*arange(1000) - 33.0
     >>> m = LinearModel(data)
@@ -93,7 +92,6 @@ class ModelFitter (object):
     >>> pprint(info)  # doctest: +ELLIPSIS, +REPORT_UDIFF
     {'active fitted parameters': array([  6.999..., -32.889...]),
      'active parameters': array([  6.999..., -32.889...]),
-     'active scale': [0.699..., 202.071...],
      'convergence flag': 2,
      'covariance matrix': array([[  1.199...e-08,  -5.993...e-06],
            [ -5.993...e-06,   3.994...e-03]]),
@@ -107,7 +105,7 @@ class ModelFitter (object):
      'initial parameters': [6.992..., -33.0],
      'message': 'The relative error between two consecutive iterates is at most 0.000...',
      'rescaled': False,
-     'scale': [0.699..., 202.071...]}
+     'scale': [0.100..., 6.123...]}
 
     We round the outputs to protect the doctest against differences in
     machine rounding during computation.  We expect the values to be close
@@ -156,6 +154,33 @@ class ModelFitter (object):
         return []
 
     def guess_scale(self, params, outqueue=None):
+        """Guess the problem length scale in each parameter dimension.
+
+        Notes
+        -----
+        From the :func:`scipy.optimize.leastsq` documentation, `diag`
+        (which we refer to as `scale`, sets `factor * || diag * x||`
+        as the initial step.  If `x == 0`, then `factor` is used
+        instead (from `lmdif.f`_)::
+
+                      do 70 j = 1, n
+                        wa3(j) = diag(j)*x(j)
+               70       continue
+                      xnorm = enorm(n,wa3)
+                      delta = factor*xnorm
+                      if (delta .eq. zero) delta = factor
+  
+        For most situations then, you don't need to do anything fancy.
+        The default scaling (if you don't set a scale) is::
+
+            c        on the first iteration and if mode is 1, scale according
+            c        to the norms of the columns of the initial jacobian.
+
+        (i.e. `diag(j) = acnorm(j)`, where `acnorm(j) is the norm of the `j`th column
+        of the initial Jacobian).
+
+        .. _lmdif.f: http://www.netlib.org/minpack/lmdif.f
+        """
         return None
 
     def residual(self, params):
@@ -188,18 +213,20 @@ class ModelFitter (object):
             initial_params = self.guess_initial_params(outqueue)
         if scale == None:
             scale = self.guess_scale(initial_params, outqueue)
-        assert min(scale) > 0, scale
+        if scale != None:
+            assert min(scale) > 0, scale
         if self._rescale == True:
-            self._param_scale_factors = scale
-            active_scale = [1.0 for s in scale]
+            self._param_scale_factors = initial_params
+            for i,s in enumerate(self._param_scale_factors):
+                if s == 0:
+                    self._param_scale_factors[i] = 1.0
             active_params = [p/s for p,s in zip(initial_params,
                                                 self._param_scale_factors)]
         else:
-            active_scale = scale
             active_params = initial_params
         params,cov,info,mesg,ier = leastsq(
             func=self.residual, x0=active_params, full_output=True,
-            diag=active_scale, **kwargs)
+            diag=scale, **kwargs)
         if self._rescale == True:
             active_params = params
             if len(initial_params) == 1:  # params is a float
@@ -215,7 +242,6 @@ class ModelFitter (object):
                     'initial parameters': initial_params,
                     'active parameters': active_params,
                     'scale': scale,
-                    'active scale': active_scale,
                     'data scale factor': self._data_scale_factor,
                     'active fitted parameters': active_params,
                     'fitted parameters': params,
-- 
2.26.2