hooke/plugin/convfilt.py

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