56befcaeef7f8c9d6d757993029da4f77e14e617
[hooke.git] / hooke / plugin / multifit.py
1 # Copyright (C) 2010 Alberto Gomez-Casado
2 #                    W. Trevor King <wking@drexel.edu>
3 #
4 # This file is part of Hooke.
5 #
6 # Hooke is free software: you can redistribute it and/or
7 # modify it under the terms of the GNU Lesser General Public
8 # License as published by the Free Software Foundation, either
9 # version 3 of the License, or (at your option) any later version.
10 #
11 # Hooke is distributed in the hope that it will be useful,
12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 # GNU Lesser General Public License for more details.
15 #
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/>.
19
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.
23 """
24
25 from ..command import Command, Argument, Failure
26 from ..playlist import FilePlaylist
27 from ..plugin import Plugin
28 from ..plugin.playlist import PlaylistArgument
29
30
31 class MultifitPlugin (Plugin):
32     def __init__(self):
33         super(MultifitPlugin, self).__init__(name='multifit')
34         self._commands = [MultifitCommand(self)]
35
36
37 class MultifitCommand (Command):
38     """Presents curves for interactive manual analysis.
39
40     Obtains contour length, persistance length, rupture force and 
41     slope (loading rate).
42
43     See :class:`hooke.plugin.fit.FitCommand` help for more information
44     on the options and fit procedures.
45     """
46     #NOTE: centerzero plot modifier should be activated (set centerzero
47     #1).
48     def __init__(self, plugin):
49         super(MultifitCommand, self).__init__(
50             name='multifit',
51             arguments=[
52                 PlaylistArgument,
53                 Argument(name='persistence length', type='float',
54                     help="""
55 Persistence length in meters for WLC fits.  If not given, the fit will
56 be a 2-variable fit.
57 """.strip()),
58                 Argument(name='Kuhn length', type='float',
59                     help="""
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
62 (citation?).
63 """.strip()),
64                 Argument(name='temperature', type='float', default=293.0,
65                     help="""
66 Temperature in Kelvin.  Defaults to 293 K.
67 """.strip()),
68                 Argument(name='slope width', type='int',
69                     help="""
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
72 slope.
73 """.strip()),
74                 Argument(name='base width', type='int', default=15,
75                     help="""
76 Width in points for slope fitting (points following the rupture?).
77 """.strip()),
78             Argument(name='point slope', type='bool', default=False,
79                      help="""
80 Calculate the slope from the two selected points instead of fitting
81 data between them.
82 """.strip()),
83             Argument(name='just one', type='bool', default=False,
84                      help="""
85 Perform the fits on the current curve instead of iterating.
86 """.strip()),
87                 ],
88             help=self.__doc__, plugin=plugin)
89
90     def _run(self, hooke, inqueue, outqueue, params):
91         counter=0
92         savecounter=0
93         curveindex=0
94         
95         start_curve = params['playlist'].current()
96         curve = start_curve
97         while True:
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:
101                 break
102             params['playlist'].next()
103             next_curve = params['playlist'].current()
104             if next_curve == start_curve:
105                 break # loop through playlist complete
106             curve = next_curve
107
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
123         self.do_plot(0)
124         if contact_point_index>5.0/6.0*sizeret:
125             if params['just one']:
126                 break
127             curveindex+=1
128             continue                            
129         
130         self.wlccontact_point=contact_point
131         self.wlccontact_index=contact_point.index
132         self.wlccurrent=self.current.path
133         
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()
139         if slopew==None:
140             print 'Click a point to calculate slope'
141             slopepoint=self._measure_N_points(N=1,whatset=1)
142             
143         fitpoints=[contact_point]+wlcpoints
144         #use the currently displayed plot for the fit
145         displayed_plot=self._get_displayed_plot()
146         
147         #use both fit functions
148         try:
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 )
150             wlcerror=False      
151         except:
152             print 'WLC fit not possible'
153             wlcerror=True
154             
155         try:
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 )
157             fjcerror=False
158         except:
159             print 'eFJC fit not possible'
160             fjcerror=True
161                 
162         #Measure rupture force
163         ruptpoint=ClickedPoint()    
164         if wlcpoints[0].graph_coords[1]<wlcpoints[1].graph_coords[1]:
165             ruptpoint=wlcpoints[0]
166         else:
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]                                        
171         
172         #Measure the slope - loading rate
173             
174         if slopew==None:
175             if twops:
176                 xint=ruptpoint.graph_coords[0]-slopepoint[0].graph_coords[0]
177                 yint=ruptpoint.graph_coords[1]-slopepoint[0].graph_coords[1]
178                 slope=yint/xint
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']
184                 else:
185                     slopeplot.styles+=['scatter']
186                     slopeplot.colors+=['red']
187             else:    
188                 slope=self._slope([ruptpoint]+slopepoint,slopew)
189         else:
190             slope=self._slope([ruptpoint],slopew)
191          
192         #plot results (_slope already did)
193             
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])
200             
201         #create a custom PlotObject to gracefully plot the fit along the curves
202             
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:
207                 wlcyfit[i]=toplevel    
208         for i in np.arange(0,len(fjcyfit)):
209             if fjcyfit[i] < 1.2*lowestpoint:
210                 fjcyfit[i]=toplevel
211             
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)
216
217         fitplot.styles+=[None,None,'scatter']
218         fitplot.colors+=["red","blue",None]
219             
220         self._send_plot([fitplot])
221         
222
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)
228                 
229         if fjcfit_errors:
230             fjcfit_nm=[i*(10**9) for i in fjcfit_errors]
231             if len(fjcfit_nm)==1:
232                 fjcfit_nm.append(0)
233         else:
234             fjcfit_errors=[0,0]        
235         if wlcfit_errors:
236             wlcfit_nm=[i*(10**9) for i in wlcfit_errors]
237         if len(wlcfit_nm)==1:
238             wlcfit_nm.append(0)
239         else:
240             wlcfit_errors=[0,0]
241             
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])
244             
245         print '\nRESULTS'
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)
249         print '---'
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)   
253         print '---'
254         print 'Force : '+str(1e12*force)+' pN'
255         print 'Slope : '+str(slope)+' N/m'
256
257         try:
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)
261         except:
262             print 'Loading rate : n/a'
263             lrate='n/a'
264             
265         #save accept/discard/repeat
266         print '\n\nDo you want to save these?'
267         if self.YNclick():
268                     
269         #Save file info
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
275                     self.autofile=name
276                     os.close(osf)
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')
281                     f.close()
282                     
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')
288                 f.close()
289                 
290             print 'Saving...'
291             savecounter+=1
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')
295             f.close()
296         else:
297             print '\nWould you like to try again on this same curve?'
298             if self.YNclick():
299                 continue
300         curveindex+=1