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
  11 # modify it under the terms of the GNU Lesser General Public
  12 # License as published by the Free Software Foundation, either
  13 # version 3 of the License, or (at your option) any later version.
  14 #
  15 # Hooke is distributed in the hope that it will be useful,
  16 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18 # GNU Lesser General 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.peak import find_peaks, find_peaks_arguments, Peak, _kwargs
  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', 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(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, 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         for col in ['surface z piezo (m)', 'deflection (N)']:
 156             if col not in curve.data[0].info['columns']:
 157                 scale(curve)
 158         data = None
 159         for block in curve.data:
 160             if block.info['name'].startswith('retract'):
 161                 data = block
 162                 break
 163         if data == None:
 164             raise Failure('No retraction blocks in %s.' % curve)
 165         z_data = data[:,data.info['columns'].index('surface z piezo (m)')]
 166         if 'flattened deflection (N)' in data.info['columns']:
 167             params['deflection column name'] = 'flattened deflection (N)'
 168         else:
 169             params['deflection column name'] = 'deflection (N)'
 170         d_data = data[:,data.info['columns'].index(
 171                 params['deflection column name'])]
 172         for key,value in params.items():
 173             if value == None: # Use configured default value.
 174                 params[key] = self.plugin.config[key]
 175         return z_data,d_data,params
 176
 177
 178 class ConvolutionFilterCommand (FilterCommand):
 179     u"""Create a subset playlist of curves with enough convolution peaks.
 180
 181     Notes
 182     -----
 183     This filter can reduce a dataset like the one in [#brucale2009] to
 184     analyze by hand by about 80-90% (depending on the overall
 185     cleanliness of the data set). Thousands of curves can be
 186     automatically filtered this way in a few minutes on a standard PC,
 187     but the algorithm could still be optimized.
 188
 189     .. [#brucale2009] M. Brucale, M. Sandal, S. Di Maio, A. Rampioni,
 190       I. Tessari, L. Tosatto, M. Bisaglia, L. Bubacco, B. Samorì.
 191       "Pathogenic mutations shift the equilibria of
 192       α-Synuclein single molecules towards structured
 193       conformers."
 194       Chembiochem., 2009.
 195       doi: `10.1002/cbic.200800581 <http://dx.doi.org/10.1002/cbic.200800581>`_
 196
 197     See Also
 198     --------
 199     ConvolutionCommand : Underlying convolution-based peak detection.
 200     """
 201     def __init__(self, plugin):
 202         super(ConvolutionFilterCommand, self).__init__(
 203             plugin, name='convolution filter playlist')
 204         self.arguments.extend(plugin._arguments)
 205
 206     def filter(self, curve, hooke, inqueue, outqueue, params):
 207         inq = Queue()
 208         outq = Queue()
 209         conv_command = [c for c in self.hooke.commands
 210                         if c.name=='convolution peaks'][0]
 211         conv_command.run(hooke, inq, outq, **params)
 212         peaks = outq.get()
 213         if not (isinstance(peaks, list) and (len(peaks) == 0
 214                                              or isinstance(peaks[0], Peak))):
 215             raise Failure('Expected a list of Peaks, not %s' % peaks)
 216         ret = outq.get()
 217         if not isinstance(ret, Success):
 218             raise ret
 219         return len(peaks) >= params['min peaks']