autopeak.py

   1 #!/usr/bin/env python
   2
   3 from libhooke import WX_GOOD, ClickedPoint
   4 import wxversion
   5 wxversion.select(WX_GOOD)
   6 from wx import PostEvent
   7 import numpy as np
   8 import scipy as sp
   9 import copy
  10 import os.path
  11 import time
  12
  13 import warnings
  14 warnings.simplefilter('ignore',np.RankWarning)
  15
  16
  17 class autopeakCommands:
  18
  19     def do_autopeak(self,args):
  20         '''
  21         AUTOPEAK
  22         (autopeak.py)
  23         Automatically performs a number of analyses on the peaks of the given curve.
  24         Currently it automatically:
  25         - fits peaks with WLC function
  26         - measures peak maximum forces with a baseline
  27         - measures slope in proximity of peak maximum
  28         Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
  29
  30         Syntax:
  31         autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
  32
  33         rebase : Re-asks baseline interval
  34
  35         pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
  36                      the fit will be a 2-variable
  37                      fit. DO NOT put spaces between 'pl', '=' and the value.
  38                      The value must be in meters.
  39                      Scientific notation like 0.35e-9 is fine.
  40
  41         t=[value] : Use a user-defined temperature. The value must be in
  42                     kelvins; by default it is 293 K.
  43                     DO NOT put spaces between 't', '=' and the value.
  44
  45         noauto : allows for clicking the contact point by
  46                  hand (otherwise it is automatically estimated) the first time.
  47                  If subsequent measurements are made, the same contact point
  48                  clicked the first time is used
  49
  50         reclick : redefines by hand the contact point, if noauto has been used before
  51                   but the user is unsatisfied of the previously choosen contact point.
  52
  53         usepoints : fit interval by number of points instead than by nanometers
  54
  55         When you first issue the command, it will ask for the filename. If you are giving the filename
  56         of an existing file, autopeak will resume it and append measurements to it. If you are giving
  57         a new filename, it will create the file and append to it until you close Hooke.
  58
  59
  60         Useful variables (to set with SET command):
  61         ---
  62         temperature= temperature of the system for wlc fit (in K)
  63
  64         auto_slope_span = number of points on which measure the slope, for slope
  65
  66         auto_fit_nm = number of nm to fit before the peak maximum, for WLC (if usepoints false)
  67         auto_fit_points = number of points to fit before the peak maximum, for WLC (if usepoints true)
  68
  69         baseline_clicks = 0: automatic baseline
  70                           1: decide baseline with a single click and length defined in auto_left_baseline
  71                           2: let user click points of baseline
  72         auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
  73         auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
  74         '''
  75
  76         #MACROS.
  77         #FIXME: to move outside function
  78         def fit_interval_nm(start_index,plot,nm,backwards):
  79             '''
  80             Calculates the number of points to fit, given a fit interval in nm
  81             start_index: index of point
  82             plot: plot to use
  83             backwards: if true, finds a point backwards.
  84             '''
  85             x_vect=plot.vectors[1][0]
  86
  87             c=0
  88             i=start_index
  89             start=x_vect[start_index]
  90             maxlen=len(x_vect)
  91             while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
  92                 if i==0 or i==maxlen-1: #we reached boundaries of vector!
  93                     return c
  94
  95                 if backwards:
  96                     i-=1
  97                 else:
  98                     i+=1
  99                 c+=1
 100             return c
 101
 102         def pickup_contact_point():
 103             '''macro to pick up the contact point by clicking'''
 104             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 105             contact_point_index=contact_point.index
 106             self.wlccontact_point=contact_point
 107             self.wlccontact_index=contact_point.index
 108             self.wlccurrent=self.current.path
 109             return contact_point, contact_point_index
 110
 111         def find_current_peaks():
 112             #Find peaks.
 113             defplot=self.current.curve.default_plots()[0]
 114             flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
 115             defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
 116             peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
 117             return peak_location, peak_size
 118
 119         #default fit etc. variables
 120         pl_value=None
 121         T=self.config['temperature']
 122
 123         slope_span=int(self.config['auto_slope_span'])
 124         delta_force=10
 125         rebase=False #if true=we select rebase
 126
 127         #initialize output data vectors
 128         c_lengths=[]
 129         p_lengths=[]
 130         sigma_c_lengths=[]
 131         sigma_p_lengths=[]
 132         forces=[]
 133         slopes=[]
 134
 135         #pick up plot
 136         displayed_plot=self._get_displayed_plot(0)
 137
 138         #COMMAND LINE PARSING
 139         #--Using points instead of nm interval
 140         if 'usepoints' in args.split():
 141             fit_points=int(self.config['auto_fit_points'])
 142             usepoints=True
 143         else:
 144             fit_points=None
 145             usepoints=False
 146         #--Recalculate baseline
 147         if 'rebase' in args or (self.basecurrent != self.current.path):
 148             rebase=True
 149
 150         #--Custom persistent length / custom temperature
 151         for arg in args.split():
 152             #look for a persistent length argument.
 153             if 'pl=' in arg:
 154                 pl_expression=arg.split('=')
 155                 pl_value=float(pl_expression[1]) #actual value
 156             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 157             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 158                 t_expression=arg.split('=')
 159                 T=float(t_expression[1])
 160         #--Contact point arguments
 161         if 'reclick' in args.split():
 162             print 'Click contact point'
 163             contact_point, contact_point_index = pickup_contact_point()
 164         elif 'noauto' in args.split():
 165             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 166                 print 'Click contact point'
 167                 contact_point , contact_point_index = pickup_contact_point()
 168             else:
 169                 contact_point=self.wlccontact_point
 170                 contact_point_index=self.wlccontact_index
 171         else:
 172             #Automatically find contact point
 173             cindex=self.find_contact_point()
 174             contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
 175         #--END COMMAND LINE PARSING--
 176
 177
 178         peak_location, peak_size = find_current_peaks()
 179
 180         fitplot=copy.deepcopy(displayed_plot)
 181
 182         #Pick up force baseline
 183         if rebase:
 184             clicks=self.config['baseline_clicks']
 185             if clicks==0:
 186                 self.basepoints=[]
 187                 base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
 188                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
 189                 base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
 190                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 191             elif clicks>0:
 192                 print 'Select baseline'
 193                 if clicks==1:
 194                     self.basepoints=self._measure_N_points(N=1, whatset=whatset)
 195                     base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
 196                     self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 197                 else:
 198                     self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 199
 200             self.basecurrent=self.current.path
 201
 202         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 203         boundaries.sort()
 204         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 205         avg=np.mean(to_average)
 206
 207
 208         for peak in peak_location:
 209             #WLC FITTING
 210             #define fit interval
 211             if not usepoints:
 212                 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
 213             peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
 214             other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
 215
 216             #points for the fit
 217             points=[contact_point, peak_point, other_fit_point]
 218
 219             if abs(peak_point.index-other_fit_point.index) < 2:
 220                 continue
 221
 222             params, yfit, xfit, fit_errors = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], pl_value, T, return_errors=True)
 223
 224
 225             #Measure forces
 226             delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
 227             y=min(delta_to_measure)
 228             #save force values (pN)
 229             #Measure slopes
 230             slope=self.linefit_between(peak-slope_span,peak)[0]
 231
 232
 233             #check fitted data and, if right, add peak to the measurement
 234             #FIXME: code duplication
 235             if len(params)==1: #if we did choose 1-value fit
 236                 p_lengths.append(pl_value)
 237                 c_lengths.append(params[0]*(1.0e+9))
 238                 sigma_p_lengths.append(0)
 239                 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 240                 forces.append(abs(y-avg)*(1.0e+12))
 241                 slopes.append(slope)
 242                 #Add WLC fit lines to plot
 243                 fitplot.add_set(xfit,yfit)
 244                 if len(fitplot.styles)==0:
 245                     fitplot.styles=[]
 246                 else:
 247                     fitplot.styles.append(None)
 248             else: #2-value fit
 249                 p_leng=params[1]*(1.0e+9)
 250                 #check if persistent length makes sense. otherwise, discard peak.
 251                 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
 252                     p_lengths.append(p_leng)
 253                     c_lengths.append(params[0]*(1.0e+9))
 254                     sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 255                     sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
 256                     forces.append(abs(y-avg)*(1.0e+12))
 257                     slopes.append(slope)
 258
 259                     #Add WLC fit lines to plot
 260                     fitplot.add_set(xfit,yfit)
 261                     if len(fitplot.styles)==0:
 262                         fitplot.styles=[]
 263                     else:
 264                         fitplot.styles.append(None)
 265                 else:
 266                     pass
 267
 268
 269         #add basepoints to fitplot
 270         fitplot.add_set([self.basepoints[0].graph_coords[0],self.basepoints[1].graph_coords[0]],[self.basepoints[0].graph_coords[1],self.basepoints[1].graph_coords[1]])
 271         fitplot.styles.append('scatter')
 272
 273
 274         #Show wlc fits and peak locations
 275         self._send_plot([fitplot])
 276         #self.do_peaks('')
 277
 278
 279         #Ask the user what peaks to ignore from analysis.
 280         print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 281         print 'N to discard measurement'
 282         exclude_raw=raw_input('Input:')
 283         if exclude_raw=='N':
 284             print 'Discarded.'
 285             return
 286         if not exclude_raw=='':
 287             exclude=exclude_raw.split(',')
 288             try:
 289                 exclude=[int(item) for item in exclude]
 290                 for i in exclude:
 291                     c_lengths[i]=None
 292                     p_lengths[i]=None
 293                     forces[i]=None
 294                     slopes[i]=None
 295                     sigma_c_lengths[i]=None
 296                     sigma_p_lengths[i]=None
 297             except:
 298                  print 'Bad input, taking all...'
 299         #Clean data vectors from ignored peaks
 300         #FIXME:code duplication
 301         c_lengths=[item for item in c_lengths if item != None]
 302         p_lengths=[item for item in p_lengths if item != None]
 303         forces=[item for item in forces if item != None]
 304         slopes=[item for item in slopes if item != None]
 305         sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
 306         sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
 307
 308         print 'contour (nm)',c_lengths
 309         print 'sigma contour (nm)',sigma_c_lengths
 310         print 'p (nm)',p_lengths
 311         print 'sigma p (nm)',sigma_p_lengths
 312         print 'forces (pN)',forces
 313         print 'slopes (N/m)',slopes
 314
 315         #Save file info
 316         if self.autofile=='':
 317             self.autofile=raw_input('Autopeak filename? (return to ignore) ')
 318             if self.autofile=='':
 319                 print 'Not saved.'
 320                 return
 321
 322         if not os.path.exists(self.autofile):
 323             f=open(self.autofile,'w+')
 324             f.write('Analysis started '+time.asctime()+'\n')
 325             f.write('----------------------------------------\n')
 326             f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) ;  Sigma contour (nm) ; Sigma persistence (nm)\n')
 327             f.close()
 328
 329         print 'Saving...'
 330         f=open(self.autofile,'a+')
 331
 332         f.write(self.current.path+'\n')
 333         for i in range(len(c_lengths)):
 334             f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+' ; '+str(sigma_c_lengths[i])+' ; '+str(sigma_p_lengths[i])+'\n')
 335
 336         f.close()
 337         self.do_note('autopeak')
 338