1 # Copyright (C) 2010 Alberto Gomez-Casado
2 # W. Trevor King <wking@drexel.edu>
4 # This file is part of Hooke.
6 # Hooke is free software: you can redistribute it and/or modify it
7 # under the terms of the GNU Lesser General Public License as
8 # published by the Free Software Foundation, either version 3 of the
9 # License, or (at your option) any later version.
11 # Hooke is distributed in the hope that it will be useful, but WITHOUT
12 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
13 # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
14 # Public License for more details.
16 # You should have received a copy of the GNU Lesser General Public
17 # License along with Hooke. If not, see
18 # <http://www.gnu.org/licenses/>.
20 """The ``multifit`` module provides :class:`MultifitPlugin` and
21 several associated :class:`hooke.command.Command`\s for fitting
22 all sorts of :class:`hooke.experiment.VelocityClamp` parameters.
25 from ..command import Command, Argument, Failure
26 from ..playlist import FilePlaylist
28 from .playlist import PlaylistArgument
31 class MultifitPlugin (Plugin):
33 super(MultifitPlugin, self).__init__(name='multifit')
34 self._commands = [MultifitCommand(self)]
37 class MultifitCommand (Command):
38 """Presents curves for interactive manual analysis.
40 Obtains contour length, persistance length, rupture force and
43 See :class:`hooke.plugin.fit.FitCommand` help for more information
44 on the options and fit procedures.
46 #NOTE: centerzero plot modifier should be activated (set centerzero
48 def __init__(self, plugin):
49 super(MultifitCommand, self).__init__(
53 Argument(name='persistence length', type='float',
55 Persistence length in meters for WLC fits. If not given, the fit will
58 Argument(name='Kuhn length', type='float',
60 Kuhn length in meters for eFJC fits. If not given, the fit will be a
61 2-variable fit. The eFJC fit works better with a fixed Kuhn length
64 Argument(name='temperature', type='float', default=293.0,
66 Temperature in Kelvin. Defaults to 293 K.
68 Argument(name='slope width', type='int',
70 Width in points for slope fitting (points leading up to the rupture).
71 If not given, the routine will prompt for an extra point to fit the
74 Argument(name='base width', type='int', default=15,
76 Width in points for slope fitting (points following the rupture?).
78 Argument(name='point slope', type='bool', default=False,
80 Calculate the slope from the two selected points instead of fitting
83 Argument(name='just one', type='bool', default=False,
85 Perform the fits on the current curve instead of iterating.
88 help=self.__doc__, plugin=plugin)
90 def _run(self, hooke, inqueue, outqueue, params):
95 start_curve = params['playlist'].current()
98 outqueue.put('TODO to end %s execution' % self.name)
99 self._fit_curve(hooke, inqueue, outqueue, params, curve)
100 if params['just one'] == True:
102 params['playlist'].next()
103 next_curve = params['playlist'].current()
104 if next_curve == start_curve:
105 break # loop through playlist complete
108 def _fit_curve(self, hooke, inqueue, outqueue, params, curve):
109 contact = self._get_point('Mark the contact point of the curve.')
110 #hightlights an area dedicated to reject current curve
111 sizeret=len(self.plots[0].vectors[1][0])
112 noarea=self.plots[0].vectors[1][0][-sizeret/6:-1]
113 noareaplot=copy.deepcopy(self._get_displayed_plot())
114 #noise dependent slight shift to improve visibility
115 noareaplot.add_set(noarea,np.ones_like(noarea)*5.0*np.std(self.plots[0].vectors[1][1][-sizeret/6:-1]))
116 noareaplot.styles+=['scatter']
117 noareaplot.colors+=["red"]
118 self._send_plot([noareaplot])
119 contact_point=self._measure_N_points(N=1, whatset=1)[0]
120 contact_point_index=contact_point.index
121 noareaplot.remove_set(-1)
122 #remove set leaves some styles disturbing the next graph, fixed by doing do_plot
124 if contact_point_index>5.0/6.0*sizeret:
125 if params['just one']:
130 self.wlccontact_point=contact_point
131 self.wlccontact_index=contact_point.index
132 self.wlccurrent=self.current.path
134 print 'Now click two points for the fitting area (one should be the rupture point)'
135 wlcpoints=self._measure_N_points(N=2,whatset=1)
136 print 'And one point of the top of the jump'
137 toppoint=self._measure_N_points(N=1,whatset=1)
138 slopepoint=ClickedPoint()
140 print 'Click a point to calculate slope'
141 slopepoint=self._measure_N_points(N=1,whatset=1)
143 fitpoints=[contact_point]+wlcpoints
144 #use the currently displayed plot for the fit
145 displayed_plot=self._get_displayed_plot()
147 #use both fit functions
149 wlcparams, wlcyfit, wlcxfit, wlcfit_errors,wlc_qstd = self.wlc_fit(fitpoints, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T, return_errors=True )
152 print 'WLC fit not possible'
156 fjcparams, fjcyfit, fjcxfit, fjcfit_errors,fjc_qstd = self.efjc_fit(fitpoints, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],kl_value,T, return_errors=True )
159 print 'eFJC fit not possible'
162 #Measure rupture force
163 ruptpoint=ClickedPoint()
164 if wlcpoints[0].graph_coords[1]<wlcpoints[1].graph_coords[1]:
165 ruptpoint=wlcpoints[0]
167 ruptpoint=wlcpoints[1]
168 tpindex=toppoint[0].index
169 toplevel=np.average(displayed_plot.vectors[1][1][tpindex:tpindex+basew])
170 force=toplevel-ruptpoint.graph_coords[1]
172 #Measure the slope - loading rate
176 xint=ruptpoint.graph_coords[0]-slopepoint[0].graph_coords[0]
177 yint=ruptpoint.graph_coords[1]-slopepoint[0].graph_coords[1]
179 slopeplot=self._get_displayed_plot(0) #get topmost displayed plot
180 slopeplot.add_set([ruptpoint.graph_coords[0],slopepoint[0].graph_coords[0]],[ruptpoint.graph_coords[1],slopepoint[0].graph_coords[1]])
181 if slopeplot.styles==[]:
182 slopeplot.styles=[None,None,'scatter']
183 slopeplot.colors=[None,None,'red']
185 slopeplot.styles+=['scatter']
186 slopeplot.colors+=['red']
188 slope=self._slope([ruptpoint]+slopepoint,slopew)
190 slope=self._slope([ruptpoint],slopew)
192 #plot results (_slope already did)
194 #now we have the fit, we can plot it
195 #add the clicked points in the final PlotObject
196 clickvector_x, clickvector_y=[], []
197 for item in wlcpoints:
198 clickvector_x.append(item.graph_coords[0])
199 clickvector_y.append(item.graph_coords[1])
201 #create a custom PlotObject to gracefully plot the fit along the curves
203 #first fix those irritating zoom-destroying fits
204 lowestpoint=min(displayed_plot.vectors[1][1])
205 for i in np.arange(0,len(wlcyfit)):
206 if wlcyfit[i] < 1.2*lowestpoint:
208 for i in np.arange(0,len(fjcyfit)):
209 if fjcyfit[i] < 1.2*lowestpoint:
212 fitplot=copy.deepcopy(displayed_plot)
213 fitplot.add_set(wlcxfit,wlcyfit)
214 fitplot.add_set(fjcxfit,fjcyfit)
215 fitplot.add_set(clickvector_x,clickvector_y)
217 fitplot.styles+=[None,None,'scatter']
218 fitplot.colors+=["red","blue",None]
220 self._send_plot([fitplot])
223 #present results of measurement
224 if len(wlcparams)==1:
225 wlcparams.append(pl_value*1e-9)
226 if len(fjcparams)==1:
227 fjcparams.append(kl_value*1e-9)
230 fjcfit_nm=[i*(10**9) for i in fjcfit_errors]
231 if len(fjcfit_nm)==1:
236 wlcfit_nm=[i*(10**9) for i in wlcfit_errors]
237 if len(wlcfit_nm)==1:
242 wlc_fitq=wlc_qstd/np.std(displayed_plot.vectors[1][1][-20:-1])
243 fjc_fitq=fjc_qstd/np.std(displayed_plot.vectors[1][1][-20:-1])
246 print 'WLC contour : '+str(1e9*wlcparams[0])+u' \u00b1 '+str(wlcfit_nm[0])+' nm'
247 print 'Per. length : '+str(1e9*wlcparams[1])+u' \u00b1 '+str(wlcfit_nm[1])+' nm'
248 print 'Quality :'+str(wlc_fitq)
250 print 'eFJC contour : '+str(1e9*fjcparams[0])+u' \u00b1 '+str(fjcfit_nm[0])+' nm'
251 print 'Kuhn length (e): '+str(1e9*fjcparams[1])+u' \u00b1 '+str(fjcfit_nm[1])+' nm'
252 print 'Quality :'+str(fjc_fitq)
254 print 'Force : '+str(1e12*force)+' pN'
255 print 'Slope : '+str(slope)+' N/m'
258 #FIXME all drivers should provide retract velocity, in SI or homogeneous units
259 lrate=slope*self.current.curve.retract_velocity
260 print 'Loading rate (UNTESTED):'+str(lrate)
262 print 'Loading rate : n/a'
265 #save accept/discard/repeat
266 print '\n\nDo you want to save these?'
270 if self.autofile=='':
271 self.autofile=raw_input('Results filename? (press enter for a random generated name)')
272 if self.autofile=='':
273 [osf,name]=tempfile.mkstemp(prefix='analysis-',suffix='.txt',text=True,dir=self.config['hookedir'])
274 print 'saving in '+name
277 f=open(self.autofile,'a+')
278 f.write('Analysis started '+time.asctime()+'\n')
279 f.write('----------------------------------------\n')
280 f.write(' File ; WLC Cont. length (nm) ; Sigma WLC cl; Per. Length (nm) ; Sigma pl ; WLC-Q ; eFJC Cont. length (nm) ; Sigma eFJC cl ; Kuhn length (nm); Sigma kl ; FJC-Q ;Force (pN) ; Slope (N/m) ; (untested) Loading rate (pN/s)\n')
283 if not os.path.exists(self.autofile):
284 f=open(self.autofile,'a+')
285 f.write('Analysis started '+time.asctime()+'\n')
286 f.write('----------------------------------------\n')
287 f.write(' File ; WLC Cont. length (nm) ; Sigma WLC cl; Per. Length (nm) ; Sigma pl; WLC-Q ;eFJC Cont. length (nm) ; Sigma eFJC cl ; Kuhn length (nm); Sigma kl ; FJC-Q ;Force (pN) ; Slope (N/m) ; (untested) Loading rate (pN/s)\n')
292 f=open(self.autofile,'a+')
293 f.write(self.current.path)
294 f.write(' ; '+str(1e9*wlcparams[0])+' ; '+str(wlcfit_nm[0])+' ; '+str(1e9*wlcparams[1])+' ; '+str(wlcfit_nm[1])+' ; '+ str(wlc_fitq)+' ; '+str(1e9*fjcparams[0])+' ; '+str(fjcfit_nm[0])+' ; '+str(1e9*fjcparams[1])+' ; '+str(fjcfit_nm[1])+' ; '+str(fjc_fitq)+' ; '+str(1e12*force)+' ; '+ str(slope)+' ; '+str(lrate)+'\n')
297 print '\nWould you like to try again on this same curve?'