1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2010-2012 W. Trevor King <wking@tremily.us>
5 # This file is part of Hooke.
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
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
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/>.
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.
27 :mod:`~hooke.plugin.convfilt` for a convolution-based filter.
31 from Queue import Queue
33 from numpy import absolute, diff
34 from scipy.signal.signaltools import medfilt
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
48 class FlatFiltPlugin (Plugin):
49 """Standard-devitiation-based peak recognition and filtering.
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).
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.
61 Argument('min peaks', type='int', default=4, help="""
62 Minimum number of peaks for curve acceptance.
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'),
69 ('min deviations', 9.0),
71 ('see double', 10), # TODO: points vs. meters. 10e-9),
73 argument = [a for a in self._arguments if a.name == key][0]
74 argument.default = value
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
84 self._commands.append(FlatFilterCommand(self))
86 def dependencies(self):
89 def default_settings(self):
93 class FlatPeaksCommand (ColumnAddingCommand):
94 """Detect peaks in velocity clamp data using noise statistics.
98 Noise analysis on the retraction curve:
100 1) A median window filter (using
101 :func:`scipy.signal.signaltools.medfilt`) smooths the
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`.
108 The algorithm was originally Francesco Musiani's idea.
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',
118 ('distance column', 'surface distance (m)', """
119 Name of the column to use as the surface position input.
121 ('deflection column', 'surface deflection (m)', """
122 Name of the column to use as the deflection input.
126 ('output peak column', 'flat filter peaks', """
127 Name of the column (without units) to use as the peak output.
131 Argument(name='peak info name', type='string',
132 default='flat filter peaks',
134 Name for storing the list of peaks in the `.info` dictionary.
136 ] + plugin_arguments,
137 help=self.__doc__, plugin=plugin)
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'])
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
157 self._set_column(hooke=hooke, params=params,
158 column_name='output peak column',
159 values=peaks_to_mask(def_data, peaks) * def_data)
162 def _setup_params(self, hooke, params):
163 """Setup `params` from config and return the z piezo and
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)
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)
182 class FlatFilterCommand (FilterCommand):
183 u"""Create a subset playlist of curves with enough flat peaks.
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%.
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."
198 doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
202 FlatPeaksCommand : Underlying flat-based peak detection.
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)
217 def filter(self, curve, hooke, inqueue, outqueue, params):
218 params['curve'] = curve
221 filt_command = hooke.command_by_name['flat filter peaks']
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: %s'
230 % (type(peaks), peaks))
232 if not isinstance(ret, Success):
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'])