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