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