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 diff
39 from scipy.signal.signaltools import medfilt
41 from ..command import Command, 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 CurveArgument
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._arguments.extend([ # Non-configurable arguments
88 Argument(name='retract block', type='string',
91 Name of the retracting data block.
93 Argument(name='input distance column', type='string',
94 default='surface distance (m)',
96 Name of the column to use as the surface position input.
98 Argument(name='input deflection column', type='string',
99 default='deflection (N)',
101 Name of the column to use as the deflection input.
103 Argument(name='output peak column', type='string',
104 default='peak deflection',
106 Name of the column (without units) to use as the peak-maske deflection output.
108 Argument(name='peak info name', type='string',
109 default='flat filter peaks',
111 Name for storing the list of peaks in the `.info` dictionary.
114 self._commands = [FlatPeaksCommand(self), FlatFilterCommand(self)]
116 def dependencies(self):
119 def default_settings(self):
120 return self._settings
123 class FlatPeaksCommand (Command):
124 """Detect peaks in velocity clamp data using noise statistics.
128 Noise analysis on the retraction curve:
130 1) A median window filter (using
131 :func:`scipy.signal.signaltools.medfilt`) smooths the
133 2) The deflection derivative is calculated (using
134 :func:`numpy.diff` which uses forward differencing).
135 3) Peaks in the derivative curve are extracted with
136 :func:`~hooke.plugins.peak.find_peaks`.
138 The algorithm was originally Francesco Musiani's idea.
140 def __init__(self, plugin):
141 plugin_arguments = [a for a in plugin._arguments
142 if a.name != 'min peaks']
143 # Drop min peaks, since we're not filtering with this
144 # function, just detecting peaks.
145 super(FlatPeaksCommand, self).__init__(
146 name='flat filter peaks',
149 ] + plugin_arguments,
150 help=self.__doc__, plugin=plugin)
152 def _run(self, hooke, inqueue, outqueue, params):
153 data,block_index,z_data,d_data,params = self._setup(hooke, params)
155 while (start_index < len(z_data)
156 and z_data[start_index] < params['blind window']):
158 median = medfilt(d_data[start_index:], params['median window'])
160 peaks = find_peaks(deriv, **_kwargs(params, find_peaks_arguments,
161 argument_input_keys=True))
162 d_name,d_unit = split_data_label(params['input deflection column'])
163 for i,peak in enumerate(peaks):
164 peak.name = 'flat filter peak %d of %s' % (i, d_name)
165 peak.index += start_index
166 data.info[params['peak info name']] = peaks
168 # Add a peak-masked deflection column.
169 new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
170 new.info = copy.deepcopy(data.info)
172 d_name,d_unit = split_data_label(params['input deflection column'])
173 new.info['columns'].append(
174 join_data_label(params['output peak column'], d_unit))
175 new[:,-1] = peaks_to_mask(d_data, peaks) * d_data
177 params['curve'].data[block_index] = new
179 def _setup(self, hooke, params):
180 """Setup `params` from config and return the z piezo and
183 curve = params['curve']
184 if curve.info['experiment'] != VelocityClamp:
185 raise Failure('%s operates on VelocityClamp experiments, not %s'
186 % (self.name, curve.info['experiment']))
188 for i,block in enumerate(curve.data):
189 if block.info['name'].startswith(params['retract block']):
194 raise Failure('No retraction blocks in %s.' % curve)
195 z_data = data[:,data.info['columns'].index(
196 params['input distance column'])]
197 d_data = data[:,data.info['columns'].index(
198 params['input deflection column'])]
199 for key,value in params.items():
200 if value == None: # Use configured default value.
201 params[key] = self.plugin.config[key]
202 # TODO: convert 'see double' from nm to points
203 return (data, block_index, z_data, d_data, params)
206 class FlatFilterCommand (FilterCommand):
207 u"""Create a subset playlist of curves with enough flat peaks.
211 This type of filter is of course very raw, and requires relatively
212 conservative settings to safely avoid false negatives (that is, to
213 avoid discarding interesting curves). Using it on the protein
214 unfolding experiments described by Sandal [#sandal2008] it has
215 been found to reduce the data set to analyze by hand by 60-80%.
217 .. [#sandal2008] M. Sandal, F. Valle, I. Tessari, S. Mammi, E. Bergantino,
218 F. Musiani, M. Brucale, L. Bubacco, B. Samorì.
219 "Conformational equilibria in monomeric α-Synuclein at the
220 single molecule level."
222 doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
226 FlatPeaksCommand : Underlying flat-based peak detection.
228 def __init__(self, plugin):
229 super(FlatFilterCommand, self).__init__(
230 plugin, name='flat filter playlist')
231 self.arguments.extend(plugin._arguments)
233 def filter(self, curve, hooke, inqueue, outqueue, params):
234 params['curve'] = curve
237 filt_command = [c for c in hooke.commands
238 if c.name=='flat filter peaks'][0]
239 filt_command.run(hooke, inq, outq, **params)
241 if isinstance(peaks, UncaughtException) \
242 and isinstance(peaks.exception, PoorFit):
244 if not (isinstance(peaks, list) and (len(peaks) == 0
245 or isinstance(peaks[0], Peak))):
246 raise Failure('Expected a list of Peaks, not %s' % peaks)
248 if not isinstance(ret, Success):
250 if params['min peaks'] == None: # Use configured default value.
251 params['min peaks'] = self.plugin.config['min peaks']
252 return len(peaks) >= int(params['min peaks'])