-# Copyright (C) 2010 Fibrin's Benedetti
-# W. Trevor King <wking@drexel.edu>
+# Copyright (C) 2008-2010 Alberto Gomez-Casado
+# Fabrizio Benedetti
+# Massimo Sandal <devicerandom@gmail.com>
+# W. Trevor King <wking@drexel.edu>
#
# This file is part of Hooke.
#
"""
from ..command import Command, Argument, Failure
+from ..curve import Data
from ..plugin import Builtin
from ..plugin.playlist import current_playlist_callback
+from ..util.calculus import derivative
class CurvePlugin (Builtin):
super(CurvePlugin, self).__init__(name='curve')
def commands(self):
- return [InfoCommand(), ]
+ return [InfoCommand(), ExportCommand()]
# Define common or complicated arguments
def _get_block_sizes(self, curve):
return [block.shape for block in curve.data]
-
class ExportCommand (Command):
"""Export a :class:`hooke.curve.Curve` data block as TAB-delimeted
ASCII text.
"""
def __init__(self):
- super(InfoCommand, self).__init__(
- name='curve info',
+ super(ExportCommand, self).__init__(
+ name='export block',
arguments=[
CurveArgument,
Argument(name='block', aliases=['set'], type='int', default=0,
- help="""
+ help="""
Data block to save. For an approach/retract force curve, `0` selects
the approacing curve and `1` selects the retracting curve.
""".strip()),
help=self.__doc__)
def _run(self, hooke, inqueue, outqueue, params):
- data = params['curve'].data[params['index']]
+ data = params['curve'].data[params['block']]
f = open(params['output'], 'w')
data.tofile(f, sep='\t')
f.close()
+
+class DifferenceCommand (Command):
+ """Calculate the derivative (actually, the discrete differentiation)
+ of a curve data block.
+
+ See :func:`hooke.util.calculus.derivative` for implementation
+ details.
+ """
+ def __init__(self):
+ super(DifferenceCommand, self).__init__(
+ name='block difference',
+ arguments=[
+ CurveArgument,
+ Argument(name='block one', aliases=['set one'], type='int',
+ default=1,
+ help="""
+Block A in A-B. For an approach/retract force curve, `0` selects the
+approacing curve and `1` selects the retracting curve.
+""".strip()),
+ Argument(name='block two', aliases=['set two'], type='int',
+ default=0,
+ help='Block B in A-B.'),
+ Argument(name='x column', type='int', default=0,
+ help="""
+Column of data block to differentiate with respect to.
+""".strip()),
+ Argument(name='y column', type='int', default=1,
+ help="""
+Column of data block to differentiate.
+""".strip()),
+ ],
+ help=self.__doc__)
+
+ def _run(self, hooke, inqueue, outqueue, params):
+ a = params['curve'].data[params['block one']]
+ b = params['curve'].data[params['block two']]
+ assert a[:,params['x column']] == b[:,params['x column']]:
+ out = Data((a.shape[0],2), dtype=a.dtype)
+ out[:,0] = a[:,params['x column']]
+ out[:,1] = a[:,params['y column']] - b[:,params['y column']]:
+ outqueue.put(out)
+
+class DerivativeCommand (Command):
+ """Calculate the difference between two blocks of data.
+ """
+ def __init__(self):
+ super(DerivativeCommand, self).__init__(
+ name='block derivative',
+ arguments=[
+ CurveArgument,
+ Argument(name='block', aliases=['set'], type='int', default=0,
+ help="""
+Data block to differentiate. For an approach/retract force curve, `0`
+selects the approacing curve and `1` selects the retracting curve.
+""".strip()),
+ Argument(name='x column', type='int', default=0,
+ help="""
+Column of data block to differentiate with respect to.
+""".strip()),
+ Argument(name='y column', type='int', default=1,
+ help="""
+Column of data block to differentiate.
+""".strip()),
+ Argument(name='weights', type='dict', default={-1:-0.5, 1:0.5},
+ help="""
+Weighting scheme dictionary for finite differencing. Defaults to
+central differencing.
+""".strip()),
+ ],
+ help=self.__doc__)
+
+ def _run(self, hooke, inqueue, outqueue, params):
+ data = params['curve'].data[params['block']]
+ outqueue.put(derivative(
+ block, x_col=params['x column'], y_col=params['y column'],
+ weights=params['weights']))
class generalvclampCommands(object):
+ def do_subtplot(self, args):
+ '''
+ SUBTPLOT
+ (procplots.py plugin)
+ Plots the difference between ret and ext current curve
+ -------
+ Syntax: subtplot
+ '''
+ #FIXME: sub_filter and sub_order must be args
+
+ if len(self.plots[0].vectors) != 2:
+ print 'This command only works on a curve with two different plots.'
+ pass
+
+ outplot=self.subtract_curves(sub_order=1)
+
+ plot_graph=self.list_of_events['plot_graph']
+ wx.PostEvent(self.frame,plot_graph(plots=[outplot]))
+
def _plug_init(self):
self.basecurrent=None
self.basepoints=None
--- /dev/null
+
+class Plotmanip (object):
+#-----PLOT MANIPULATORS
+ def plotmanip_median(self, plot, current, customvalue=None):
+ '''
+ does the median of the y values of a plot
+ '''
+ if customvalue:
+ median_filter=customvalue
+ else:
+ median_filter=self.config['medfilt']
+
+ if median_filter==0:
+ return plot
+
+ if float(median_filter)/2 == int(median_filter)/2:
+ median_filter+=1
+
+ nplots=len(plot.vectors)
+ c=0
+ while c<nplots:
+ plot.vectors[c][1]=scipy.signal.medfilt(plot.vectors[c][1],median_filter)
+ c+=1
+
+ return plot
+
+
+ def plotmanip_correct(self, plot, current, customvalue=None):
+ '''
+ does the correction for the deflection for a force spectroscopy curve.
+ Assumes that:
+ - the current plot has a deflection() method that returns a vector of values
+ - the deflection() vector is as long as the X of extension + the X of retraction
+ - plot.vectors[0][0] is the X of extension curve
+ - plot.vectors[1][0] is the X of retraction curve
+
+ FIXME: both this method and the picoforce driver have to be updated, deflection() must return
+ a more senseful data structure!
+ '''
+ #use only for force spectroscopy experiments!
+ if current.curve.experiment != 'smfs':
+ return plot
+
+ if customvalue != None:
+ execute_me=customvalue
+ else:
+ execute_me=self.config['correct']
+ if not execute_me:
+ return plot
+
+ defl_ext,defl_ret=current.curve.deflection()
+ #halflen=len(deflall)/2
+
+ plot.vectors[0][0]=[(zpoint-deflpoint) for zpoint,deflpoint in zip(plot.vectors[0][0],defl_ext)]
+ plot.vectors[1][0]=[(zpoint-deflpoint) for zpoint,deflpoint in zip(plot.vectors[1][0],defl_ret)]
+
+ return plot
+
+
+ def plotmanip_centerzero(self, plot, current, customvalue=None):
+ '''
+ Centers the force curve so the median (the free level) corresponds to 0 N
+ Assumes that:
+ - plot.vectors[0][1] is the Y of extension curve
+ - plot.vectors[1][1] is the Y of retraction curve
+
+
+ '''
+ #use only for force spectroscopy experiments!
+ if current.curve.experiment != 'smfs':
+ return plot
+
+ if customvalue != None:
+ execute_me=customvalue
+ else:
+ execute_me=self.config['centerzero']
+ if not execute_me:
+ return plot
+
+
+
+ #levelapp=float(np.median(plot.vectors[0][1]))
+ levelret=float(np.median(plot.vectors[1][1][-300:-1]))
+
+ level=levelret
+
+ approach=[i-level for i in plot.vectors[0][1]]
+ retract=[i-level for i in plot.vectors[1][1]]
+
+ plot.vectors[0][1]=approach
+ plot.vectors[1][1]=retract
+ return plot
+++ /dev/null
-# Copyright (C) 2008-2010 Alberto Gomez-Casado
-# Massimo Sandal <devicerandom@gmail.com>
-# W. Trevor King <wking@drexel.edu>
-#
-# This file is part of Hooke.
-#
-# Hooke is free software: you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation, either
-# version 3 of the License, or (at your option) any later version.
-#
-# Hooke is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with Hooke. If not, see
-# <http://www.gnu.org/licenses/>.
-
-"""Processed plots plugin for force curves.
-"""
-
-from ..libhooke import WX_GOOD
-import wxversion
-wxversion.select(WX_GOOD)
-
-import wx
-import numpy as np
-import scipy as sp
-import scipy.signal
-import copy
-
-from .. import curve as lhc
-
-
-class procplotsCommands(object):
-
- def _plug_init(self):
- pass
-
- def do_derivplot(self,args):
- '''
- DERIVPLOT
- (procplots.py plugin)
- Plots the derivate (actually, the discrete differentiation) of the current force curve
- ---------
- Syntax: derivplot
- '''
- dplot=self.derivplot_curves()
- plot_graph=self.list_of_events['plot_graph']
- wx.PostEvent(self.frame,plot_graph(plots=[dplot]))
-
- def derivplot_curves(self):
- '''
- do_derivplot helper function
-
- create derivate plot curves for force curves.
- '''
- dplot=lhc.PlotObject()
- dplot.vectors=[]
-
- for vector in self.plots[0].vectors:
- dplot.vectors.append([])
- dplot.vectors[-1].append(vector[0][:-1])
- dplot.vectors[-1].append(np.diff(vector[1]))
-
- dplot.destination=1
- dplot.units=self.plots[0].units
-
- return dplot
-
- def do_subtplot(self, args):
- '''
- SUBTPLOT
- (procplots.py plugin)
- Plots the difference between ret and ext current curve
- -------
- Syntax: subtplot
- '''
- #FIXME: sub_filter and sub_order must be args
-
- if len(self.plots[0].vectors) != 2:
- print 'This command only works on a curve with two different plots.'
- pass
-
- outplot=self.subtract_curves(sub_order=1)
-
- plot_graph=self.list_of_events['plot_graph']
- wx.PostEvent(self.frame,plot_graph(plots=[outplot]))
-
- def subtract_curves(self, sub_order):
- '''
- subtracts the extension from the retraction
- '''
- xext=self.plots[0].vectors[0][0]
- yext=self.plots[0].vectors[0][1]
- xret=self.plots[0].vectors[1][0]
- yret=self.plots[0].vectors[1][1]
-
- #we want the same number of points
- maxpoints_tot=min(len(xext),len(xret))
- xext=xext[0:maxpoints_tot]
- yext=yext[0:maxpoints_tot]
- xret=xret[0:maxpoints_tot]
- yret=yret[0:maxpoints_tot]
-
- if sub_order:
- ydiff=[yretval-yextval for yretval,yextval in zip(yret,yext)]
- else: #reverse subtraction (not sure it's useful, but...)
- ydiff=[yextval-yretval for yextval,yretval in zip(yext,yret)]
-
- outplot=copy.deepcopy(self.plots[0])
- outplot.vectors[0][0], outplot.vectors[1][0] = xext,xret #FIXME: if I use xret, it is not correct!
- outplot.vectors[1][1]=ydiff
- outplot.vectors[0][1]=[0 for item in outplot.vectors[1][0]]
-
- return outplot
-
-
-#-----PLOT MANIPULATORS
- def plotmanip_median(self, plot, current, customvalue=None):
- '''
- does the median of the y values of a plot
- '''
- if customvalue:
- median_filter=customvalue
- else:
- median_filter=self.config['medfilt']
-
- if median_filter==0:
- return plot
-
- if float(median_filter)/2 == int(median_filter)/2:
- median_filter+=1
-
- nplots=len(plot.vectors)
- c=0
- while c<nplots:
- plot.vectors[c][1]=scipy.signal.medfilt(plot.vectors[c][1],median_filter)
- c+=1
-
- return plot
-
-
- def plotmanip_correct(self, plot, current, customvalue=None):
- '''
- does the correction for the deflection for a force spectroscopy curve.
- Assumes that:
- - the current plot has a deflection() method that returns a vector of values
- - the deflection() vector is as long as the X of extension + the X of retraction
- - plot.vectors[0][0] is the X of extension curve
- - plot.vectors[1][0] is the X of retraction curve
-
- FIXME: both this method and the picoforce driver have to be updated, deflection() must return
- a more senseful data structure!
- '''
- #use only for force spectroscopy experiments!
- if current.curve.experiment != 'smfs':
- return plot
-
- if customvalue != None:
- execute_me=customvalue
- else:
- execute_me=self.config['correct']
- if not execute_me:
- return plot
-
- defl_ext,defl_ret=current.curve.deflection()
- #halflen=len(deflall)/2
-
- plot.vectors[0][0]=[(zpoint-deflpoint) for zpoint,deflpoint in zip(plot.vectors[0][0],defl_ext)]
- plot.vectors[1][0]=[(zpoint-deflpoint) for zpoint,deflpoint in zip(plot.vectors[1][0],defl_ret)]
-
- return plot
-
-
- def plotmanip_centerzero(self, plot, current, customvalue=None):
- '''
- Centers the force curve so the median (the free level) corresponds to 0 N
- Assumes that:
- - plot.vectors[0][1] is the Y of extension curve
- - plot.vectors[1][1] is the Y of retraction curve
-
-
- '''
- #use only for force spectroscopy experiments!
- if current.curve.experiment != 'smfs':
- return plot
-
- if customvalue != None:
- execute_me=customvalue
- else:
- execute_me=self.config['centerzero']
- if not execute_me:
- return plot
-
-
-
- #levelapp=float(np.median(plot.vectors[0][1]))
- levelret=float(np.median(plot.vectors[1][1][-300:-1]))
-
- level=levelret
-
- approach=[i-level for i in plot.vectors[0][1]]
- retract=[i-level for i in plot.vectors[1][1]]
-
- plot.vectors[0][1]=approach
- plot.vectors[1][1]=retract
- return plot
-
- '''
- def plotmanip_detriggerize(self, plot, current, customvalue=None):
- #DEPRECATED
- if self.config['detrigger']==0:
- return plot
-
- cutindex=2
- startvalue=plot.vectors[0][0][0]
-
- for index in range(len(plot.vectors[0][0])-1,2,-2):
- if plot.vectors[0][0][index]>startvalue:
- cutindex=index
- else:
- break
-
- plot.vectors[0][0]=plot.vectors[0][0][:cutindex]
- plot.vectors[0][1]=plot.vectors[0][1][:cutindex]
-
- return plot
- '''
-
-
-
-#FFT---------------------------
- def fft_plot(self, vector):
- '''
- calculates the fast Fourier transform for the selected vector in the plot
- '''
- fftplot=lhc.PlotObject()
- fftplot.vectors=[[]]
-
- fftlen=len(vector)/2 #need just 1/2 of length
- fftplot.vectors[-1].append(np.arange(1,fftlen).tolist())
-
- try:
- fftplot.vectors[-1].append(abs(np.fft(vector)[1:fftlen]).tolist())
- except TypeError: #we take care of newer NumPy (1.0.x)
- fftplot.vectors[-1].append(abs(np.fft.fft(vector)[1:fftlen]).tolist())
-
-
- fftplot.destination=1
-
-
- return fftplot
-
-
- def do_fft(self,args):
- '''
- FFT
- (procplots.py plugin)
- Plots the fast Fourier transform of the selected plot
- ---
- Syntax: fft [top,bottom] [select] [0,1...]
-
- By default, fft performs the Fourier transform on all the 0-th data set on the
- top plot.
-
- [top,bottom]: which plot is the data set to fft (default: top)
- [select]: you pick up two points on the plot and fft only the segment between
- [0,1,...]: which data set on the selected plot you want to fft (default: 0)
- '''
-
- #args parsing
- #whatplot = plot to fft
- #whatset = set to fft in the plot
- select=('select' in args)
- if 'top' in args:
- whatplot=0
- elif 'bottom' in args:
- whatplot=1
- else:
- whatplot=0
- whatset=0
- for arg in args:
- try:
- whatset=int(arg)
- except ValueError:
- pass
-
- if select:
- points=self._measure_N_points(N=2, whatset=whatset)
- boundaries=[points[0].index, points[1].index]
- boundaries.sort()
- y_to_fft=self.plots[whatplot].vectors[whatset][1][boundaries[0]:boundaries[1]] #y
- else:
- y_to_fft=self.plots[whatplot].vectors[whatset][1] #y
-
- fftplot=self.fft_plot(y_to_fft)
- fftplot.units=['frequency', 'power']
- plot_graph=self.list_of_events['plot_graph']
- wx.PostEvent(self.frame,plot_graph(plots=[fftplot]))
--- /dev/null
+# Copyright
+
+"""The `calculus` module provides functions for calculating
+derivatives and integrals of discrete data.
+"""
+
+import copy
+
+import numpy
+
+from ..curve import Data
+
+
+def derivative(data, x_col=0, f_col=1, weights={-1:-0.5, 1:0.5}):
+ """Calculate the discrete derivative (finite difference) of
+ data[:,f_col] with respect to data[:,x_col].
+
+ Examples
+ --------
+
+ >>> import pprint
+ >>> d = Data((5,2), dtype=numpy.float,
+ ... info={'columns':['x', 'x**2']})
+ >>> for i in range(5):
+ ... d[i,0] = i
+ ... d[i,1] = i**2
+ >>> d
+ Data([[ 0., 0.],
+ [ 1., 1.],
+ [ 2., 4.],
+ [ 3., 9.],
+ [ 4., 16.]])
+ >>> dd = derivative(d)
+ >>> dd
+ Data([[ 0., 1.],
+ [ 1., 2.],
+ [ 2., 4.],
+ [ 3., 6.],
+ [ 4., 7.]])
+ >>> pprint.pprint(dd.info)
+ {'columns': ['x', 'deriv x**2 with respect to x']}
+
+ Notes
+ -----
+
+ Weights
+ ~~~~~~~
+
+ The returned :class:`Data` block shares its x vector with the
+ input data. The ith df/dx value in the returned data is
+ caclulated with::
+
+ (df/dx)[i] = (SUM_j w[j] f[i+j]) / h
+
+ where ``h = x[i+1]-x[i]`` is the x coordinate spacing (assumed
+ constant) and ``j`` ranges over the keys of `weights`.
+
+ There standard schemes translate as follows:
+
+ ======== ====================== ===================
+ scheme formula weights
+ ======== ====================== ===================
+ forward ``(f[i+1]-f[i])/h`` ``{0:-1,1:1}``
+ backward ``(f[i]-f[i-1])/h`` ``{0:1,-1:-1}``
+ central ``(f[i+1]-f[i-1])/2h`` ``{-1:-0.5,1:0.5}``
+ ======== ====================== ===================
+
+ The default scheme is central differencing.
+
+ Boundary conditions
+ ~~~~~~~~~~~~~~~~~~~
+
+ These could be configurable in principle. The current scheme just
+ extrapolates virtual points out to negative `i` following::
+
+ f[i<0] = 2*f[0] - f[-i]
+
+ With analogous treatment for `i > data.shape[0]`. This ensures that
+ `f[i]-f[0]` is odd about `i=0`, which keeps derivatives smooth.::
+
+ f[i] - f[0] = f[0] - f[-i] == -(f[-i] - f[0])
+ """
+ output = Data((data.shape[0],2), dtype=data.dtype)
+ output.info = copy.copy(data.info)
+ output.info['columns'] = [
+ data.info['columns'][x_col],
+ 'deriv %s with respect to %s' \
+ % (data.info['columns'][f_col], data.info['columns'][x_col]),
+ ]
+ h = data[1,x_col] - data[0,x_col]
+ chunks = []
+ for i,w in weights.items():
+ chunk = numpy.roll(w*data[:,f_col], -i)
+ if i > 0: # chunk shifted down, replace the high `i`s
+ zero = len(chunk) - 1 - i
+ for j in range(1,i+1):
+ chunk[zero+j] = 2*chunk[zero] - chunk[zero-j]
+ elif i < 0: # chunk shifted up, replace the low `i`s
+ zero = -i
+ for j in range(1,zero+1):
+ chunk[zero-j] = 2*chunk[zero] - chunk[zero+j]
+ chunks.append(chunk)
+ print chunks
+ output[:,0] = data[:,x_col]
+ output[:,1] = sum(chunks)
+ return output