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 ..util.fit import PoorFit
  44 from ..util.peak import find_peaks, find_peaks_arguments, Peak, _kwargs
  45 from . import Plugin, argument_to_setting
  46 from .curve import CurveArgument
  47 from .playlist import FilterCommand
  48
  49
  50 class ConvFiltPlugin (Plugin):
  51     """Convolution-based peak recognition and filtering.
  52     """
  53     def __init__(self):
  54         super(ConvFiltPlugin, self).__init__(name='convfilt')
  55         self._arguments = [ # for Command initialization
  56             Argument('convolution', type='float', count=-1,
  57                      default=[11.0]+[-1.0]*11, help="""
  58 Convolution vector roughly matching post-peak cantilever rebound.
  59 This should roughly match the shape of the feature you're looking for.
  60 """.strip()), # TODO: per-curve convolution vector + interpolation, to
  61               # take advantage of the known spring constant.
  62             Argument('blind window', type='float', default=20e-9, help="""
  63 Meters after the contact point where we do not count peaks to avoid
  64 non-specific surface interaction.
  65 """.strip()),
  66             Argument('min peaks', type='int', default=5, help="""
  67 Minimum number of peaks for curve acceptance.
  68 """.strip()),
  69             ] + copy.deepcopy(find_peaks_arguments)
  70         # Set convolution-specific defaults for the fit_peak_arguments.
  71         for key,value in [('cut side', 'positive'),
  72                           ('stable', 0.005),
  73                           ('max cut', 0.2),
  74                           ('min deviations', 5.0),
  75                           ('min points', 1),
  76                           ('see double', 10), # TODO: points vs. meters. 10e-9),
  77                           ]:
  78             argument = [a for a in self._arguments if a.name == key][0]
  79             argument.default = value
  80         self._settings = [
  81             Setting(section=self.setting_section, help=self.__doc__)]
  82         for argument in self._arguments:
  83             self._settings.append(argument_to_setting(
  84                     self.setting_section, argument))
  85             argument.default = None # if argument isn't given, use the config.
  86         self._commands = [ConvolutionPeaksCommand(self),
  87                           ConvolutionFilterCommand(self)]
  88
  89     def dependencies(self):
  90         return ['vclamp']
  91
  92     def default_settings(self):
  93         return self._settings
  94
  95
  96 class ConvolutionPeaksCommand (Command):
  97     """Detect peaks in velocity clamp data with a convolution.
  98
  99     Notes
 100     -----
 101     A simplified version of Kasas' filter [#kasas2000].
 102     For any given retracting curve:
 103
 104     1) The contact point is found.
 105     2) The off-surface data is convolved (using :func:`numpy.convolve`)
 106       with a vector that encodes the approximate shape of the target
 107       feature.
 108     3) Peaks in the convolved curve are extracted with
 109       :func:`~hooke.plugins.peak.find_peaks`.
 110
 111     The convolution algorithm, with appropriate thresholds, usually
 112     recognizes peaks well more than 95% of the time.
 113
 114     .. [#kasas2000] S. Kasas, B.M. Riederer, S. Catsicas, B. Cappella,
 115       G. Dietler.
 116       "Fuzzy logic algorithm to extract specific interaction forces
 117       from atomic force microscopy data"
 118       Rev. Sci. Instrum., 2000.
 119       doi: `10.1063/1.1150583 <http://dx.doi.org/10.1063/1.1150583>`_
 120     """
 121     def __init__(self, plugin):
 122         config_arguments = [a for a in plugin._arguments
 123                             if a.name != 'min peaks']
 124         # Drop min peaks, since we're not filtering with this
 125         # function, just detecting peaks.
 126         super(ConvolutionPeaksCommand, self).__init__(
 127             name='convolution peaks',
 128             arguments=[CurveArgument] + config_arguments,
 129             help=self.__doc__, plugin=plugin)
 130
 131     def _run(self, hooke, inqueue, outqueue, params):
 132         z_data,d_data,params = self._setup(hooke, params)
 133         start_index = 0
 134         while z_data[start_index] < params['bind window']:
 135             start_index += 1
 136         conv = numpy.convolve(d_data[start_index:], params['convolution'],
 137                               mode='valid')
 138         peaks = find_peaks(conv, **_kwargs(params, find_peaks_arguments,
 139                                             argument_input_keys=True))
 140         for peak in peaks:
 141             peak.name = 'convolution of %s with %s' \
 142                 % (params['deflection column name'], params['convolution'])
 143             peak.index += start_index
 144         outqueue.put(peaks)
 145
 146     def _setup(self, hooke, params):
 147         """Setup `params` from config and return the z piezo and
 148         deflection arrays.
 149         """
 150         curve = params['curve']
 151         data = None
 152         for block in curve.data:
 153             if block.info['name'].startswith('retract'):
 154                 data = block
 155                 break
 156         if data == None:
 157             raise Failure('No retraction blocks in %s.' % curve)
 158         z_data = data[:,data.info['columns'].index('surface distance (m)')]
 159         if 'flattened deflection (N)' in data.info['columns']:
 160             params['deflection column name'] = 'flattened deflection (N)'
 161         else:
 162             params['deflection column name'] = 'deflection (N)'
 163         d_data = data[:,data.info['columns'].index(
 164                 params['deflection column name'])]
 165         for key,value in params.items():
 166             if value == None: # Use configured default value.
 167                 params[key] = self.plugin.config[key]
 168         return z_data,d_data,params
 169
 170
 171 class ConvolutionFilterCommand (FilterCommand):
 172     u"""Create a subset playlist of curves with enough convolution peaks.
 173
 174     Notes
 175     -----
 176     This filter can reduce a dataset like the one in [#brucale2009] to
 177     analyze by hand by about 80-90% (depending on the overall
 178     cleanliness of the data set). Thousands of curves can be
 179     automatically filtered this way in a few minutes on a standard PC,
 180     but the algorithm could still be optimized.
 181
 182     .. [#brucale2009] M. Brucale, M. Sandal, S. Di Maio, A. Rampioni,
 183       I. Tessari, L. Tosatto, M. Bisaglia, L. Bubacco, B. Samorì.
 184       "Pathogenic mutations shift the equilibria of
 185       α-Synuclein single molecules towards structured
 186       conformers."
 187       Chembiochem., 2009.
 188       doi: `10.1002/cbic.200800581 <http://dx.doi.org/10.1002/cbic.200800581>`_
 189
 190     See Also
 191     --------
 192     ConvolutionCommand : Underlying convolution-based peak detection.
 193     """
 194     def __init__(self, plugin):
 195         super(ConvolutionFilterCommand, self).__init__(
 196             plugin, name='convolution filter playlist')
 197         self.arguments.extend(plugin._arguments)
 198
 199     def filter(self, curve, hooke, inqueue, outqueue, params):
 200         params['curve'] = curve
 201         inq = Queue()
 202         outq = Queue()
 203         conv_command = [c for c in self.hooke.commands
 204                         if c.name=='convolution peaks'][0]
 205         conv_command.run(hooke, inq, outq, **params)
 206         peaks = outq.get()
 207         if isinstance(peaks, UncaughtException) \
 208                 and isinstance(peaks.exception, PoorFit):
 209             return False
 210         if not (isinstance(peaks, list) and (len(peaks) == 0
 211                                              or isinstance(peaks[0], Peak))):
 212             raise Failure('Expected a list of Peaks, not %s' % peaks)
 213         ret = outq.get()
 214         if not isinstance(ret, Success):
 215             raise ret
 216         if params['min peaks'] == None: # Use configured default value.
 217             params['min peaks'] = self.plugin.config['min peaks']
 218         return len(peaks) >= params['min peaks']