hooke/plugin/convfilt.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2010 Alberto Gomez-Casado
   4 #                         Fabrizio Benedetti
   5 #                         Massimo Sandal <devicerandom@gmail.com>
   6 #                         W. Trevor King <wking@drexel.edu>
   7 #
   8 # This file is part of Hooke.
   9 #
  10 # Hooke is free software: you can redistribute it and/or modify it
  11 # under the terms of the GNU Lesser General Public License as
  12 # published by the Free Software Foundation, either version 3 of the
  13 # License, or (at your option) any later version.
  14 #
  15 # Hooke is distributed in the hope that it will be useful, but WITHOUT
  16 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  17 # or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  18 # Public License for more details.
  19 #
  20 # You should have received a copy of the GNU Lesser General Public
  21 # License along with Hooke.  If not, see
  22 # <http://www.gnu.org/licenses/>.
  23
  24 """The ``convfilt`` module provides :class:`ConvFiltPlugin` and
  25 for convolution-based filtering of :mod:`hooke.curve.Curve`\s.
  26
  27 See Also
  28 --------
  29 :mod:`~hooke.plugin.flatfilt` for a collection of additional filters.
  30
  31 :class:`ConfFiltPlugin` was broken out of
  32 :class:`~hooke.plugin.flatfilt` to separate its large number of
  33 configuration settings.
  34 """
  35
  36 import copy
  37 from multiprocessing import Queue
  38
  39 import numpy
  40
  41 from ..command import Command, Argument, Success, Failure
  42 from ..config import Setting
  43 from ..experiment import VelocityClamp
  44 from ..plugin import Plugin, argument_to_setting
  45 from ..plugin.curve import CurveArgument
  46 from ..plugin.playlist import FilterCommand
  47 from ..plugin.vclamp import scale
  48 from ..util.fit import PoorFit
  49 from ..util.peak import find_peaks, find_peaks_arguments, Peak, _kwargs
  50
  51
  52 class ConvFiltPlugin (Plugin):
  53     """Convolution-based peak recognition and filtering.
  54     """
  55     def __init__(self):
  56         super(ConvFiltPlugin, self).__init__(name='convfilt')
  57         self._arguments = [ # for Command initialization
  58             Argument('convolution', type='float', count=-1,
  59                      default=[11.0]+[-1.0]*11, help="""
  60 Convolution vector roughly matching post-peak cantilever rebound.
  61 This should roughly match the shape of the feature you're looking for.
  62 """.strip()), # TODO: per-curve convolution vector + interpolation, to
  63               # take advantage of the known spring constant.
  64             Argument('blind window', type='float', default=20e-9, help="""
  65 Meters after the contact point where we do not count peaks to avoid
  66 non-specific surface interaction.
  67 """.strip()),
  68             Argument('min peaks', type='int', default=5, help="""
  69 Minimum number of peaks for curve acceptance.
  70 """.strip()),
  71             ] + copy.deepcopy(find_peaks_arguments)
  72         # Set convolution-specific defaults for the fit_peak_arguments.
  73         for key,value in [('cut side', 'positive'),
  74                           ('stable', 0.005),
  75                           ('max cut', 0.2),
  76                           ('min deviations', 5.0),
  77                           ('min points', 1),
  78                           ('see double', 10), # TODO: points vs. meters. 10e-9),
  79                           ]:
  80             argument = [a for a in self._arguments if a.name == key][0]
  81             argument.default = value
  82         self._settings = [
  83             Setting(section=self.setting_section, help=self.__doc__)]
  84         for argument in self._arguments:
  85             self._settings.append(argument_to_setting(
  86                     self.setting_section, argument))
  87             argument.default = None # if argument isn't given, use the config.
  88         self._commands = [ConvolutionPeaksCommand(self),
  89                           ConvolutionFilterCommand(self)]
  90
  91     def dependencies(self):
  92         return ['vclamp']
  93
  94     def default_settings(self):
  95         return self._settings
  96
  97
  98 class ConvolutionPeaksCommand (Command):
  99     """Detect peaks in velocity clamp data with a convolution.
 100
 101     Notes
 102     -----
 103     A simplified version of Kasas' filter [#kasas2000].
 104     For any given retracting curve:
 105
 106     1) The contact point is found.
 107     2) The off-surface data is convolved (using :func:`numpy.convolve`)
 108       with a vector that encodes the approximate shape of the target
 109       feature.
 110     3) Peaks in the convolved curve are extracted with
 111       :func:`~hooke.plugins.peak.find_peaks`.
 112
 113     The convolution algorithm, with appropriate thresholds, usually
 114     recognizes peaks well more than 95% of the time.
 115
 116     .. [#kasas2000] S. Kasas, B.M. Riederer, S. Catsicas, B. Cappella,
 117       G. Dietler.
 118       "Fuzzy logic algorithm to extract specific interaction forces
 119       from atomic force microscopy data"
 120       Rev. Sci. Instrum., 2000.
 121       doi: `10.1063/1.1150583 <http://dx.doi.org/10.1063/1.1150583>`_
 122     """
 123     def __init__(self, plugin):
 124         config_arguments = [a for a in plugin._arguments
 125                             if a.name != 'min peaks']
 126         # Drop min peaks, since we're not filtering with this
 127         # function, just detecting peaks.
 128         super(ConvolutionPeaksCommand, self).__init__(
 129             name='convolution peaks',
 130             arguments=[CurveArgument] + config_arguments,
 131             help=self.__doc__, plugin=plugin)
 132
 133     def _run(self, hooke, inqueue, outqueue, params):
 134         z_data,d_data,params = self._setup(hooke, params)
 135         start_index = 0
 136         while z_data[start_index] < params['bind window']:
 137             start_index += 1
 138         conv = numpy.convolve(d_data[start_index:], params['convolution'],
 139                               mode='valid')
 140         peaks = find_peaks(conv, **_kwargs(params, find_peaks_arguments,
 141                                             argument_input_keys=True))
 142         for peak in peaks:
 143             peak.name = 'convolution of %s with %s' \
 144                 % (params['deflection column name'], params['convolution'])
 145             peak.index += start_index
 146         outqueue.put(peaks)
 147
 148     def _setup(self, hooke, params):
 149         """Setup `params` from config and return the z piezo and
 150         deflection arrays.
 151         """
 152         curve = params['curve']
 153         if curve.info['experiment'] != VelocityClamp:
 154             raise Failure('%s operates on VelocityClamp experiments, not %s'
 155                           % (self.name, curve.info['experiment']))
 156         for col in ['surface distance (m)', 'deflection (N)']:
 157             if col not in curve.data[0].info['columns']:
 158                 scale(hooke, curve)
 159         data = None
 160         for block in curve.data:
 161             if block.info['name'].startswith('retract'):
 162                 data = block
 163                 break
 164         if data == None:
 165             raise Failure('No retraction blocks in %s.' % curve)
 166         z_data = data[:,data.info['columns'].index('surface distance (m)')]
 167         if 'flattened deflection (N)' in data.info['columns']:
 168             params['deflection column name'] = 'flattened deflection (N)'
 169         else:
 170             params['deflection column name'] = 'deflection (N)'
 171         d_data = data[:,data.info['columns'].index(
 172                 params['deflection column name'])]
 173         for key,value in params.items():
 174             if value == None: # Use configured default value.
 175                 params[key] = self.plugin.config[key]
 176         return z_data,d_data,params
 177
 178
 179 class ConvolutionFilterCommand (FilterCommand):
 180     u"""Create a subset playlist of curves with enough convolution peaks.
 181
 182     Notes
 183     -----
 184     This filter can reduce a dataset like the one in [#brucale2009] to
 185     analyze by hand by about 80-90% (depending on the overall
 186     cleanliness of the data set). Thousands of curves can be
 187     automatically filtered this way in a few minutes on a standard PC,
 188     but the algorithm could still be optimized.
 189
 190     .. [#brucale2009] M. Brucale, M. Sandal, S. Di Maio, A. Rampioni,
 191       I. Tessari, L. Tosatto, M. Bisaglia, L. Bubacco, B. Samorì.
 192       "Pathogenic mutations shift the equilibria of
 193       α-Synuclein single molecules towards structured
 194       conformers."
 195       Chembiochem., 2009.
 196       doi: `10.1002/cbic.200800581 <http://dx.doi.org/10.1002/cbic.200800581>`_
 197
 198     See Also
 199     --------
 200     ConvolutionCommand : Underlying convolution-based peak detection.
 201     """
 202     def __init__(self, plugin):
 203         super(ConvolutionFilterCommand, self).__init__(
 204             plugin, name='convolution filter playlist')
 205         self.arguments.extend(plugin._arguments)
 206
 207     def filter(self, curve, hooke, inqueue, outqueue, params):
 208         params['curve'] = curve
 209         inq = Queue()
 210         outq = Queue()
 211         conv_command = [c for c in self.hooke.commands
 212                         if c.name=='convolution peaks'][0]
 213         conv_command.run(hooke, inq, outq, **params)
 214         peaks = outq.get()
 215         if isinstance(peaks, UncaughtException) \
 216                 and isinstance(peaks.exception, PoorFit):
 217             return False
 218         if not (isinstance(peaks, list) and (len(peaks) == 0
 219                                              or isinstance(peaks[0], Peak))):
 220             raise Failure('Expected a list of Peaks, not %s' % peaks)
 221         ret = outq.get()
 222         if not isinstance(ret, Success):
 223             raise ret
 224         if params['min peaks'] == None: # Use configured default value.
 225             params['min peaks'] = self.plugin.config['min peaks']
 226         return len(peaks) >= params['min peaks']