hooke/plugin/flatfilt.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 ``flatfilt`` module provides :class:`FlatFiltPlugin` and
  22 several associated :class:`~hooke.command.Command`\s for removing flat
  23 (featureless) :mod:`~hooke.curve.Curve`\s from
  24 :class:`~hooke.playlist.Playlist`\s.
  25
  26 See Also
  27 --------
  28 :mod:`~hooke.plugin.convfilt` for a convolution-based filter.
  29 """
  30
  31 import copy
  32 from Queue import Queue
  33
  34 from numpy import absolute, diff
  35 from scipy.signal.signaltools import medfilt
  36
  37 from ..command import Argument, Success, Failure, UncaughtException
  38 from ..config import Setting
  39 from ..curve import Data
  40 from ..util.fit import PoorFit
  41 from ..util.peak import (find_peaks, peaks_to_mask,
  42                          find_peaks_arguments, Peak, _kwargs)
  43 from ..util.si import join_data_label, split_data_label
  44 from . import Plugin, argument_to_setting
  45 from .curve import ColumnAddingCommand
  46 from .playlist import FilterCommand
  47
  48
  49 class FlatFiltPlugin (Plugin):
  50     """Standard-devitiation-based peak recognition and filtering.
  51     """
  52     def __init__(self):
  53         super(FlatFiltPlugin, self).__init__(name='flatfilt')
  54         self._arguments = [ # For Command initialization
  55             Argument('median window', type='int', default=7, help="""
  56 Median window filter size (in points).
  57 """.strip()),
  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=4, help="""
  63 Minimum number of peaks for curve acceptance.
  64 """.strip()),
  65             ] + copy.deepcopy(find_peaks_arguments)
  66         # Set flat-filter-specific defaults for the fit_peak_arguments.
  67         for key,value in [('cut side', 'both'),
  68                           ('stable', 0.005),
  69                           ('max cut', 0.2),
  70                           ('min deviations', 9.0),
  71                           ('min points', 4),
  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 = [FlatPeaksCommand(self)]
  83         # append FlatFilterCommand so it can steal arguments from
  84         # FlatPeaksCommand.
  85         self._commands.append(FlatFilterCommand(self))
  86
  87     def dependencies(self):
  88         return ['vclamp']
  89
  90     def default_settings(self):
  91         return self._settings
  92
  93
  94 class FlatPeaksCommand (ColumnAddingCommand):
  95     """Detect peaks in velocity clamp data using noise statistics.
  96
  97     Notes
  98     -----
  99     Noise analysis on the retraction curve:
 100
 101     1) A median window filter (using
 102       :func:`scipy.signal.signaltools.medfilt`) smooths the
 103       deflection.
 104     2) The deflection derivative is calculated (using
 105       :func:`numpy.diff` which uses forward differencing).
 106     3) Peaks in the derivative curve are extracted with
 107       :func:`~hooke.plugins.peak.find_peaks`.
 108
 109     The algorithm was originally Francesco Musiani's idea.
 110     """
 111     def __init__(self, plugin):
 112         plugin_arguments = [a for a in plugin._arguments
 113                             if a.name != 'min peaks']
 114         # Drop min peaks, since we're not filtering with this
 115         # function, just detecting peaks.
 116         super(FlatPeaksCommand, self).__init__(
 117             name='flat filter peaks',
 118             columns=[
 119                 ('distance column', 'surface distance (m)', """
 120 Name of the column to use as the surface position input.
 121 """.strip()),
 122                 ('deflection column', 'surface deflection (m)', """
 123 Name of the column to use as the deflection input.
 124 """.strip()),
 125                 ],
 126             new_columns=[
 127                 ('output peak column', 'flat filter peaks', """
 128 Name of the column (without units) to use as the peak output.
 129 """.strip()),
 130                 ],
 131             arguments=[
 132                 Argument(name='peak info name', type='string',
 133                          default='flat filter peaks',
 134                          help="""
 135 Name for storing the list of peaks in the `.info` dictionary.
 136 """.strip()),
 137                 ] + plugin_arguments,
 138             help=self.__doc__, plugin=plugin)
 139
 140     def _run(self, hooke, inqueue, outqueue, params):
 141         self._add_to_command_stack(params)
 142         params = self._setup_params(hooke=hooke, params=params)
 143         block = self._block(hooke=hooke, params=params)
 144         dist_data = self._get_column(hooke=hooke, params=params,
 145                                      column_name='distance column')
 146         def_data = self._get_column(hooke=hooke, params=params,
 147                                      column_name='deflection column')
 148         start_index = absolute(dist_data-params['blind window']).argmin()
 149         median = medfilt(def_data[start_index:], params['median window'])
 150         deriv = diff(median)
 151         peaks = find_peaks(deriv, **_kwargs(params, find_peaks_arguments,
 152                                             argument_input_keys=True))
 153         for i,peak in enumerate(peaks):
 154             peak.name = self._peak_name(params, i)
 155             peak.index += start_index
 156         block.info[params['peak info name']] = peaks
 157
 158         self._set_column(hooke=hooke, params=params,
 159                          column_name='output peak column',
 160                          values=peaks_to_mask(def_data, peaks) * def_data)
 161         outqueue.put(peaks)
 162
 163     def _setup_params(self, hooke, params):
 164         """Setup `params` from config and return the z piezo and
 165         deflection arrays.
 166         """
 167         curve = self._curve(hooke=hooke, params=params)
 168         for key,value in params.items():
 169             if value == None and key in self.plugin.config:
 170                 # Use configured default value.
 171                 params[key] = self.plugin.config[key]
 172         # TODO: convert 'see double' from nm to points
 173         name,def_unit = split_data_label(params['deflection column'])
 174         params['output peak column'] = join_data_label(
 175             params['output peak column'], def_unit)
 176         return params
 177
 178     def _peak_name(self, params, index):
 179         d_name,d_unit = split_data_label(params['deflection column'])
 180         return 'flat filter peak %d of %s' % (index, d_name)
 181
 182
 183 class FlatFilterCommand (FilterCommand):
 184     u"""Create a subset playlist of curves with enough flat peaks.
 185
 186     Notes
 187     -----
 188     This type of filter is of course very raw, and requires relatively
 189     conservative settings to safely avoid false negatives (that is, to
 190     avoid discarding interesting curves).  Using it on the protein
 191     unfolding experiments described by Sandal [#sandal2008] it has
 192     been found to reduce the data set to analyze by hand by 60-80%.
 193
 194     .. [#sandal2008] M. Sandal, F. Valle, I. Tessari, S. Mammi, E. Bergantino,
 195       F. Musiani, M. Brucale, L. Bubacco, B. Samorì.
 196       "Conformational equilibria in monomeric α-Synuclein at the
 197       single molecule level."
 198       PLOS Biology, 2009.
 199       doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
 200
 201     See Also
 202     --------
 203     FlatPeaksCommand : Underlying flat-based peak detection.
 204     """
 205     def __init__(self, plugin):
 206         flat_peaks = [c for c in plugin._commands
 207                       if c.name == 'flat filter peaks'][0]
 208         flat_peaks_arguments = [a for a in flat_peaks.arguments
 209                                 if a.name not in ['help', 'stack']]
 210         flat_peaks_arg_names = [a.name for a in flat_peaks_arguments]
 211         plugin_arguments = [a for a in plugin._arguments
 212                             if a.name not in flat_peaks_arg_names]
 213         arguments = flat_peaks_arguments + plugin_arguments
 214         super(FlatFilterCommand, self).__init__(
 215             plugin, name='flat filter playlist')
 216         self.arguments.extend(arguments)
 217
 218     def filter(self, curve, hooke, inqueue, outqueue, params):
 219         params['curve'] = curve
 220         inq = Queue()
 221         outq = Queue()
 222         filt_command = hooke.command_by_name['flat filter peaks']
 223         filt_command.run(hooke, inq, outq, **params)
 224         peaks = outq.get()
 225         if isinstance(peaks, UncaughtException) \
 226                 and isinstance(peaks.exception, PoorFit):
 227             return False
 228         if not (isinstance(peaks, list) and (len(peaks) == 0
 229                                              or isinstance(peaks[0], Peak))):
 230             raise Failure('Expected a list of Peaks, not %s: %s'
 231                           % (type(peaks), peaks))
 232         ret = outq.get()
 233         if not isinstance(ret, Success):
 234             raise ret
 235         if params['min peaks'] == None: # Use configured default value.
 236             params['min peaks'] = self.plugin.config['min peaks']
 237         return len(peaks) >= int(params['min peaks'])