hooke/plugin/autopeak.py

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