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 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.
  10 #
  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.
  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