From: albertogomcas <devnull@localhost>
Date: Thu, 11 Mar 2010 12:02:18 +0000 (+0000)
Subject: multifit.py :	   new plugin, mouse input based
X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=d2b868bb350d1ddf8c5d46343c74fdd0114060d1;p=hooke.git

multifit.py :	   new plugin, mouse input based
generalvclamp.py : separated code from do_slope into _slope allowing use in other funtions
hooke.conf:        activate multifit plugin by default
---

diff --git a/generalvclamp.py b/generalvclamp.py
index 85a4bce..3e534a9 100644
--- a/generalvclamp.py
+++ b/generalvclamp.py
@@ -253,30 +253,34 @@ class generalvclampCommands:
         if fitspan == 0:
             # Gets the Xs of two clicked points as indexes on the current curve vector
             print 'Click twice to delimit chunk'
-            clickedpoints=[]
             points=self._measure_N_points(N=2,whatset=1)
-            clickedpoints=[points[0].index,points[1].index]
-            clickedpoints.sort()
         else:
             print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
-            clickedpoints=[]
             points=self._measure_N_points(N=1,whatset=1)
-            clickedpoints=[points[0].index-fitspan,points[0].index]
-        
-        # Calls the function linefit_between
+            
+        slope=self._slope(points,fitspan)
+
+        # Outputs the relevant slope parameter
+        print 'Slope:'
+        print str(slope)
+        to_dump='slope '+self.current.path+' '+str(slope)
+        self.outlet.push(to_dump)
+
+    def _slope(self,points,fitspan):
+		# Calls the function linefit_between
         parameters=[0,0,[],[]]
+        try:
+            clickedpoints=[points[0].index,points[1].index]
+            clickedpoints.sort()
+        except:
+            clickedpoints=[points[0].index-fitspan,points[0].index]        
+
         try:
             parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
         except:
             print 'Cannot fit. Did you click twice the same point?'
             return
-          
-        # Outputs the relevant slope parameter
-        print 'Slope:'
-        print str(parameters[0])
-        to_dump='slope '+self.current.path+' '+str(parameters[0])
-        self.outlet.push(to_dump)
-                
+             
         # Makes a vector with the fitted parameters and sends it to the GUI
         xtoplot=parameters[2]
         ytoplot=[]
@@ -307,6 +311,9 @@ class generalvclampCommands:
         
         self._send_plot([lineplot])
 
+        return parameters[0]	
+
+
     def linefit_between(self,index1,index2,whatset=1):
         '''
         Creates two vectors (xtofit,ytofit) slicing out from the
diff --git a/hooke.conf b/hooke.conf
index 6723440..3a49693 100755
--- a/hooke.conf
+++ b/hooke.conf
@@ -39,7 +39,8 @@ This section defines which plugins have to be loaded by Hooke.
     <generaltccd/>
     <multidistance/>
     <jumpstat/>
-    <review/>	
+    <review/>
+    <multifit/>		
 </plugins>
 
 <!--
diff --git a/multifit.py b/multifit.py
new file mode 100644
index 0000000..3f66641
--- /dev/null
+++ b/multifit.py
@@ -0,0 +1,330 @@
+#!/usr/bin/env python
+
+'''
+multifit.py
+
+Alberto Gomez-Casado, (c) 2010, University of Twente (The Netherlands)
+Licensed under GNU GPL v2
+'''
+
+#FIXME clean this, probably some dependencies are not needed
+
+from libhooke import WX_GOOD, ClickedPoint
+import wxversion
+wxversion.select(WX_GOOD)
+from wx import PostEvent
+import numpy as np
+import scipy as sp
+import copy
+import os.path
+import time
+import tempfile
+import warnings
+warnings.simplefilter('ignore',np.RankWarning)
+import Queue
+
+global measure_wlc
+global EVT_MEASURE_WLC
+
+global events_from_fit
+events_from_fit=Queue.Queue() #GUI ---> CLI COMMUNICATION
+
+
+class multifitCommands:
+
+    def do_multifit(self,args):
+        '''
+        MULTIFIT
+        (multifit.py)
+        Presents curves for manual analysis in a comfortable mouse-only fashion. Obtains
+        contour length, persistance length, rupture force and slope - loading rate.
+        WLC is shown in red, FJC in blue.
+        -------------
+        Syntax:
+        multifit [pl=value] [kl=value] [t=value] [slopew=value] [basew=value] [justone]
+
+		pl=[value] and kl=[value]: Use a fixed persistent length (WLC) or Kuhn length (FJC) for
+			the fit. If pl is not given, the fit will be a 2-variable fit. 
+			DO NOT put spaces between 'pl', '=' and the value.
+        The value must be in nanometers. 
+        
+        t=[value] : Use a user-defined temperature. The value must be in kelvins; 
+            by default it is 293 K.
+            DO NOT put spaces between 't', '=' and the value.
+        
+        slopew and basew : width in points for slope fitting (points to the right of clicked rupture) 
+            and base level fitting(points to the left of clicked top of rupture), default is 15.
+            DO NOT put spaces between 'slopew' or 'basew', '=' and the value.
+        
+        justone : performs the fits over current curve instead of iterating
+
+        see fit command help for more information on the options and fit procedures.
+
+		NOTE: centerzero plot modifier should be activated (set centerzero 1).
+		'''
+
+		#NOTE duplicates a lot of code from do_fit in fit.py, a call to it could be used directly
+		#but it is easier to control the program flow bypassing it
+	
+        pl_value=None
+        kl_value=None
+        T=self.config['temperature']
+        slopew=15
+        basew=15
+        justone=False
+        
+        #FIXME spaces are not allowed between pl or t and value
+        for arg in args.split():
+            #look for a persistent length argument.
+            if 'pl=' in arg:
+                pl_expression=arg.split('=')
+                pl_value=float(pl_expression[1]) #actual value
+            if 'kl=' in arg:
+                kl_expression=arg.split('=')
+                kl_value=float(kl_expression[1]) #actual value
+            #look for a T argument.
+            if ('t=' in arg) or ('T=' in arg):
+                t_expression=arg.split('=')
+                T=float(t_expression[1])
+            #look for a basew argument.
+            if ('basew=' in arg):
+                basew_expression=arg.split('=')
+                basew=int(basew_expression[1])
+            #look for a slopew argument.
+            if ('slopew=' in arg):
+                slopew_expression=arg.split('=')
+                slopew=int(slopew_expression[1])
+            if('justone' in arg):
+                justone=True
+		
+        counter=0
+        savecounter=0
+        curveindex=0
+        
+        if not justone:
+            print 'What curve no. you would like to start? (enter for ignoring)'
+            skip=raw_input()
+
+            if skip.isdigit()==False:
+                skip=0
+            else:
+                skip=int(skip)-1
+                print 'Skipping '+str(skip)+ ' curves'
+        else:
+            skip=0
+        #begin presenting curves for analysis
+        while curveindex <len(self.current_list):
+            if not justone:
+                counter+=1
+                curve=self.current_list[curveindex]
+                if skip>curveindex:
+                    curveindex+=1
+                    continue	
+
+            #give periodically the opportunity to stop the analysis
+                if counter%20==0 and counter>0:
+                    print '\nYou checked '+str(counter)+' curves. Do you want to continue?'
+                    self.current=curve
+                    self.do_plot(0)
+                    if self.YNclick():
+                        print 'Going on...'
+                    else:
+                        break
+                else:
+                    self.current=curve
+                    self.do_plot(0)
+            else:
+                curve=self.current
+                self.do_plot(0)
+            if not justone:
+                print '\nCurve '+str(curveindex+1)+' of '+str(len(self.current_list))       
+            print 'Click contact point or left end of the curve to skip'
+            #FIXME "left half" is a bit ad hoc, and "set correct 1" makes
+            #the 3/4s point not very reliable. 
+            #Anyway clicking the end should be safe
+            contact_point=self._measure_N_points(N=1, whatset=1)[0]
+            contact_point_index=contact_point.index
+
+            retract=self.plots[0].vectors[1][0]
+            
+            #some fixing for x data that is negative (depends on driver)
+            if min(retract)<0:
+                cppoint=contact_point.graph_coords[0]+abs(min(retract))
+                retract=retract+abs(min(retract))
+            else:
+                cppoint=contact_point.graph_coords[0]
+            threequarters=3*(max(retract)-min(retract))/4
+                
+  
+            if cppoint < threequarters:
+                if justone:
+                    break
+                curveindex+=1
+                continue				
+                
+            self.wlccontact_point=contact_point
+            self.wlccontact_index=contact_point.index
+            self.wlccurrent=self.current.path
+                
+            print 'Now click two points for the fitting area (one should be the rupture point)'
+            wlcpoints=self._measure_N_points(N=2,whatset=1)
+            print 'And one point of the top of the jump'
+            toppoint=self._measure_N_points(N=1,whatset=1)
+
+            fitpoints=[contact_point]+wlcpoints
+            #use the currently displayed plot for the fit
+            displayed_plot=self._get_displayed_plot()
+
+            #use both fit functions
+            try:
+                wlcparams, wlcyfit, wlcxfit, wlcfit_errors = self.wlc_fit(fitpoints, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T, return_errors=True )
+                wlcerror=False	
+            except:
+                print 'WLC fit not possible'
+                wlcerror=True
+
+            try:
+                fjcparams, fjcyfit, fjcxfit, fjcfit_errors = self.fjc_fit(fitpoints, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],kl_value,T, return_errors=True )
+                fjcerror=False
+            except:
+                print 'FJC fit not possible'
+                fjcerror=True
+                
+            #Measure rupture force
+            ruptpoint=ClickedPoint()    
+            if wlcpoints[0].graph_coords[1]<wlcpoints[1].graph_coords[1]:
+                ruptpoint=wlcpoints[0]
+            else:
+                ruptpoint=wlcpoints[1]
+            tpindex=toppoint[0].index
+            toplevel=np.average(displayed_plot.vectors[1][1][tpindex:tpindex+basew])
+            force=toplevel-ruptpoint.graph_coords[1]					
+         
+            #Measure the slope - loading rate
+            slope=self._slope([ruptpoint],slopew)
+         
+            #plot results (_slope already did)
+            
+            #now we have the fit, we can plot it
+            #add the clicked points in the final PlotObject
+            clickvector_x, clickvector_y=[], []
+            for item in wlcpoints:
+                clickvector_x.append(item.graph_coords[0])
+                clickvector_y.append(item.graph_coords[1])
+            
+            #create a custom PlotObject to gracefully plot the fit along the curves
+            
+            #first fix those irritating zoom-destroying fits
+            lowestpoint=min(displayed_plot.vectors[1][1])
+            for i in np.arange(0,len(wlcyfit)):
+                if wlcyfit[i] < 1.2*lowestpoint:
+                    wlcyfit[i]=toplevel    
+            for i in np.arange(0,len(fjcyfit)):
+                if fjcyfit[i] < 1.2*lowestpoint:
+                    fjcyfit[i]=toplevel
+            
+            fitplot=copy.deepcopy(displayed_plot)
+            fitplot.add_set(wlcxfit,wlcyfit)
+            fitplot.add_set(fjcxfit,fjcyfit) 
+            fitplot.add_set(clickvector_x,clickvector_y)
+
+            fitplot.styles+=[None,None,'scatter']
+            fitplot.colors+=["red","blue",None]
+            
+            self._send_plot([fitplot])
+            
+
+             #present results of measurement
+            if len(wlcparams)==1:
+                wlcparams.append(pl_value*1e-9)
+            if len(fjcparams)==1:
+                fjcparams.append(kl_value*1e-9)
+                
+            if fjcfit_errors:
+                fjcfit_nm=[i*(10**9) for i in fjcfit_errors]
+                if len(fjcfit_nm)==1:
+                    fjcfit_nm.append(0)
+            else:
+                fjcfit_errors=[0,0]        
+            if wlcfit_errors:
+                wlcfit_nm=[i*(10**9) for i in wlcfit_errors]
+                if len(wlcfit_nm)==1:
+                    wlcfit_nm.append(0)
+            else:
+                wlcfit_errors=[0,0]
+            
+            
+            print '\nRESULTS'
+            print 'WLC contour : '+str(1e9*wlcparams[0])+u' \u00b1 '+str(wlcfit_nm[0])+' nm'
+            print 'Per. length : '+str(1e9*wlcparams[1])+u' \u00b1 '+str(wlcfit_nm[1])+' nm'
+            print '---'
+            print 'FJC contour : '+str(1e9*fjcparams[0])+u' \u00b1 '+str(fjcfit_nm[0])+' nm'
+            print 'Kuhn length : '+str(1e9*fjcparams[1])+u' \u00b1 '+str(fjcfit_nm[1])+' nm'    
+            print '---'
+            print 'Force : '+str(1e12*force)+' pN'
+            print 'Slope : '+str(slope)+' N/m'
+            try:
+                #FIXME all drivers should provide retract velocity, in SI or homogeneous units    
+                lrate=slope*self.current.curve.retract_velocity
+                print 'Loading rate (UNTESTED):'+str(lrate)
+            except:
+                print 'Loading rate : n/a'
+                lrate='n/a'
+            
+            if justone:
+                return
+            
+            #save accept/discard/repeat
+            print '\n\nDo you want to save these?'
+            if self.YNclick():
+                    
+                #Save file info
+                if self.autofile=='':
+                    self.autofile=raw_input('Results filename? (press enter for a random generated name)')
+                    if self.autofile=='':
+                        [osf,name]=tempfile.mkstemp(prefix='analysis-',suffix='.txt',text=True,dir=self.config['hookedir'])
+                        print 'saving in '+name
+                        self.autofile=name
+                        os.close(osf)
+                        f=open(self.autofile,'a+')
+                        f.write('Analysis started '+time.asctime()+'\n')
+                        f.write('----------------------------------------\n')
+                        f.write(' File ; WLC Cont. length (nm) ; Sigma WLC cl;  Per. Length (nm) ; Sigma pl; FJC Cont. length (nm) ; Sigma FJC cl ; Kuhn length (nm); Sigma kl ; Force (pN) ; Slope (N/m) ; (untested) Loading rate (pN/s)\n')
+                        f.close()
+        
+                if not os.path.exists(self.autofile):
+                    f=open(self.autofile,'a+')
+                    f.write('Analysis started '+time.asctime()+'\n')
+                    f.write('----------------------------------------\n')
+                    f.write(' File ; WLC Cont. length (nm) ; Sigma WLC cl;  Per. Length (nm) ; Sigma pl; FJC Cont. length (nm) ; Sigma FJC cl ; Kuhn length (nm); Sigma kl ; Force (pN) ; Slope (N/m) ; (untested) Loading rate (pN/s)\n')
+                    f.close()
+                
+                print 'Saving...'
+                savecountr+=1
+                f=open(self.autofile,'a+')
+                f.write(self.current.path)
+                f.write(' ; '+str(1e9*wlcparams[0])+' ; '+str(wlcfit_nm[0])+' ; '+str(1e9*wlcparams[1])+' ; '+str(wlcfit_nm[1])+' ; '+str(1e9*fjcparams[0])+' ; '+str(fjcfit_nm[0])+' ; '+str(1e9*fjcparams[1])+' ; '+str(fjcfit_nm[1])+' ; '+str(1e12*force)+' ; '+ str(slope)+' ; '+str(lrate)+'\n')
+                f.close()
+            else:
+                print '\nWould you like to try again on this same curve?'
+                if self.YNclick():
+                    continue
+            curveindex+=1
+
+        if not justone:
+            print 'You saved '+str(savecounter)+' out of '+str(len(self.current_list))+' total curves.'
+            
+            
+            
+    def YNclick(self):
+        print 'Click any point over y=0 for Yes, under it for No'
+        point=self._measure_N_points(N=1,whatset=1)[0]
+        value=point.absolute_coords[1]
+        if value<0:
+            return False
+        else:
+            return True
+
+
+
+