hooke/plugin/multifit.py

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