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         forces=[]
 131         slopes=[]
 132
 133         #pick up plot
 134         displayed_plot=self._get_displayed_plot(0)
 135
 136         #COMMAND LINE PARSING
 137         #--Using points instead of nm interval
 138         if 'usepoints' in args.split():
 139             fit_points=int(self.config['auto_fit_points'])
 140             usepoints=True
 141         else:
 142             fit_points=None
 143             usepoints=False
 144         #--Recalculate baseline
 145         if 'rebase' in args or (self.basecurrent != self.current.path):
 146             rebase=True
 147
 148         #--Custom persistent length / custom temperature
 149         for arg in args.split():
 150             #look for a persistent length argument.
 151             if 'pl=' in arg:
 152                 pl_expression=arg.split('=')
 153                 pl_value=float(pl_expression[1]) #actual value
 154             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 155             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 156                 t_expression=arg.split('=')
 157                 T=float(t_expression[1])
 158         #--Contact point arguments
 159         if 'reclick' in args.split():
 160             print 'Click contact point'
 161             contact_point, contact_point_index = pickup_contact_point()
 162         elif 'noauto' in args.split():
 163             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 164                 print 'Click contact point'
 165                 contact_point , contact_point_index = pickup_contact_point()
 166             else:
 167                 contact_point=self.wlccontact_point
 168                 contact_point_index=self.wlccontact_index
 169         else:
 170             #Automatically find contact point
 171             cindex=self.find_contact_point()
 172             contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
 173         #--END COMMAND LINE PARSING--
 174
 175
 176         peak_location, peak_size = find_current_peaks()
 177
 178         fitplot=copy.deepcopy(displayed_plot)
 179
 180         #Pick up force baseline
 181         if rebase:
 182             clicks=self.config['baseline_clicks']
 183             if clicks==0:
 184                 self.basepoints=[]
 185                 base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
 186                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
 187                 base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
 188                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 189             elif clicks>0:
 190                 print 'Select baseline'
 191                 if clicks==1:
 192                     self.basepoints=self._measure_N_points(N=1, whatset=whatset)
 193                     base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
 194                     self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 195                 else:
 196                     self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 197
 198             self.basecurrent=self.current.path
 199
 200         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 201         boundaries.sort()
 202         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 203         avg=np.mean(to_average)
 204
 205
 206         for peak in peak_location:
 207             #WLC FITTING
 208             #define fit interval
 209             if not usepoints:
 210                 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
 211             peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
 212             other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
 213
 214             #points for the fit
 215             points=[contact_point, peak_point, other_fit_point]
 216
 217             if abs(peak_point.index-other_fit_point.index) < 2:
 218                 continue
 219
 220             params, yfit, xfit, fit_errors = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], pl_value, T)
 221
 222
 223             #Measure forces
 224             delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
 225             y=min(delta_to_measure)
 226             #save force values (pN)
 227             #Measure slopes
 228             slope=self.linefit_between(peak-slope_span,peak)[0]
 229
 230
 231             #check fitted data and, if right, add peak to the measurement
 232             #FIXME: code duplication
 233             if len(params)==1: #if we did choose 1-value fit
 234                 p_lengths.append(pl_value)
 235                 c_lengths.append(params[0]*(1.0e+9))
 236                 forces.append(abs(y-avg)*(1.0e+12))
 237                 slopes.append(slope)
 238                 #Add WLC fit lines to plot
 239                 fitplot.add_set(xfit,yfit)
 240                 if len(fitplot.styles)==0:
 241                     fitplot.styles=[]
 242                 else:
 243                     fitplot.styles.append(None)
 244             else:
 245                 p_leng=params[1]*(1.0e+9)
 246                 #check if persistent length makes sense. otherwise, discard peak.
 247                 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
 248                     p_lengths.append(p_leng)
 249                     c_lengths.append(params[0]*(1.0e+9))
 250                     forces.append(abs(y-avg)*(1.0e+12))
 251                     slopes.append(slope)
 252
 253                     #Add WLC fit lines to plot
 254                     fitplot.add_set(xfit,yfit)
 255                     if len(fitplot.styles)==0:
 256                         fitplot.styles=[]
 257                     else:
 258                         fitplot.styles.append(None)
 259                 else:
 260                     pass
 261
 262
 263         #add basepoints to fitplot
 264         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]])
 265         fitplot.styles.append('scatter')
 266
 267
 268         #Show wlc fits and peak locations
 269         self._send_plot([fitplot])
 270         #self.do_peaks('')
 271
 272
 273         #Ask the user what peaks to ignore from analysis.
 274         print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 275         print 'N to discard measurement'
 276         exclude_raw=raw_input('Input:')
 277         if exclude_raw=='N':
 278             print 'Discarded.'
 279             return
 280         if not exclude_raw=='':
 281             exclude=exclude_raw.split(',')
 282             try:
 283                 exclude=[int(item) for item in exclude]
 284                 for i in exclude:
 285                     c_lengths[i]=None
 286                     p_lengths[i]=None
 287                     forces[i]=None
 288                     slopes[i]=None
 289             except:
 290                  print 'Bad input, taking all...'
 291         #Clean data vectors from ignored peaks
 292         c_lengths=[item for item in c_lengths if item != None]
 293         p_lengths=[item for item in p_lengths if item != None]
 294         forces=[item for item in forces if item != None]
 295         slopes=[item for item in slopes if item != None]
 296         print 'contour (nm)',c_lengths
 297         print 'p (nm)',p_lengths
 298         print 'forces (pN)',forces
 299         print 'slopes (N/m)',slopes
 300
 301         #Save file info
 302         if self.autofile=='':
 303             self.autofile=raw_input('Autopeak filename? (return to ignore) ')
 304             if self.autofile=='':
 305                 print 'Not saved.'
 306                 return
 307
 308         if not os.path.exists(self.autofile):
 309             f=open(self.autofile,'w+')
 310             f.write('Analysis started '+time.asctime()+'\n')
 311             f.write('----------------------------------------\n')
 312             f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) \n')
 313             f.close()
 314
 315         print 'Saving...'
 316         f=open(self.autofile,'a+')
 317
 318         f.write(self.current.path+'\n')
 319         for i in range(len(c_lengths)):
 320             f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+'\n')
 321
 322         f.close()
 323         self.do_note('autopeak')
 324