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