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