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