1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2008-2010 Alberto Gomez-Casado
5 # Massimo Sandal <devicerandom@gmail.com>
6 # W. Trevor King <wking@drexel.edu>
8 # This file is part of Hooke.
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.
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.
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/>.
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.
31 :mod:`~hooke.plugin.convfilt` for a convolution-based filter for
32 :class:`~hooke.experiment.VelocityClamp` experiments.
36 from Queue import Queue
38 from numpy import absolute, diff
39 from scipy.signal.signaltools import medfilt
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 ..plugin import Plugin, argument_to_setting
46 from ..plugin.curve import ColumnAddingCommand
47 from ..plugin.playlist import FilterCommand
48 from ..util.fit import PoorFit
49 from ..util.peak import (find_peaks, peaks_to_mask,
50 find_peaks_arguments, Peak, _kwargs)
51 from ..util.si import join_data_label, split_data_label
54 class FlatFiltPlugin (Plugin):
55 """Standard-devitiation-based peak recognition and filtering.
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).
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.
67 Argument('min peaks', type='int', default=4, help="""
68 Minimum number of peaks for curve acceptance.
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'),
75 ('min deviations', 9.0),
77 ('see double', 10), # TODO: points vs. meters. 10e-9),
79 argument = [a for a in self._arguments if a.name == key][0]
80 argument.default = value
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), FlatFilterCommand(self)]
89 def dependencies(self):
92 def default_settings(self):
96 class FlatPeaksCommand (ColumnAddingCommand):
97 """Detect peaks in velocity clamp data using noise statistics.
101 Noise analysis on the retraction curve:
103 1) A median window filter (using
104 :func:`scipy.signal.signaltools.medfilt`) smooths the
106 2) The deflection derivative is calculated (using
107 :func:`numpy.diff` which uses forward differencing).
108 3) Peaks in the derivative curve are extracted with
109 :func:`~hooke.plugins.peak.find_peaks`.
111 The algorithm was originally Francesco Musiani's idea.
113 def __init__(self, plugin):
114 plugin_arguments = [a for a in plugin._arguments
115 if a.name != 'min peaks']
116 # Drop min peaks, since we're not filtering with this
117 # function, just detecting peaks.
118 super(FlatPeaksCommand, self).__init__(
119 name='flat filter peaks',
121 ('distance column', 'surface distance (m)', """
122 Name of the column to use as the surface position input.
124 ('deflection column', 'surface deflection (m)', """
125 Name of the column to use as the deflection input.
129 ('output peak column', 'flat filter peaks', """
130 Name of the column (without units) to use as the peak output.
134 Argument(name='peak info name', type='string',
135 default='flat filter peaks',
137 Name (without units) for storing the list of peaks in the `.info`
140 ] + plugin_arguments,
141 help=self.__doc__, plugin=plugin)
143 def _run(self, hooke, inqueue, outqueue, params):
144 params = self.__setup_params(hooke=hooke, params=params)
145 block = self._block(hooke=hooke, params=params)
146 dist_data = self._get_column(hooke=hooke, params=params,
147 column_name='distance column')
148 def_data = self._get_column(hooke=hooke, params=params,
149 column_name='deflection column')
150 start_index = absolute(dist_data-params['blind window']).argmin()
151 median = medfilt(def_data[start_index:], params['median window'])
153 peaks = find_peaks(deriv, **_kwargs(params, find_peaks_arguments,
154 argument_input_keys=True))
155 for i,peak in enumerate(peaks):
156 peak.name = self._peak_name(params, i)
157 peak.index += start_index
158 block.info[params['peak info name']] = peaks
160 self._set_column(hooke=hooke, params=params,
161 column_name='output peak column',
162 values=peaks_to_mask(def_data, peaks) * def_data)
165 def __setup_params(self, hooke, params):
166 """Setup `params` from config and return the z piezo and
169 curve = self._curve(hooke=hooke, params=params)
170 if curve.info['experiment'] != VelocityClamp:
171 raise Failure('%s operates on VelocityClamp experiments, not %s'
172 % (self.name, curve.info['experiment']))
173 for key,value in params.items():
174 if value == None: # Use configured default value.
175 params[key] = self.plugin.config[key]
176 # TODO: convert 'see double' from nm to points
177 name,def_unit = split_data_label(params['deflection column'])
178 params['output peak column'] = join_data_label(
179 params['output peak column'], def_unit)
180 params['peak info name'] = join_data_label(
181 params['peak info name'], def_unit)
184 def _peak_name(self, params, index):
185 d_name,d_unit = split_data_label(params['deflection column'])
186 return 'flat filter peak %d of %s' % (index, d_name)
189 class FlatFilterCommand (FilterCommand):
190 u"""Create a subset playlist of curves with enough flat peaks.
194 This type of filter is of course very raw, and requires relatively
195 conservative settings to safely avoid false negatives (that is, to
196 avoid discarding interesting curves). Using it on the protein
197 unfolding experiments described by Sandal [#sandal2008] it has
198 been found to reduce the data set to analyze by hand by 60-80%.
200 .. [#sandal2008] M. Sandal, F. Valle, I. Tessari, S. Mammi, E. Bergantino,
201 F. Musiani, M. Brucale, L. Bubacco, B. Samorì.
202 "Conformational equilibria in monomeric α-Synuclein at the
203 single molecule level."
205 doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
209 FlatPeaksCommand : Underlying flat-based peak detection.
211 def __init__(self, plugin):
212 super(FlatFilterCommand, self).__init__(
213 plugin, name='flat filter playlist')
214 self.arguments.extend(plugin._arguments) # TODO: curve & block arguments
216 def filter(self, curve, hooke, inqueue, outqueue, params):
217 params['curve'] = curve
220 filt_command = [c for c in hooke.commands
221 if c.name=='flat filter peaks'][0]
222 filt_command.run(hooke, inq, outq, **params)
224 if isinstance(peaks, UncaughtException) \
225 and isinstance(peaks.exception, PoorFit):
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' % peaks)
231 if not isinstance(ret, Success):
233 if params['min peaks'] == None: # Use configured default value.
234 params['min peaks'] = self.plugin.config['min peaks']
235 return len(peaks) >= int(params['min peaks'])