hooke/plugin/jumpstat.py

   1 # Copyright (C) 2010-2012 W. Trevor King <wking@drexel.edu>
   2 #
   3 # This file is part of Hooke.
   4 #
   5 # Hooke is free software: you can redistribute it and/or modify it
   6 # under the terms of the GNU Lesser General Public License as
   7 # published by the Free Software Foundation, either version 3 of the
   8 # License, or (at your option) any later version.
   9 #
  10 # Hooke is distributed in the hope that it will be useful, but WITHOUT
  11 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12 # or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  13 # Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with Hooke.  If not, see
  17 # <http://www.gnu.org/licenses/>.
  18
  19 from ..libhooke import WX_GOOD, ClickedPoint
  20 import wxversion
  21 wxversion.select(WX_GOOD)
  22 from wx import PostEvent
  23 import numpy as np
  24 import scipy as sp
  25 import copy
  26 import os.path
  27 import time
  28 import sys
  29 import warnings
  30 warnings.simplefilter('ignore',np.RankWarning)
  31
  32
  33 class jumpstatCommands():
  34
  35     def do_jumpstat(self,args):
  36      '''
  37      JUMPSTAT
  38      jumpstat.py
  39      Based on the convolution recognition automatically give:
  40      - the delta distance between the peaks,
  41      - the delta-force from the top of the peaks and subsequent relaxation,
  42      - the delta-force from the top of the peaks and the baseline
  43      The command allow also to remove the unwanted peaks that can be due to interference.
  44      When you first issue the command, it will ask for the filename. If you are giving the filename
  45      of an existing file, autopeak will resume it and append measurements to it. If you are giving
  46      a new filename, it will create the file and append to it until you close Hooke.
  47      You can also define a minimun deviation of the peaks.
  48
  49      Syntax:
  50      jumpstat [deviation]
  51      deviation = number of times the convolution signal is above the noise absolute deviation.
  52      '''
  53
  54
  55      #finding the max and the minimum positions for all the peaks
  56      noflatten=False
  57      #we use if else only to avoid a "bad input" message from find_current_peaks
  58      if (len(args)==0):
  59              max_peaks_location, peak_size=self.find_current_peaks(noflatten)
  60              min_peaks_location, pks2=self.find_current_peaks(noflatten, True, False)
  61      else:
  62              max_peaks_location, peak_size=self.find_current_peaks(noflatten, args)
  63              min_peaks_location, pks2=self.find_current_peaks(noflatten, args, False)
  64
  65
  66      #print "max_peaks_location:  "+str(len(max_peaks_location))
  67      #print "min_peaks_location:  "+str(len(min_peaks_location))
  68
  69      #if no peaks, we have nothing to plot. exit.
  70      if len(max_peaks_location)==0:
  71             print "No peaks on this curve."
  72             return
  73
  74      if len(max_peaks_location)!=len(min_peaks_location):
  75             print "Something went wrong in peaks recognition, number of minima is different from number of maxima. Exiting."
  76             return
  77
  78      #otherwise, we plot the peak locations.
  79      xplotted_ret=self.plots[0].vectors[1][0]
  80      yplotted_ret=self.plots[0].vectors[1][1]
  81      xgood=[xplotted_ret[index] for index in max_peaks_location]
  82      ygood=[yplotted_ret[index] for index in max_peaks_location]
  83
  84      xafter=[xplotted_ret[index] for index in min_peaks_location]
  85      yafter=[yplotted_ret[index] for index in min_peaks_location]
  86
  87      recplot=self._get_displayed_plot()
  88      recplot2=self._get_displayed_plot()
  89      recplot.vectors.append([xgood,ygood])
  90      recplot2.vectors.append([xafter,yafter])
  91
  92      if recplot.styles==[]:
  93          recplot.styles=[None,None,'scatter']
  94          recplot.colors=[None,None,None]
  95      else:
  96          recplot.styles+=['scatter']
  97          recplot.colors+=[None]
  98
  99      if recplot2.styles==[]:
 100          recplot2.styles=[None,None,None]
 101          recplot2.colors=[None,'1.0',None]
 102      else:
 103          recplot2.styles+=['scatter']
 104          recplot2.colors+=['0.5']
 105
 106      self._send_plot([recplot])
 107      self._send_plot([recplot2])
 108
 109
 110      #finding the baseline
 111      self.basepoints=self.baseline_points(max_peaks_location, recplot)
 112      boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 113      boundaries.sort()
 114      to_average=recplot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 115      avg=np.mean(to_average)
 116
 117
 118      dist=[]
 119      jumpforce=[]
 120      force=[]
 121
 122      #we calculate the distance vector
 123      for g in range(len(max_peaks_location)-1):
 124        dist.append((10**9)*(xplotted_ret[max_peaks_location[g]]-xplotted_ret[max_peaks_location[g+1]]))
 125      print "Distance values for the peaks in nm:"
 126      print dist
 127
 128      #the jump-force vector
 129      for g in range(len(max_peaks_location)):
 130       jumpforce.append((10**12)     *(yplotted_ret[min_peaks_location[g]] -yplotted_ret[max_peaks_location[g]])   )
 131      print "Force values for the jumps of the peaks in pN:"
 132      print jumpforce
 133
 134      #the force from baseline vector
 135      for g in range(len(max_peaks_location)):
 136       force.append((10**12)*(avg-yplotted_ret[max_peaks_location[g]]))
 137      print "Force values for the peaks in pN:"
 138      print force
 139
 140
 141
 142      #Now ask for the peaks that we don't want
 143      print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 144      print 'N to discard measurement'
 145      exclude_raw=raw_input('Input:')
 146      if exclude_raw=='N':
 147         print 'Discarded.'
 148         return
 149
 150      if not exclude_raw=='':
 151         exclude=exclude_raw.split(',')
 152         #we convert in numbers the input
 153         try:
 154             exclude=[int(item) for item in exclude]
 155         except:
 156             print 'Bad input, taking nothing.'
 157             return
 158
 159 #    we remove the peaks that we don't want from the list, we need a counter beacause if we remove
 160 #    a peaks the other peaks in the list are shifted by one at each step
 161         count=0
 162         for a in exclude:
 163           if (a==0):
 164              max_peaks_location=max_peaks_location[1:]
 165              min_peaks_location=min_peaks_location[1:]
 166           else:
 167              new_a=a-count
 168              max_peaks_location=  max_peaks_location[0:new_a]+max_peaks_location[new_a+1:]
 169              min_peaks_location=  min_peaks_location[0:new_a]+min_peaks_location[new_a+1:]
 170              peak_size=            peak_size[0:new_a]+peak_size[new_a+1:]
 171           count+=1
 172
 173
 174      #print "max_peaks_location:  "+str(len(max_peaks_location))
 175      #print "min_peaks_location:  "+str(len(min_peaks_location))
 176
 177
 178      dist=[]
 179      jumpforce=[]
 180      force=[]
 181      #we recalculate the distances and the forces after the removing of the unwanted peaks
 182      for g in range(len(max_peaks_location)-1):
 183          dist.append(xplotted_ret[max_peaks_location[g]]-xplotted_ret[max_peaks_location[g+1]])
 184      for g in range(len(max_peaks_location)):
 185          jumpforce.append( yplotted_ret[min_peaks_location[g]] - yplotted_ret[max_peaks_location[g]]   )
 186      for g in range(len(max_peaks_location)):
 187          force.append(avg  -  yplotted_ret[max_peaks_location[g]])
 188
 189
 190
 191
 192
 193         #Save file info
 194      if self.autofile=='':
 195             self.autofile=raw_input('Jumpstat filename? (return to ignore) ')
 196             if self.autofile=='':
 197                 print 'Not saved.'
 198                 return
 199
 200      if not os.path.exists(self.autofile):
 201             f=open(self.autofile,'w+')
 202             f.write('Analysis started '+time.asctime()+'\n')
 203             f.write('----------------------------------------\n')
 204             f.write('; Delta Distance length (m);  Jump Force pN;  Standard Force pN\n')
 205             f.write(self.current.path+'\n')
 206             for k in range(len(dist)):
 207                f.write(";")
 208                f.write(str(dist[k])+";"+str(jumpforce[k])+";"+str(force[k])+"\n"   )
 209             f.write("\n")
 210             f.close()
 211
 212      else:
 213             f=open(self.autofile,'a+')
 214             f.write(self.current.path+'\n')
 215             for k in range(len(dist)):
 216               f.write(";")
 217               f.write(str(dist[k])+";"+str(jumpforce[k])+";"+str(force[k])+"\n"   )
 218             f.write("\n")
 219             f.close()
 220
 221      print 'Saving...'