hooke/plugin/multifit.py

   1 # Copyright (C) 2010-2012 W. Trevor King <wking@drexel.edu>
   2 #
   3 # This file is part of Hooke.
   4 #
   5 # Hooke is free software: you can redistribute it and/or modify it under the
   6 # terms of the GNU Lesser General Public License as published by the Free
   7 # Software Foundation, either version 3 of the License, or (at your option) any
   8 # later version.
   9 #
  10 # Hooke is distributed in the hope that it will be useful, but WITHOUT ANY
  11 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  12 # A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
  13 # details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public License
  16 # along with Hooke.  If not, see <http://www.gnu.org/licenses/>.
  17
  18 """The ``multifit`` module provides :class:`MultifitPlugin` and
  19 several associated :class:`hooke.command.Command`\s for fitting
  20 all sorts of :class:`hooke.experiment.VelocityClamp` parameters.
  21 """
  22
  23 from ..command import Command, Argument, Failure
  24 from ..playlist import FilePlaylist
  25 from . import Plugin
  26 from .playlist import PlaylistArgument
  27
  28
  29 class MultifitPlugin (Plugin):
  30     def __init__(self):
  31         super(MultifitPlugin, self).__init__(name='multifit')
  32         self._commands = [MultifitCommand(self)]
  33
  34
  35 class MultifitCommand (Command):
  36     """Presents curves for interactive manual analysis.
  37
  38     Obtains contour length, persistance length, rupture force and
  39     slope (loading rate).
  40
  41     See :class:`hooke.plugin.fit.FitCommand` help for more information
  42     on the options and fit procedures.
  43     """
  44     #NOTE: centerzero plot modifier should be activated (set centerzero
  45     #1).
  46     def __init__(self, plugin):
  47         super(MultifitCommand, self).__init__(
  48             name='multifit',
  49             arguments=[
  50                 PlaylistArgument,
  51                 Argument(name='persistence length', type='float',
  52                     help="""
  53 Persistence length in meters for WLC fits.  If not given, the fit will
  54 be a 2-variable fit.
  55 """.strip()),
  56                 Argument(name='Kuhn length', type='float',
  57                     help="""
  58 Kuhn length in meters for eFJC fits.  If not given, the fit will be a
  59 2-variable fit.  The eFJC fit works better with a fixed Kuhn length
  60 (citation?).
  61 """.strip()),
  62                 Argument(name='temperature', type='float', default=293.0,
  63                     help="""
  64 Temperature in Kelvin.  Defaults to 293 K.
  65 """.strip()),
  66                 Argument(name='slope width', type='int',
  67                     help="""
  68 Width in points for slope fitting (points leading up to the rupture).
  69 If not given, the routine will prompt for an extra point to fit the
  70 slope.
  71 """.strip()),
  72                 Argument(name='base width', type='int', default=15,
  73                     help="""
  74 Width in points for slope fitting (points following the rupture?).
  75 """.strip()),
  76             Argument(name='point slope', type='bool', default=False,
  77                      help="""
  78 Calculate the slope from the two selected points instead of fitting
  79 data between them.
  80 """.strip()),
  81             Argument(name='just one', type='bool', default=False,
  82                      help="""
  83 Perform the fits on the current curve instead of iterating.
  84 """.strip()),
  85                 ],
  86             help=self.__doc__, plugin=plugin)
  87
  88     def _run(self, hooke, inqueue, outqueue, params):
  89         counter=0
  90         savecounter=0
  91         curveindex=0
  92
  93         start_curve = params['playlist'].current()
  94         curve = start_curve
  95         while True:
  96             outqueue.put('TODO to end %s execution' % self.name)
  97             self._fit_curve(hooke, inqueue, outqueue, params, curve)
  98             if params['just one'] == True:
  99                 break
 100             params['playlist'].next()
 101             next_curve = params['playlist'].current()
 102             if next_curve == start_curve:
 103                 break # loop through playlist complete
 104             curve = next_curve
 105
 106     def _fit_curve(self, hooke, inqueue, outqueue, params, curve):
 107         contact = self._get_point('Mark the contact point of the curve.')
 108         #hightlights an area dedicated to reject current curve
 109         sizeret=len(self.plots[0].vectors[1][0])
 110         noarea=self.plots[0].vectors[1][0][-sizeret/6:-1]
 111         noareaplot=copy.deepcopy(self._get_displayed_plot())
 112         #noise dependent slight shift to improve visibility
 113         noareaplot.add_set(noarea,np.ones_like(noarea)*5.0*np.std(self.plots[0].vectors[1][1][-sizeret/6:-1]))
 114         noareaplot.styles+=['scatter']
 115         noareaplot.colors+=["red"]
 116         self._send_plot([noareaplot])
 117         contact_point=self._measure_N_points(N=1, whatset=1)[0]
 118         contact_point_index=contact_point.index
 119         noareaplot.remove_set(-1)
 120         #remove set leaves some styles disturbing the next graph, fixed by doing do_plot
 121         self.do_plot(0)
 122         if contact_point_index>5.0/6.0*sizeret:
 123             return
 124
 125         self.wlccontact_point=contact_point
 126         self.wlccontact_index=contact_point.index
 127         self.wlccurrent=self.current.path
 128
 129         print 'Now click two points for the fitting area (one should be the rupture point)'
 130         wlcpoints=self._measure_N_points(N=2,whatset=1)
 131         print 'And one point of the top of the jump'
 132         toppoint=self._measure_N_points(N=1,whatset=1)
 133         slopepoint=ClickedPoint()
 134         if slopew==None:
 135             print 'Click a point to calculate slope'
 136             slopepoint=self._measure_N_points(N=1,whatset=1)
 137
 138         fitpoints=[contact_point]+wlcpoints
 139         #use the currently displayed plot for the fit
 140         displayed_plot=self._get_displayed_plot()
 141
 142         #use both fit functions
 143         try:
 144             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 )
 145             wlcerror=False
 146         except:
 147             print 'WLC fit not possible'
 148             wlcerror=True
 149
 150         try:
 151             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 )
 152             fjcerror=False
 153         except:
 154             print 'eFJC fit not possible'
 155             fjcerror=True
 156
 157         #Measure rupture force
 158         ruptpoint=ClickedPoint()
 159         if wlcpoints[0].graph_coords[1]<wlcpoints[1].graph_coords[1]:
 160             ruptpoint=wlcpoints[0]
 161         else:
 162             ruptpoint=wlcpoints[1]
 163         tpindex=toppoint[0].index
 164         toplevel=np.average(displayed_plot.vectors[1][1][tpindex:tpindex+basew])
 165         force=toplevel-ruptpoint.graph_coords[1]
 166
 167         #Measure the slope - loading rate
 168
 169         if slopew==None:
 170             if twops:
 171                 xint=ruptpoint.graph_coords[0]-slopepoint[0].graph_coords[0]
 172                 yint=ruptpoint.graph_coords[1]-slopepoint[0].graph_coords[1]
 173                 slope=yint/xint
 174                 slopeplot=self._get_displayed_plot(0) #get topmost displayed plot
 175                 slopeplot.add_set([ruptpoint.graph_coords[0],slopepoint[0].graph_coords[0]],[ruptpoint.graph_coords[1],slopepoint[0].graph_coords[1]])
 176                 if slopeplot.styles==[]:
 177                     slopeplot.styles=[None,None,'scatter']
 178                     slopeplot.colors=[None,None,'red']
 179                 else:
 180                     slopeplot.styles+=['scatter']
 181                     slopeplot.colors+=['red']
 182             else:
 183                 slope=self._slope([ruptpoint]+slopepoint,slopew)
 184         else:
 185             slope=self._slope([ruptpoint],slopew)
 186
 187         #plot results (_slope already did)
 188
 189         #now we have the fit, we can plot it
 190         #add the clicked points in the final PlotObject
 191         clickvector_x, clickvector_y=[], []
 192         for item in wlcpoints:
 193             clickvector_x.append(item.graph_coords[0])
 194             clickvector_y.append(item.graph_coords[1])
 195
 196         #create a custom PlotObject to gracefully plot the fit along the curves
 197
 198         #first fix those irritating zoom-destroying fits
 199         lowestpoint=min(displayed_plot.vectors[1][1])
 200         for i in np.arange(0,len(wlcyfit)):
 201             if wlcyfit[i] < 1.2*lowestpoint:
 202                 wlcyfit[i]=toplevel
 203         for i in np.arange(0,len(fjcyfit)):
 204             if fjcyfit[i] < 1.2*lowestpoint:
 205                 fjcyfit[i]=toplevel
 206
 207         fitplot=copy.deepcopy(displayed_plot)
 208         fitplot.add_set(wlcxfit,wlcyfit)
 209         fitplot.add_set(fjcxfit,fjcyfit)
 210         fitplot.add_set(clickvector_x,clickvector_y)
 211
 212         fitplot.styles+=[None,None,'scatter']
 213         fitplot.colors+=["red","blue",None]
 214
 215         self._send_plot([fitplot])
 216
 217
 218         #present results of measurement
 219         if len(wlcparams)==1:
 220             wlcparams.append(pl_value*1e-9)
 221         if len(fjcparams)==1:
 222             fjcparams.append(kl_value*1e-9)
 223
 224         if fjcfit_errors:
 225             fjcfit_nm=[i*(10**9) for i in fjcfit_errors]
 226             if len(fjcfit_nm)==1:
 227                 fjcfit_nm.append(0)
 228         else:
 229             fjcfit_errors=[0,0]
 230         if wlcfit_errors:
 231             wlcfit_nm=[i*(10**9) for i in wlcfit_errors]
 232         if len(wlcfit_nm)==1:
 233             wlcfit_nm.append(0)
 234         else:
 235             wlcfit_errors=[0,0]
 236
 237         wlc_fitq=wlc_qstd/np.std(displayed_plot.vectors[1][1][-20:-1])
 238         fjc_fitq=fjc_qstd/np.std(displayed_plot.vectors[1][1][-20:-1])
 239
 240         print '\nRESULTS'
 241         print 'WLC contour : '+str(1e9*wlcparams[0])+u' \u00b1 '+str(wlcfit_nm[0])+' nm'
 242         print 'Per. length : '+str(1e9*wlcparams[1])+u' \u00b1 '+str(wlcfit_nm[1])+' nm'
 243         print 'Quality :'+str(wlc_fitq)
 244         print '---'
 245         print 'eFJC contour : '+str(1e9*fjcparams[0])+u' \u00b1 '+str(fjcfit_nm[0])+' nm'
 246         print 'Kuhn length (e): '+str(1e9*fjcparams[1])+u' \u00b1 '+str(fjcfit_nm[1])+' nm'
 247         print 'Quality :'+str(fjc_fitq)
 248         print '---'
 249         print 'Force : '+str(1e12*force)+' pN'
 250         print 'Slope : '+str(slope)+' N/m'
 251
 252         try:
 253             #FIXME all drivers should provide retract velocity, in SI or homogeneous units
 254             lrate=slope*self.current.curve.retract_velocity
 255             print 'Loading rate (UNTESTED):'+str(lrate)
 256         except:
 257             print 'Loading rate : n/a'
 258             lrate='n/a'
 259
 260         #save accept/discard/repeat
 261         print '\n\nDo you want to save these?'
 262         if self.YNclick():
 263
 264         #Save file info
 265             if self.autofile=='':
 266                 self.autofile=raw_input('Results filename? (press enter for a random generated name)')
 267                 if self.autofile=='':
 268                     [osf,name]=tempfile.mkstemp(prefix='analysis-',suffix='.txt',text=True,dir=self.config['hookedir'])
 269                     print 'saving in '+name
 270                     self.autofile=name
 271                     os.close(osf)
 272                     f=open(self.autofile,'a+')
 273                     f.write('Analysis started '+time.asctime()+'\n')
 274                     f.write('----------------------------------------\n')
 275                     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')
 276                     f.close()
 277
 278             if not os.path.exists(self.autofile):
 279                 f=open(self.autofile,'a+')
 280                 f.write('Analysis started '+time.asctime()+'\n')
 281                 f.write('----------------------------------------\n')
 282                 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                 f.close()
 284
 285             print 'Saving...'
 286             savecounter+=1
 287             f=open(self.autofile,'a+')
 288             f.write(self.current.path)
 289             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')
 290             f.close()
 291         else:
 292             print '\nWould you like to try again on this same curve?'
 293             if self.YNclick():
 294                 return