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