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