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 #FIXME clean this, probably some dependencies are not needed
  21
  22 from ..libhooke import WX_GOOD, ClickedPoint
  23 import wxversion
  24 wxversion.select(WX_GOOD)
  25 from wx import PostEvent
  26 import numpy as np
  27 import scipy as sp
  28 import copy
  29 import os.path
  30 import time
  31 import tempfile
  32 import warnings
  33 warnings.simplefilter('ignore',np.RankWarning)
  34 import Queue
  35
  36 global measure_wlc
  37 global EVT_MEASURE_WLC
  38
  39 global events_from_fit
  40 events_from_fit=Queue.Queue() #GUI ---> CLI COMMUNICATION
  41
  42
  43 class multifitCommands:
  44
  45     def do_multifit(self,args):
  46         '''
  47         MULTIFIT
  48         (multifit.py)
  49         Presents curves for manual analysis in a comfortable mouse-only fashion.
  50         Obtains contour length, persistance length, rupture force and
  51         slope - loading rate.
  52         WLC is shown in red, eFJC in blue.
  53         -------------
  54         Syntax:
  55         multifit [pl=value] [kl=value] [t=value] [slopew=value] [basew=value]
  56                 [slope2p] [justone]
  57
  58         pl=[value] and kl=[value]: Use a fixed persistent length (WLC) or Kuhn
  59                 length (eFJC) for the fit. If pl is not given, the fit will be
  60                 a 2-variable fit.
  61                 DO NOT put spaces between 'pl', '=' and the value.
  62                 eFJC fit works better with a fixed kl
  63                 The value must be in nanometers.
  64
  65         t=[value] : Use a user-defined temperature. The value must be in
  66                 kelvins; by default it is 293 K.
  67                 DO NOT put spaces between 't', '=' and the value.
  68
  69         slopew and basew : width in points for slope fitting (points to the
  70                 right of clicked rupture) and base level fitting (points to
  71                 the left of clicked top of rupture).
  72                 If slopew is not provided, the routine will ask for a 5th
  73                 point to fit the slope.
  74                 DO NOT put spaces between 'slopew' or 'basew', '=' value.
  75
  76         slope2p : calculates the slope from the two clicked points,
  77                 instead of fitting data between them
  78
  79         justone : performs the fits over current curve instead of iterating
  80
  81         See fit command help for more information on the options and fit
  82         procedures.
  83         NOTE: centerzero plot modifier should be activated (set centerzero 1).
  84         '''
  85
  86                 #NOTE duplicates a lot of code from do_fit in fit.py, a call to it could be used directly
  87                 #but it is easier to control the program flow bypassing it
  88
  89         pl_value=None
  90         kl_value=None
  91         T=self.config['temperature']
  92         slopew=None
  93         basew=15
  94         justone=False
  95         twops=False
  96
  97         #FIXME spaces are not allowed between pl or t and value
  98         for arg in args.split():
  99             #look for a persistent length argument.
 100             if 'pl=' in arg:
 101                 pl_expression=arg.split('=')
 102                 pl_value=float(pl_expression[1]) #actual value
 103             if 'kl=' in arg:
 104                 kl_expression=arg.split('=')
 105                 kl_value=float(kl_expression[1]) #actual value
 106             #look for a T argument.
 107             if ('t=' in arg) or ('T=' in arg):
 108                 t_expression=arg.split('=')
 109                 T=float(t_expression[1])
 110             #look for a basew argument.
 111             if ('basew=' in arg):
 112                 basew_expression=arg.split('=')
 113                 basew=int(basew_expression[1])
 114             #look for a slopew argument.
 115             if ('slopew=' in arg):
 116                 slopew_expression=arg.split('=')
 117                 slopew=int(slopew_expression[1])
 118             if('justone' in arg):
 119                 justone=True
 120             if('slope2p' in arg):
 121                 twops=True
 122
 123         counter=0
 124         savecounter=0
 125         curveindex=0
 126
 127         if not justone:
 128             print 'What curve no. you would like to start? (enter for ignoring)'
 129             skip=raw_input()
 130
 131             if skip.isdigit()==False:
 132                 skip=0
 133             else:
 134                 skip=int(skip)-1
 135                 print 'Skipping '+str(skip)+ ' curves'
 136         else:
 137             skip=0
 138         #begin presenting curves for analysis
 139         while curveindex <len(self.current_list):
 140             if not justone:
 141                 counter+=1
 142                 curve=self.current_list[curveindex]
 143                 if skip>curveindex:
 144                     curveindex+=1
 145                     continue
 146
 147             #give periodically the opportunity to stop the analysis
 148                 if counter%20==0 and counter>0:
 149                     print '\nYou checked '+str(counter)+' curves. Do you want to continue?'
 150                     self.current=curve
 151                     self.do_plot(0)
 152                     if self.YNclick():
 153                         print 'Going on...'
 154                     else:
 155                         break
 156                 else:
 157                     self.current=curve
 158                     self.do_plot(0)
 159             else:
 160                 curve=self.current
 161                 self.do_plot(0)
 162             if not justone:
 163                 print '\nCurve '+str(curveindex+1)+' of '+str(len(self.current_list))
 164             print 'Click contact point or left end (red area)of the curve to skip'
 165             #hightlights an area dedicated to reject current curve
 166             sizeret=len(self.plots[0].vectors[1][0])
 167             noarea=self.plots[0].vectors[1][0][-sizeret/6:-1]
 168             noareaplot=copy.deepcopy(self._get_displayed_plot())
 169             #noise dependent slight shift to improve visibility
 170             noareaplot.add_set(noarea,np.ones_like(noarea)*5.0*np.std(self.plots[0].vectors[1][1][-sizeret/6:-1]))
 171             noareaplot.styles+=['scatter']
 172             noareaplot.colors+=["red"]
 173             self._send_plot([noareaplot])
 174             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 175             contact_point_index=contact_point.index
 176             noareaplot.remove_set(-1)
 177             #remove set leaves some styles disturbing the next graph, fixed by doing do_plot
 178             self.do_plot(0)
 179             if contact_point_index>5.0/6.0*sizeret:
 180                 if justone:
 181                     break
 182                 curveindex+=1
 183                 continue
 184
 185             self.wlccontact_point=contact_point
 186             self.wlccontact_index=contact_point.index
 187             self.wlccurrent=self.current.path
 188
 189             print 'Now click two points for the fitting area (one should be the rupture point)'
 190             wlcpoints=self._measure_N_points(N=2,whatset=1)
 191             print 'And one point of the top of the jump'
 192             toppoint=self._measure_N_points(N=1,whatset=1)
 193             slopepoint=ClickedPoint()
 194             if slopew==None:
 195                 print 'Click a point to calculate slope'
 196                 slopepoint=self._measure_N_points(N=1,whatset=1)
 197
 198             fitpoints=[contact_point]+wlcpoints
 199             #use the currently displayed plot for the fit
 200             displayed_plot=self._get_displayed_plot()
 201
 202             #use both fit functions
 203             try:
 204                 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 )
 205                 wlcerror=False
 206             except:
 207                 print 'WLC fit not possible'
 208                 wlcerror=True
 209
 210             try:
 211                 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 )
 212                 fjcerror=False
 213             except:
 214                 print 'eFJC fit not possible'
 215                 fjcerror=True
 216
 217             #Measure rupture force
 218             ruptpoint=ClickedPoint()
 219             if wlcpoints[0].graph_coords[1]<wlcpoints[1].graph_coords[1]:
 220                 ruptpoint=wlcpoints[0]
 221             else:
 222                 ruptpoint=wlcpoints[1]
 223             tpindex=toppoint[0].index
 224             toplevel=np.average(displayed_plot.vectors[1][1][tpindex:tpindex+basew])
 225             force=toplevel-ruptpoint.graph_coords[1]
 226
 227             #Measure the slope - loading rate
 228
 229             if slopew==None:
 230                 if twops:
 231                     xint=ruptpoint.graph_coords[0]-slopepoint[0].graph_coords[0]
 232                     yint=ruptpoint.graph_coords[1]-slopepoint[0].graph_coords[1]
 233                     slope=yint/xint
 234                     slopeplot=self._get_displayed_plot(0) #get topmost displayed plot
 235                     slopeplot.add_set([ruptpoint.graph_coords[0],slopepoint[0].graph_coords[0]],[ruptpoint.graph_coords[1],slopepoint[0].graph_coords[1]])
 236                     if slopeplot.styles==[]:
 237                         slopeplot.styles=[None,None,'scatter']
 238                         slopeplot.colors=[None,None,'red']
 239                     else:
 240                         slopeplot.styles+=['scatter']
 241                         slopeplot.colors+=['red']
 242                 else:
 243                     slope=self._slope([ruptpoint]+slopepoint,slopew)
 244             else:
 245                 slope=self._slope([ruptpoint],slopew)
 246
 247             #plot results (_slope already did)
 248
 249             #now we have the fit, we can plot it
 250             #add the clicked points in the final PlotObject
 251             clickvector_x, clickvector_y=[], []
 252             for item in wlcpoints:
 253                 clickvector_x.append(item.graph_coords[0])
 254                 clickvector_y.append(item.graph_coords[1])
 255
 256             #create a custom PlotObject to gracefully plot the fit along the curves
 257
 258             #first fix those irritating zoom-destroying fits
 259             lowestpoint=min(displayed_plot.vectors[1][1])
 260             for i in np.arange(0,len(wlcyfit)):
 261                 if wlcyfit[i] < 1.2*lowestpoint:
 262                     wlcyfit[i]=toplevel
 263             for i in np.arange(0,len(fjcyfit)):
 264                 if fjcyfit[i] < 1.2*lowestpoint:
 265                     fjcyfit[i]=toplevel
 266
 267             fitplot=copy.deepcopy(displayed_plot)
 268             fitplot.add_set(wlcxfit,wlcyfit)
 269             fitplot.add_set(fjcxfit,fjcyfit)
 270             fitplot.add_set(clickvector_x,clickvector_y)
 271
 272             fitplot.styles+=[None,None,'scatter']
 273             fitplot.colors+=["red","blue",None]
 274
 275             self._send_plot([fitplot])
 276
 277
 278              #present results of measurement
 279             if len(wlcparams)==1:
 280                 wlcparams.append(pl_value*1e-9)
 281             if len(fjcparams)==1:
 282                 fjcparams.append(kl_value*1e-9)
 283
 284             if fjcfit_errors:
 285                 fjcfit_nm=[i*(10**9) for i in fjcfit_errors]
 286                 if len(fjcfit_nm)==1:
 287                     fjcfit_nm.append(0)
 288             else:
 289                 fjcfit_errors=[0,0]
 290             if wlcfit_errors:
 291                 wlcfit_nm=[i*(10**9) for i in wlcfit_errors]
 292                 if len(wlcfit_nm)==1:
 293                     wlcfit_nm.append(0)
 294             else:
 295                 wlcfit_errors=[0,0]
 296
 297             wlc_fitq=wlc_qstd/np.std(displayed_plot.vectors[1][1][-20:-1])
 298             fjc_fitq=fjc_qstd/np.std(displayed_plot.vectors[1][1][-20:-1])
 299
 300             print '\nRESULTS'
 301             print 'WLC contour : '+str(1e9*wlcparams[0])+u' \u00b1 '+str(wlcfit_nm[0])+' nm'
 302             print 'Per. length : '+str(1e9*wlcparams[1])+u' \u00b1 '+str(wlcfit_nm[1])+' nm'
 303             print 'Quality :'+str(wlc_fitq)
 304             print '---'
 305             print 'eFJC contour : '+str(1e9*fjcparams[0])+u' \u00b1 '+str(fjcfit_nm[0])+' nm'
 306             print 'Kuhn length (e): '+str(1e9*fjcparams[1])+u' \u00b1 '+str(fjcfit_nm[1])+' nm'
 307             print 'Quality :'+str(fjc_fitq)
 308             print '---'
 309             print 'Force : '+str(1e12*force)+' pN'
 310             print 'Slope : '+str(slope)+' N/m'
 311
 312             try:
 313                 #FIXME all drivers should provide retract velocity, in SI or homogeneous units
 314                 lrate=slope*self.current.curve.retract_velocity
 315                 print 'Loading rate (UNTESTED):'+str(lrate)
 316             except:
 317                 print 'Loading rate : n/a'
 318                 lrate='n/a'
 319
 320             if justone:
 321                 return
 322
 323             #save accept/discard/repeat
 324             print '\n\nDo you want to save these?'
 325             if self.YNclick():
 326
 327                 #Save file info
 328                 if self.autofile=='':
 329                     self.autofile=raw_input('Results filename? (press enter for a random generated name)')
 330                     if self.autofile=='':
 331                         [osf,name]=tempfile.mkstemp(prefix='analysis-',suffix='.txt',text=True,dir=self.config['hookedir'])
 332                         print 'saving in '+name
 333                         self.autofile=name
 334                         os.close(osf)
 335                         f=open(self.autofile,'a+')
 336                         f.write('Analysis started '+time.asctime()+'\n')
 337                         f.write('----------------------------------------\n')
 338                         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')
 339                         f.close()
 340
 341                 if not os.path.exists(self.autofile):
 342                     f=open(self.autofile,'a+')
 343                     f.write('Analysis started '+time.asctime()+'\n')
 344                     f.write('----------------------------------------\n')
 345                     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')
 346                     f.close()
 347
 348                 print 'Saving...'
 349                 savecounter+=1
 350                 f=open(self.autofile,'a+')
 351                 f.write(self.current.path)
 352                 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')
 353                 f.close()
 354             else:
 355                 print '\nWould you like to try again on this same curve?'
 356                 if self.YNclick():
 357                     continue
 358             curveindex+=1
 359
 360         if not justone:
 361             print 'You saved '+str(savecounter)+' out of '+str(len(self.current_list))+' total curves.'
 362
 363
 364
 365     def YNclick(self):
 366         #shows something in the graph to indicate we expect YNclick
 367         #FIXME is there an easy way to write in the graph?
 368         askplot=self._get_displayed_plot()
 369         center=np.min(askplot.vectors[1][0])+15e-9
 370         scale=np.std(askplot.vectors[1][1][-50:-1])
 371         hscale=np.mean(np.diff(askplot.vectors[1][0][-10:-1]))
 372         xhline=[center-20*hscale,center-10*hscale,center+10*hscale,center+20*hscale]
 373         ytopline=[15*scale,20*scale,20*scale,15*scale]
 374         ybotline=[-15*scale,-20*scale,-20*scale,-15*scale]
 375         yvline=np.arange(-25*scale,26*scale,scale/2)
 376         xvline=np.ones_like(yvline)*center
 377         xshow=np.concatenate((xvline,xhline,xhline))
 378         yshow=np.concatenate((yvline,ytopline,ybotline))
 379         askplot.add_set(xshow,yshow)
 380         askplot.styles+=['scatter']
 381         askplot.colors+=["black"]
 382         self._send_plot([askplot])
 383         print 'Click any point over y=0 for Yes, under it for No'
 384         point=self._measure_N_points(N=1,whatset=1)[0]
 385         value=point.absolute_coords[1]
 386         askplot.remove_set(-1)
 387         self._send_plot([askplot])
 388         if value<0:
 389             return False
 390         else:
 391             return True
 392
 393
 394
 395