hooke/plugin/multifit.py

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