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. Obtains
  40         contour length, persistance length, rupture force and slope - loading rate.
  41         WLC is shown in red, FJC in blue.
  42         -------------
  43         Syntax:
  44         multifit [pl=value] [kl=value] [t=value] [slopew=value] [basew=value] [justone]
  45
  46                 pl=[value] and kl=[value]: Use a fixed persistent length (WLC) or Kuhn length (FJC) for
  47                         the fit. If pl is not given, the fit will be a 2-variable fit.
  48                         DO NOT put spaces between 'pl', '=' and the value.
  49         The value must be in nanometers.
  50
  51         t=[value] : Use a user-defined temperature. The value must be in kelvins;
  52             by default it is 293 K.
  53             DO NOT put spaces between 't', '=' and the value.
  54
  55         slopew and basew : width in points for slope fitting (points to the right of clicked rupture)
  56             and base level fitting(points to the left of clicked top of rupture), default is 15.
  57             DO NOT put spaces between 'slopew' or 'basew', '=' and the value.
  58
  59         justone : performs the fits over current curve instead of iterating
  60
  61         see fit command help for more information on the options and fit procedures.
  62
  63                 NOTE: centerzero plot modifier should be activated (set centerzero 1).
  64                 '''
  65
  66                 #NOTE duplicates a lot of code from do_fit in fit.py, a call to it could be used directly
  67                 #but it is easier to control the program flow bypassing it
  68
  69         pl_value=None
  70         kl_value=None
  71         T=self.config['temperature']
  72         slopew=15
  73         basew=15
  74         justone=False
  75
  76         #FIXME spaces are not allowed between pl or t and value
  77         for arg in args.split():
  78             #look for a persistent length argument.
  79             if 'pl=' in arg:
  80                 pl_expression=arg.split('=')
  81                 pl_value=float(pl_expression[1]) #actual value
  82             if 'kl=' in arg:
  83                 kl_expression=arg.split('=')
  84                 kl_value=float(kl_expression[1]) #actual value
  85             #look for a T argument.
  86             if ('t=' in arg) or ('T=' in arg):
  87                 t_expression=arg.split('=')
  88                 T=float(t_expression[1])
  89             #look for a basew argument.
  90             if ('basew=' in arg):
  91                 basew_expression=arg.split('=')
  92                 basew=int(basew_expression[1])
  93             #look for a slopew argument.
  94             if ('slopew=' in arg):
  95                 slopew_expression=arg.split('=')
  96                 slopew=int(slopew_expression[1])
  97             if('justone' in arg):
  98                 justone=True
  99
 100         counter=0
 101         savecounter=0
 102         curveindex=0
 103
 104         if not justone:
 105             print 'What curve no. you would like to start? (enter for ignoring)'
 106             skip=raw_input()
 107
 108             if skip.isdigit()==False:
 109                 skip=0
 110             else:
 111                 skip=int(skip)-1
 112                 print 'Skipping '+str(skip)+ ' curves'
 113         else:
 114             skip=0
 115         #begin presenting curves for analysis
 116         while curveindex <len(self.current_list):
 117             if not justone:
 118                 counter+=1
 119                 curve=self.current_list[curveindex]
 120                 if skip>curveindex:
 121                     curveindex+=1
 122                     continue
 123
 124             #give periodically the opportunity to stop the analysis
 125                 if counter%20==0 and counter>0:
 126                     print '\nYou checked '+str(counter)+' curves. Do you want to continue?'
 127                     self.current=curve
 128                     self.do_plot(0)
 129                     if self.YNclick():
 130                         print 'Going on...'
 131                     else:
 132                         break
 133                 else:
 134                     self.current=curve
 135                     self.do_plot(0)
 136             else:
 137                 curve=self.current
 138                 self.do_plot(0)
 139             if not justone:
 140                 print '\nCurve '+str(curveindex+1)+' of '+str(len(self.current_list))
 141             print 'Click contact point or left end of the curve to skip'
 142             #FIXME "left half" is a bit ad hoc, and "set correct 1" makes
 143             #the 3/4s point not very reliable.
 144             #Anyway clicking the end should be safe
 145             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 146             contact_point_index=contact_point.index
 147
 148             retract=self.plots[0].vectors[1][0]
 149
 150             #some fixing for x data that is negative (depends on driver)
 151             if min(retract)<0:
 152                 cppoint=contact_point.graph_coords[0]+abs(min(retract))
 153                 retract=retract+abs(min(retract))
 154             else:
 155                 cppoint=contact_point.graph_coords[0]
 156             threequarters=3*(max(retract)-min(retract))/4
 157
 158
 159             if cppoint < threequarters:
 160                 if justone:
 161                     break
 162                 curveindex+=1
 163                 continue
 164
 165             self.wlccontact_point=contact_point
 166             self.wlccontact_index=contact_point.index
 167             self.wlccurrent=self.current.path
 168
 169             print 'Now click two points for the fitting area (one should be the rupture point)'
 170             wlcpoints=self._measure_N_points(N=2,whatset=1)
 171             print 'And one point of the top of the jump'
 172             toppoint=self._measure_N_points(N=1,whatset=1)
 173
 174             fitpoints=[contact_point]+wlcpoints
 175             #use the currently displayed plot for the fit
 176             displayed_plot=self._get_displayed_plot()
 177
 178             #use both fit functions
 179             try:
 180                 wlcparams, wlcyfit, wlcxfit, wlcfit_errors = self.wlc_fit(fitpoints, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T, return_errors=True )
 181                 wlcerror=False
 182             except:
 183                 print 'WLC fit not possible'
 184                 wlcerror=True
 185
 186             try:
 187                 fjcparams, fjcyfit, fjcxfit, fjcfit_errors = self.fjc_fit(fitpoints, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],kl_value,T, return_errors=True )
 188                 fjcerror=False
 189             except:
 190                 print 'FJC fit not possible'
 191                 fjcerror=True
 192
 193             #Measure rupture force
 194             ruptpoint=ClickedPoint()
 195             if wlcpoints[0].graph_coords[1]<wlcpoints[1].graph_coords[1]:
 196                 ruptpoint=wlcpoints[0]
 197             else:
 198                 ruptpoint=wlcpoints[1]
 199             tpindex=toppoint[0].index
 200             toplevel=np.average(displayed_plot.vectors[1][1][tpindex:tpindex+basew])
 201             force=toplevel-ruptpoint.graph_coords[1]
 202
 203             #Measure the slope - loading rate
 204             slope=self._slope([ruptpoint],slopew)
 205
 206             #plot results (_slope already did)
 207
 208             #now we have the fit, we can plot it
 209             #add the clicked points in the final PlotObject
 210             clickvector_x, clickvector_y=[], []
 211             for item in wlcpoints:
 212                 clickvector_x.append(item.graph_coords[0])
 213                 clickvector_y.append(item.graph_coords[1])
 214
 215             #create a custom PlotObject to gracefully plot the fit along the curves
 216
 217             #first fix those irritating zoom-destroying fits
 218             lowestpoint=min(displayed_plot.vectors[1][1])
 219             for i in np.arange(0,len(wlcyfit)):
 220                 if wlcyfit[i] < 1.2*lowestpoint:
 221                     wlcyfit[i]=toplevel
 222             for i in np.arange(0,len(fjcyfit)):
 223                 if fjcyfit[i] < 1.2*lowestpoint:
 224                     fjcyfit[i]=toplevel
 225
 226             fitplot=copy.deepcopy(displayed_plot)
 227             fitplot.add_set(wlcxfit,wlcyfit)
 228             fitplot.add_set(fjcxfit,fjcyfit)
 229             fitplot.add_set(clickvector_x,clickvector_y)
 230
 231             fitplot.styles+=[None,None,'scatter']
 232             fitplot.colors+=["red","blue",None]
 233
 234             self._send_plot([fitplot])
 235
 236
 237              #present results of measurement
 238             if len(wlcparams)==1:
 239                 wlcparams.append(pl_value*1e-9)
 240             if len(fjcparams)==1:
 241                 fjcparams.append(kl_value*1e-9)
 242
 243             if fjcfit_errors:
 244                 fjcfit_nm=[i*(10**9) for i in fjcfit_errors]
 245                 if len(fjcfit_nm)==1:
 246                     fjcfit_nm.append(0)
 247             else:
 248                 fjcfit_errors=[0,0]
 249             if wlcfit_errors:
 250                 wlcfit_nm=[i*(10**9) for i in wlcfit_errors]
 251                 if len(wlcfit_nm)==1:
 252                     wlcfit_nm.append(0)
 253             else:
 254                 wlcfit_errors=[0,0]
 255
 256
 257             print '\nRESULTS'
 258             print 'WLC contour : '+str(1e9*wlcparams[0])+u' \u00b1 '+str(wlcfit_nm[0])+' nm'
 259             print 'Per. length : '+str(1e9*wlcparams[1])+u' \u00b1 '+str(wlcfit_nm[1])+' nm'
 260             print '---'
 261             print 'FJC contour : '+str(1e9*fjcparams[0])+u' \u00b1 '+str(fjcfit_nm[0])+' nm'
 262             print 'Kuhn length : '+str(1e9*fjcparams[1])+u' \u00b1 '+str(fjcfit_nm[1])+' nm'
 263             print '---'
 264             print 'Force : '+str(1e12*force)+' pN'
 265             print 'Slope : '+str(slope)+' N/m'
 266             try:
 267                 #FIXME all drivers should provide retract velocity, in SI or homogeneous units
 268                 lrate=slope*self.current.curve.retract_velocity
 269                 print 'Loading rate (UNTESTED):'+str(lrate)
 270             except:
 271                 print 'Loading rate : n/a'
 272                 lrate='n/a'
 273
 274             if justone:
 275                 return
 276
 277             #save accept/discard/repeat
 278             print '\n\nDo you want to save these?'
 279             if self.YNclick():
 280
 281                 #Save file info
 282                 if self.autofile=='':
 283                     self.autofile=raw_input('Results filename? (press enter for a random generated name)')
 284                     if self.autofile=='':
 285                         [osf,name]=tempfile.mkstemp(prefix='analysis-',suffix='.txt',text=True,dir=self.config['hookedir'])
 286                         print 'saving in '+name
 287                         self.autofile=name
 288                         os.close(osf)
 289                         f=open(self.autofile,'a+')
 290                         f.write('Analysis started '+time.asctime()+'\n')
 291                         f.write('----------------------------------------\n')
 292                         f.write(' File ; WLC Cont. length (nm) ; Sigma WLC cl;  Per. Length (nm) ; Sigma pl; FJC Cont. length (nm) ; Sigma FJC cl ; Kuhn length (nm); Sigma kl ; Force (pN) ; Slope (N/m) ; (untested) Loading rate (pN/s)\n')
 293                         f.close()
 294
 295                 if not os.path.exists(self.autofile):
 296                     f=open(self.autofile,'a+')
 297                     f.write('Analysis started '+time.asctime()+'\n')
 298                     f.write('----------------------------------------\n')
 299                     f.write(' File ; WLC Cont. length (nm) ; Sigma WLC cl;  Per. Length (nm) ; Sigma pl; FJC Cont. length (nm) ; Sigma FJC cl ; Kuhn length (nm); Sigma kl ; Force (pN) ; Slope (N/m) ; (untested) Loading rate (pN/s)\n')
 300                     f.close()
 301
 302                 print 'Saving...'
 303                 savecountr+=1
 304                 f=open(self.autofile,'a+')
 305                 f.write(self.current.path)
 306                 f.write(' ; '+str(1e9*wlcparams[0])+' ; '+str(wlcfit_nm[0])+' ; '+str(1e9*wlcparams[1])+' ; '+str(wlcfit_nm[1])+' ; '+str(1e9*fjcparams[0])+' ; '+str(fjcfit_nm[0])+' ; '+str(1e9*fjcparams[1])+' ; '+str(fjcfit_nm[1])+' ; '+str(1e12*force)+' ; '+ str(slope)+' ; '+str(lrate)+'\n')
 307                 f.close()
 308             else:
 309                 print '\nWould you like to try again on this same curve?'
 310                 if self.YNclick():
 311                     continue
 312             curveindex+=1
 313
 314         if not justone:
 315             print 'You saved '+str(savecounter)+' out of '+str(len(self.current_list))+' total curves.'
 316
 317
 318
 319     def YNclick(self):
 320         print 'Click any point over y=0 for Yes, under it for No'
 321         point=self._measure_N_points(N=1,whatset=1)[0]
 322         value=point.absolute_coords[1]
 323         if value<0:
 324             return False
 325         else:
 326             return True
 327
 328
 329
 330