Update SurfacePositionModel to match new diag/scale interpretation.

diff --git a/hooke/plugin/vclamp.py b/hooke/plugin/vclamp.py
index 459160372eba75f03b6d52a58bccd9ec276aff4e..78a46147f96277be9fab2a45153ee28512657320 100644 (file)
--- a/hooke/plugin/vclamp.py
+++ b/hooke/plugin/vclamp.py
@@ -76,6 +76,7 @@ def scale(hooke, curve, block=None):
             cant_adjust.run(hooke, inqueue, outqueue, **params)
     return curve
 
+
 class SurfacePositionModel (ModelFitter):
     """Bilinear surface position model.
 
@@ -97,7 +98,9 @@ class SurfacePositionModel (ModelFitter):
     In order for this model to produce a satisfactory fit, there
     should be enough data in the off-surface region that interactions
     due to proteins, etc. will not seriously skew the fit in the
-    off-surface region.
+    off-surface region.  If you don't have much of a tail, you can set
+    the `info` dict's `ignore non-contact before index` parameter to
+    the index of the last surface- or protein-related feature.
     """
     def model(self, params):
         """A continuous, bilinear model.
@@ -116,7 +119,8 @@ class SurfacePositionModel (ModelFitter):
         :math:`p_3` is the slope of the second region.
         """
         p = params  # convenient alias
-        if self.info.get('force zero non-contact slope', None) == True:
+        rNC_ignore = self.info['ignore non-contact before index']
+        if self.info['force zero non-contact slope'] == True:
             p = list(p)
             p.append(0.)  # restore the non-contact slope parameter
         r2 = numpy.round(abs(p[2]))
@@ -125,19 +129,32 @@ class SurfacePositionModel (ModelFitter):
         if r2 < len(self._data)-1:
             self._model_data[r2:] = \
                 p[0] + p[1]*p[2] + p[3] * numpy.arange(len(self._data)-r2)
+            if r2 < rNC_ignore:
+                self._model_data[r2:rNC_ignore] = self._data[r2:rNC_ignore]
         return self._model_data
 
-    def set_data(self, data, info=None):
-        super(SurfacePositionModel, self).set_data(data, info)
-        if info == None:
-            info = {}
-        self.info = info
-        self.info['min position'] = 0
-        self.info['max position'] = len(data)
-        self.info['max deflection'] = data.max()
-        self.info['min deflection'] = data.min()
-        self.info['position range'] = self.info['max position'] - self.info['min position']
-        self.info['deflection range'] = self.info['max deflection'] - self.info['min deflection']
+    def set_data(self, data, info=None, *args, **kwargs):
+        super(SurfacePositionModel, self).set_data(data, info, *args, **kwargs)
+        if self.info == None:
+            self.info = {}
+        for key,value in [
+            ('force zero non-contact slope', False),
+            ('ignore non-contact before index', 6158),
+            ('min position', 0),  # Store postions etc. to avoid recalculating.
+            ('max position', len(data)),
+            ('max deflection', data.max()),
+            ('min deflection', data.min()),
+            ]:
+            if key not in self.info:
+                self.info[key] = value
+        for key,value in [
+            ('position range',
+             self.info['max position'] - self.info['min position']),
+            ('deflection range',
+             self.info['max deflection'] - self.info['min deflection']),
+            ]:
+            if key not in self.info:
+                self.info[key] = value
 
     def guess_initial_params(self, outqueue=None):
         """Guess the initial parameters.
@@ -158,7 +175,7 @@ class SurfacePositionModel (ModelFitter):
         right_slope = 0
         self.info['guessed contact slope'] = left_slope
         params = [left_offset, left_slope, kink_position, right_slope]
-        if self.info.get('force zero non-contact slope', None) == True:
+        if self.info['force zero non-contact slope'] == True:
             params = params[:-1]
         return params
 
@@ -167,27 +184,32 @@ class SurfacePositionModel (ModelFitter):
 
         Notes
         -----
-        We guess offset scale (:math:`p_0`) as one tenth of the total
-        deflection range, the kink scale (:math:`p_2`) as one tenth of
-        the total index range, the initial (contact) slope scale
-        (:math:`p_1`) as one tenth of the contact slope estimation,
-        and the final (non-contact) slope scale (:math:`p_3`) is as
-        one tenth of the initial slope scale.
+        We the scale as one tenth for each parameter.
         """
-        offset_scale = self.info['deflection range']/10.
-        left_slope_scale = abs(params[1])/10.
-        kink_scale = self.info['position range']/10.
-        right_slope_scale = left_slope_scale/10.
+        offset_scale = 0.1
+        left_slope_scale = 0.1
+        kink_scale = 0.1
+        right_slope_scale = 0.1
         scale = [offset_scale, left_slope_scale, kink_scale, right_slope_scale]
-        if self.info.get('force zero non-contact slope', None) == True:
+        if self.info['force zero non-contact slope'] == True:
             scale = scale[:-1]
         return scale
 
     def fit(self, *args, **kwargs):
+        """Fit the model to the data.
+
+        Notes
+        -----
+        We change the `epsfcn` default from :func:`scipy.optimize.leastsq`'s
+        `0` to `1e-3`, so the initial Jacobian estimate takes larger steps,
+        which helps avoid being trapped in noise-generated local minima.
+        """
         self.info['guessed contact slope'] = None
+        if 'epsfcn' not in kwargs:
+            kwargs['epsfcn'] = 1e-3  # take big steps to estimate the Jacobian
         params = super(SurfacePositionModel, self).fit(*args, **kwargs)
         params[2] = abs(params[2])
-        if self.info.get('force zero non-contact slope', None) == True:
+        if self.info['force zero non-contact slope'] == True:
             params = list(params)
             params.append(0.)  # restore the non-contact slope parameter
 
@@ -211,6 +233,7 @@ class SurfacePositionModel (ModelFitter):
                           % (params[3], self.info['guessed contact slope']))
         return params
 
+
 class VelocityClampPlugin (Plugin):
     def __init__(self):
         super(VelocityClampPlugin, self).__init__(name='vclamp')
@@ -442,13 +465,14 @@ Name (without units) for storing fit parameters in the `.info` dictionary.
 
         Notes
         -----
-        Uses :func:`analyze_surf_pos_data_wtk` internally.
+        Uses :class:`SurfacePositionModel` internally.
         """
         reverse = z_data[0] > z_data[-1]
         if reverse == True:    # approaching, contact region on the right
             d_data = d_data[::-1]
-        s = SurfacePositionModel(d_data)
-        s.info['force zero non-contact slope'] = True
+        s = SurfacePositionModel(d_data, info={
+                'force zero non-contact slope':True},
+                                 rescale=True)
         offset,contact_slope,surface_index,non_contact_slope = s.fit(
             outqueue=outqueue)
         info = {