hooke/plugin/multifit.py

   1 # Copyright (C) 2010-2011 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 . import Plugin
  28 from .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             return
 126
 127         self.wlccontact_point=contact_point
 128         self.wlccontact_index=contact_point.index
 129         self.wlccurrent=self.current.path
 130
 131         print 'Now click two points for the fitting area (one should be the rupture point)'
 132         wlcpoints=self._measure_N_points(N=2,whatset=1)
 133         print 'And one point of the top of the jump'
 134         toppoint=self._measure_N_points(N=1,whatset=1)
 135         slopepoint=ClickedPoint()
 136         if slopew==None:
 137             print 'Click a point to calculate slope'
 138             slopepoint=self._measure_N_points(N=1,whatset=1)
 139
 140         fitpoints=[contact_point]+wlcpoints
 141         #use the currently displayed plot for the fit
 142         displayed_plot=self._get_displayed_plot()
 143
 144         #use both fit functions
 145         try:
 146             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 )
 147             wlcerror=False
 148         except:
 149             print 'WLC fit not possible'
 150             wlcerror=True
 151
 152         try:
 153             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 )
 154             fjcerror=False
 155         except:
 156             print 'eFJC fit not possible'
 157             fjcerror=True
 158
 159         #Measure rupture force
 160         ruptpoint=ClickedPoint()
 161         if wlcpoints[0].graph_coords[1]<wlcpoints[1].graph_coords[1]:
 162             ruptpoint=wlcpoints[0]
 163         else:
 164             ruptpoint=wlcpoints[1]
 165         tpindex=toppoint[0].index
 166         toplevel=np.average(displayed_plot.vectors[1][1][tpindex:tpindex+basew])
 167         force=toplevel-ruptpoint.graph_coords[1]
 168
 169         #Measure the slope - loading rate
 170
 171         if slopew==None:
 172             if twops:
 173                 xint=ruptpoint.graph_coords[0]-slopepoint[0].graph_coords[0]
 174                 yint=ruptpoint.graph_coords[1]-slopepoint[0].graph_coords[1]
 175                 slope=yint/xint
 176                 slopeplot=self._get_displayed_plot(0) #get topmost displayed plot
 177                 slopeplot.add_set([ruptpoint.graph_coords[0],slopepoint[0].graph_coords[0]],[ruptpoint.graph_coords[1],slopepoint[0].graph_coords[1]])
 178                 if slopeplot.styles==[]:
 179                     slopeplot.styles=[None,None,'scatter']
 180                     slopeplot.colors=[None,None,'red']
 181                 else:
 182                     slopeplot.styles+=['scatter']
 183                     slopeplot.colors+=['red']
 184             else:
 185                 slope=self._slope([ruptpoint]+slopepoint,slopew)
 186         else:
 187             slope=self._slope([ruptpoint],slopew)
 188
 189         #plot results (_slope already did)
 190
 191         #now we have the fit, we can plot it
 192         #add the clicked points in the final PlotObject
 193         clickvector_x, clickvector_y=[], []
 194         for item in wlcpoints:
 195             clickvector_x.append(item.graph_coords[0])
 196             clickvector_y.append(item.graph_coords[1])
 197
 198         #create a custom PlotObject to gracefully plot the fit along the curves
 199
 200         #first fix those irritating zoom-destroying fits
 201         lowestpoint=min(displayed_plot.vectors[1][1])
 202         for i in np.arange(0,len(wlcyfit)):
 203             if wlcyfit[i] < 1.2*lowestpoint:
 204                 wlcyfit[i]=toplevel
 205         for i in np.arange(0,len(fjcyfit)):
 206             if fjcyfit[i] < 1.2*lowestpoint:
 207                 fjcyfit[i]=toplevel
 208
 209         fitplot=copy.deepcopy(displayed_plot)
 210         fitplot.add_set(wlcxfit,wlcyfit)
 211         fitplot.add_set(fjcxfit,fjcyfit)
 212         fitplot.add_set(clickvector_x,clickvector_y)
 213
 214         fitplot.styles+=[None,None,'scatter']
 215         fitplot.colors+=["red","blue",None]
 216
 217         self._send_plot([fitplot])
 218
 219
 220         #present results of measurement
 221         if len(wlcparams)==1:
 222             wlcparams.append(pl_value*1e-9)
 223         if len(fjcparams)==1:
 224             fjcparams.append(kl_value*1e-9)
 225
 226         if fjcfit_errors:
 227             fjcfit_nm=[i*(10**9) for i in fjcfit_errors]
 228             if len(fjcfit_nm)==1:
 229                 fjcfit_nm.append(0)
 230         else:
 231             fjcfit_errors=[0,0]
 232         if wlcfit_errors:
 233             wlcfit_nm=[i*(10**9) for i in wlcfit_errors]
 234         if len(wlcfit_nm)==1:
 235             wlcfit_nm.append(0)
 236         else:
 237             wlcfit_errors=[0,0]
 238
 239         wlc_fitq=wlc_qstd/np.std(displayed_plot.vectors[1][1][-20:-1])
 240         fjc_fitq=fjc_qstd/np.std(displayed_plot.vectors[1][1][-20:-1])
 241
 242         print '\nRESULTS'
 243         print 'WLC contour : '+str(1e9*wlcparams[0])+u' \u00b1 '+str(wlcfit_nm[0])+' nm'
 244         print 'Per. length : '+str(1e9*wlcparams[1])+u' \u00b1 '+str(wlcfit_nm[1])+' nm'
 245         print 'Quality :'+str(wlc_fitq)
 246         print '---'
 247         print 'eFJC contour : '+str(1e9*fjcparams[0])+u' \u00b1 '+str(fjcfit_nm[0])+' nm'
 248         print 'Kuhn length (e): '+str(1e9*fjcparams[1])+u' \u00b1 '+str(fjcfit_nm[1])+' nm'
 249         print 'Quality :'+str(fjc_fitq)
 250         print '---'
 251         print 'Force : '+str(1e12*force)+' pN'
 252         print 'Slope : '+str(slope)+' N/m'
 253
 254         try:
 255             #FIXME all drivers should provide retract velocity, in SI or homogeneous units
 256             lrate=slope*self.current.curve.retract_velocity
 257             print 'Loading rate (UNTESTED):'+str(lrate)
 258         except:
 259             print 'Loading rate : n/a'
 260             lrate='n/a'
 261
 262         #save accept/discard/repeat
 263         print '\n\nDo you want to save these?'
 264         if self.YNclick():
 265
 266         #Save file info
 267             if self.autofile=='':
 268                 self.autofile=raw_input('Results filename? (press enter for a random generated name)')
 269                 if self.autofile=='':
 270                     [osf,name]=tempfile.mkstemp(prefix='analysis-',suffix='.txt',text=True,dir=self.config['hookedir'])
 271                     print 'saving in '+name
 272                     self.autofile=name
 273                     os.close(osf)
 274                     f=open(self.autofile,'a+')
 275                     f.write('Analysis started '+time.asctime()+'\n')
 276                     f.write('----------------------------------------\n')
 277                     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')
 278                     f.close()
 279
 280             if not os.path.exists(self.autofile):
 281                 f=open(self.autofile,'a+')
 282                 f.write('Analysis started '+time.asctime()+'\n')
 283                 f.write('----------------------------------------\n')
 284                 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')
 285                 f.close()
 286
 287             print 'Saving...'
 288             savecounter+=1
 289             f=open(self.autofile,'a+')
 290             f.write(self.current.path)
 291             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')
 292             f.close()
 293         else:
 294             print '\nWould you like to try again on this same curve?'
 295             if self.YNclick():
 296                 continue
 297         curveindex+=1