From: W. Trevor King Date: Wed, 11 Aug 2010 13:12:27 +0000 (-0400) Subject: Update SurfacePositionModel to match new diag/scale interpretation. X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=9f844eb89787f4b089cff3c7101b3b0dc3091519;p=hooke.git Update SurfacePositionModel to match new diag/scale interpretation. --- diff --git a/hooke/plugin/vclamp.py b/hooke/plugin/vclamp.py index 4591603..78a4614 100644 --- 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 = {