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