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 ..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
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)]
88 # append FlatFilterCommand so it can steal arguments from
90 self._commands.append(FlatFilterCommand(self))
92 def dependencies(self):
95 def default_settings(self):
99 class FlatPeaksCommand (ColumnAddingCommand):
100 """Detect peaks in velocity clamp data using noise statistics.
104 Noise analysis on the retraction curve:
106 1) A median window filter (using
107 :func:`scipy.signal.signaltools.medfilt`) smooths the
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`.
114 The algorithm was originally Francesco Musiani's idea.
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',
124 ('distance column', 'surface distance (m)', """
125 Name of the column to use as the surface position input.
127 ('deflection column', 'surface deflection (m)', """
128 Name of the column to use as the deflection input.
132 ('output peak column', 'flat filter peaks', """
133 Name of the column (without units) to use as the peak output.
137 Argument(name='peak info name', type='string',
138 default='flat filter peaks',
140 Name for storing the list of peaks in the `.info` dictionary.
142 ] + plugin_arguments,
143 help=self.__doc__, plugin=plugin)
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'])
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
163 self._set_column(hooke=hooke, params=params,
164 column_name='output peak column',
165 values=peaks_to_mask(def_data, peaks) * def_data)
168 def __setup_params(self, hooke, params):
169 """Setup `params` from config and return the z piezo and
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)
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)
191 class FlatFilterCommand (FilterCommand):
192 u"""Create a subset playlist of curves with enough flat peaks.
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%.
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."
207 doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
211 FlatPeaksCommand : Underlying flat-based peak detection.
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)
226 def filter(self, curve, hooke, inqueue, outqueue, params):
227 params['curve'] = curve
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)
234 if isinstance(peaks, UncaughtException) \
235 and isinstance(peaks.exception, PoorFit):
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))
242 if not isinstance(ret, Success):
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'])