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